From 9496b6a12f13f15dc858665ae41b5cbfc7497fc9 Mon Sep 17 00:00:00 2001 From: Cyrille Valladeau <cyrille.valladeau@c-s.fr> Date: Tue, 15 Jan 2008 16:04:59 +0000 Subject: [PATCH] Correction sur les interpolateurs. --- .../otbProlateInterpolateImageFunction.h | 54 +++++++-------- .../otbProlateInterpolateImageFunction.txx | 69 +++---------------- .../otbGenericInterpolateImageFunction.h | 6 ++ .../otbGenericInterpolateImageFunction.txx | 38 +++++++--- .../otbProlateInterpolateImageFunction.cxx | 29 ++------ 5 files changed, 80 insertions(+), 116 deletions(-) diff --git a/Code/BasicFilters/otbProlateInterpolateImageFunction.h b/Code/BasicFilters/otbProlateInterpolateImageFunction.h index 05aae130df..a91918ca5d 100644 --- a/Code/BasicFilters/otbProlateInterpolateImageFunction.h +++ b/Code/BasicFilters/otbProlateInterpolateImageFunction.h @@ -35,27 +35,39 @@ class ProlateFunction { public: typedef typename std::vector<double> VectorType; - void SetRadius(unsigned int rad){ m_Radius = rad; }; - unsigned int GetRadius() const { return m_Radius; }; + // Accessors definitions - void SetProfil(VectorType vect){ m_Profil = vect; }; - VectorType GetProfil(){ return m_Profil; }; + void SetRadius(unsigned int rad){ m_Radius = rad; }; + unsigned int GetRadius() const { return m_Radius; }; + unsigned int GetOriginalProfilSize() const { return m_OriginalProfilSize; }; + VectorType GetOriginalProfil() const { return m_OriginalProfil;}; inline TOutput operator()( const TInput & A ) const - { std::cout<<"~~~~~Prolate : Function"<<std::endl; -std::cout<<"~~~~~Prolate : Function : "<<A<<" -> "<<static_cast<TOutput>(m_Profil[static_cast<unsigned int>(A)])<<std::endl; - return (static_cast<TOutput>(m_Profil[static_cast<unsigned int>(A)])); -std::cout<<"~~~~~Prolate : Function : FIN"<<std::endl; + { + TOutput val = 0.; + if (m_Radius != 0) + { + unsigned int ival = static_cast<unsigned int>(m_OriginalProfilSize*static_cast<double>(vcl_abs(A))/static_cast<double>(m_Radius)); + val = m_OriginalProfil[ival]; + } + else + { + val = 1.; + } + return val; } private: - /** Store the resampled profil.*/ - VectorType m_Profil; /** Useless, only to be compatible with the GenericInterpolateImage. */ unsigned int m_Radius; -}; + /** Length of the original profil. */ + static const unsigned int m_OriginalProfilSize; + /** Original prolate profil */ + static const double m_OriginalProfil[721]; +}; + }//namespace Function /** \class ProlateInterpolateImageFunction @@ -105,16 +117,11 @@ public GenericInterpolateImageFunction< TInputImage, typedef typename Superclass::IteratorType IteratorType; typedef typename Superclass::ContinuousIndexType ContinuousIndexType; typedef typename std::vector<double> VectorType; + + unsigned int GetOriginalProfilSize() const { return this->GetFunction().GetOriginalProfilSize; }; + VectorType GetOriginalProfil() const { return this->GetFunction().GetOriginalProfil();}; - /** Compute a resampled profil according to the window size.*/ - void ComputeResampledProlateProfil(); - - void SetRadius(unsigned int rad); - - /** Get the resampled profil.*/ - VectorType GetResampledProfil(){ return m_ResampledProfil; }; - - + protected: ProlateInterpolateImageFunction(); ~ProlateInterpolateImageFunction(); @@ -123,13 +130,6 @@ public GenericInterpolateImageFunction< TInputImage, private: ProlateInterpolateImageFunction(const Self&); //purposely not implemented void operator=(const Self&); //purposely not implemented - - /** Length of the original profil. */ - static const unsigned int m_OriginalProfilSize; - /** Original prolate profil */ - static const double m_OriginalProfil[721]; - /** Store the resampled profil.*/ - VectorType m_ResampledProfil; }; } // end namespace itk diff --git a/Code/BasicFilters/otbProlateInterpolateImageFunction.txx b/Code/BasicFilters/otbProlateInterpolateImageFunction.txx index f961387ebb..bf4f6cb756 100644 --- a/Code/BasicFilters/otbProlateInterpolateImageFunction.txx +++ b/Code/BasicFilters/otbProlateInterpolateImageFunction.txx @@ -20,14 +20,17 @@ PURPOSE. See the above copyright notices for more information. namespace otb { - template<class TInputImage, class TBoundaryCondition, class TCoordRep, class TInputInterpolator, class TOutputInterpolator> +namespace Function +{ + template<class TInput, class TOutput> const unsigned int - ProlateInterpolateImageFunction<TInputImage, TBoundaryCondition, TCoordRep, TInputInterpolator, TOutputInterpolator> + ProlateFunction<TInput, TOutput> ::m_OriginalProfilSize = 721; - template<class TInputImage, class TBoundaryCondition, class TCoordRep, class TInputInterpolator, class TOutputInterpolator> - const double - ProlateInterpolateImageFunction<TInputImage, TBoundaryCondition, TCoordRep, TInputInterpolator, TOutputInterpolator> + + template<class TInput, class TOutput> + const double + ProlateFunction<TInput, TOutput> ::m_OriginalProfil[721] = { 0.00125, 0.00124999, 0.00124998, 0.00124997, 0.00124995, 0.00124992, 0.00124989, 0.00124985, 0.0012498, 0.00124975, 0.00124969, 0.00124962, 0.00124955, 0.00124948, 0.00124939, 0.0012493, 0.00124921, 0.00124911, 0.001249, 0.00124888, 0.00124876, 0.00124864, 0.0012485, 0.00124837, @@ -91,13 +94,15 @@ namespace otb 0.000167415, 0.000165465, 0.000163517, 0.000161572, 0.000159629, 0.000157689, 0.000155752, 0.000153816, 0.000151884, 0.000149954, 0.000148026, 0.000146101, 0.000144179 }; +} /** Constructor */ template<class TInputImage, class TBoundaryCondition, class TCoordRep, class TInputInterpolator, class TOutputInterpolator> ProlateInterpolateImageFunction<TInputImage, TBoundaryCondition, TCoordRep, TInputInterpolator, TOutputInterpolator> ::ProlateInterpolateImageFunction() { - VectorType m_ResampledProfil(1, 0.); + //VectorType m_ResampledProfil(1, 0.); + this->SetNormalizeWeight(true); } /** Destructor */ @@ -107,64 +112,12 @@ ProlateInterpolateImageFunction<TInputImage, TBoundaryCondition, TCoordRep, TInp { } -// Overload method to add the construction of resampledprfil -template<class TInputImage, class TBoundaryCondition, class TCoordRep, class TInputInterpolator, class TOutputInterpolator> -void -ProlateInterpolateImageFunction<TInputImage, TBoundaryCondition, TCoordRep, TInputInterpolator, TOutputInterpolator> -::SetRadius(unsigned int rad) -{std::cout<<"*****Prolate : SetRadius"<<std::endl; - Superclass::SetRadius(rad); - VectorType temp(2*rad+1, 0.); - m_ResampledProfil = temp; - this->ComputeResampledProlateProfil(); - this->Modified(); -std::cout<<"*****Prolate : SetRadius : Fin"<<std::endl; -} - - -// Resamplation profils calculation -template<class TInputImage, class TBoundaryCondition, class TCoordRep, class TInputInterpolator, class TOutputInterpolator> -void -ProlateInterpolateImageFunction<TInputImage, TBoundaryCondition, TCoordRep, TInputInterpolator, TOutputInterpolator> -::ComputeResampledProlateProfil() -{ std::cout<<"*****Prolate : ComputeResampledProlateProfil"<<std::endl; - unsigned int ival; - double dval; - /* Initialisation sur les lignes : */ - /* Determine les profils en largeur et hauteur */ - double sum = 0.; - double mean = 0.; - for (unsigned int i=0; i<=this->GetRadius(); i++) - { - ival = static_cast<unsigned int>(m_OriginalProfilSize*static_cast<double>(i)/static_cast<double>(this->GetRadius()+1)); - dval = m_OriginalProfil[ival]; - m_ResampledProfil[this->GetRadius()+i] = dval; - m_ResampledProfil[this->GetRadius()-i] = dval; - sum += 2*dval; - } - sum -= m_ResampledProfil[this->GetRadius()];// the middled pixel was compted twice - mean = sum/(2*static_cast<double>(this->GetRadius())+1); - std::cout<<"*****Prolate : ComputeResampledProlateProfil : mean : "<<mean<<std::endl; - for (unsigned int i=0; i<m_ResampledProfil.size(); i++) - {std::cout<<"*****Prolate : ComputeResampledProlateProfil : value : "<<m_ResampledProfil[i]<<" -> "<<m_ResampledProfil[i]/mean<<std::endl; - m_ResampledProfil[i] = m_ResampledProfil[i]/mean; - } - this->GetFunction().SetProfil(m_ResampledProfil); -std::cout<<"*****Prolate : ComputeResampledProlateProfil : Fin"<<std::endl; -} - template<class TInputImage, class TBoundaryCondition, class TCoordRep, class TInputInterpolator, class TOutputInterpolator> void ProlateInterpolateImageFunction<TInputImage, TBoundaryCondition, TCoordRep, TInputInterpolator, TOutputInterpolator> ::PrintSelf(std::ostream& os, itk::Indent indent) const { Superclass::PrintSelf( os, indent ); - os << indent << "Resample prolate profil : " << std::endl; - for (unsigned int i=0; i<m_ResampledProfil.size(); i++) - { - os << indent << i+1 << " : " << m_ResampledProfil[i] << std::endl; - } - } }//namespace otb diff --git a/Code/Common/otbGenericInterpolateImageFunction.h b/Code/Common/otbGenericInterpolateImageFunction.h index a91d948027..2bea862425 100644 --- a/Code/Common/otbGenericInterpolateImageFunction.h +++ b/Code/Common/otbGenericInterpolateImageFunction.h @@ -104,6 +104,10 @@ public itk::InterpolateImageFunction<TInputImage,TCoordRep> /** Fill the weight offset table*/ void FillWeightOffsetTable() const; + /** Weights normalization accessors*/ + itkSetMacro(NormalizeWeight, bool); + itkGetMacro(NormalizeWeight, bool); + protected: GenericInterpolateImageFunction(); ~GenericInterpolateImageFunction(); @@ -138,6 +142,8 @@ public itk::InterpolateImageFunction<TInputImage,TCoordRep> mutable unsigned int **m_WeightOffsetTable; /** True if internal statistics have been generated */ mutable bool m_TablesHaveBeenGenerated; + /** Weights normalization */ + bool m_NormalizeWeight; }; } // end namespace itk diff --git a/Code/Common/otbGenericInterpolateImageFunction.txx b/Code/Common/otbGenericInterpolateImageFunction.txx index 06d9e89728..974c728d59 100644 --- a/Code/Common/otbGenericInterpolateImageFunction.txx +++ b/Code/Common/otbGenericInterpolateImageFunction.txx @@ -31,6 +31,7 @@ GenericInterpolateImageFunction<TInputImage, TFunction, TBoundaryCondition, TCoo m_OffsetTable = NULL; m_WeightOffsetTable = NULL; m_TablesHaveBeenGenerated=false; + m_NormalizeWeight = false; } /** Destructor */ @@ -261,17 +262,38 @@ GenericInterpolateImageFunction<TInputImage, TFunction, TBoundaryCondition, TCoo // i is the relative offset in dimension dim. for( unsigned int i = 0; i < m_WindowSize; i++) { - // Increment the offset, taking it through the range - // (dist + rad - 1, ..., dist - rad), i.e. all x - // such that vcl_abs(x) <= rad - x -= 1.0; - - // Compute the weight for this m - xWeight[dim][i] = m_Function(x); + // Increment the offset, taking it through the range + // (dist + rad - 1, ..., dist - rad), i.e. all x + // such that vcl_abs(x) <= rad + x -= 1.0; + // Compute the weight for this m + xWeight[dim][i] = m_Function(x); } } } - + if (m_NormalizeWeight == true) + { + double sum = 0.; + for( unsigned int dim = 0; dim < ImageDimension; dim++ ) + { + for( unsigned int i = 0; i < m_WindowSize; i++) + { + // Compute the weight for this m + sum += xWeight[dim][i]; + } + } + for( unsigned int dim = 0; dim < ImageDimension; dim++ ) + { + for( unsigned int i = 0; i < m_WindowSize; i++) + { + // Compute the weight for this m + xWeight[dim][i] = xWeight[dim][i]/sum; + } + } + } + + + // Iterate over the neighborhood, taking the correct set // of weights in each dimension double xPixelValue = 0.0; diff --git a/Testing/Code/BasicFilters/otbProlateInterpolateImageFunction.cxx b/Testing/Code/BasicFilters/otbProlateInterpolateImageFunction.cxx index 1aa3495ea6..184876775c 100644 --- a/Testing/Code/BasicFilters/otbProlateInterpolateImageFunction.cxx +++ b/Testing/Code/BasicFilters/otbProlateInterpolateImageFunction.cxx @@ -35,51 +35,34 @@ int otbProlateInterpolateImageFunction(int argc, char * argv[]) unsigned int i = 4; std::vector<ContinuousIndexType> indicesList; - std::cout<<"Prolate : Test : while varie jusqua : "<<static_cast<unsigned int>(argc)<<std::endl; while(i<static_cast<unsigned int>(argc) && (i+1)<static_cast<unsigned int>(argc)) - { std::cout<<"Prolate : Test : while i vaut : "<<i<<" -> "<<atof(argv[i])<<" , "<<atof(argv[i+1])<<std::endl; + { ContinuousIndexType idx; - std::cout<<"Prolate : Test : while"<<std::endl; idx[0]=atof(argv[i]); - idx[1]=atof(argv[i+1]); - + idx[1]=atof(argv[i+1]); indicesList.push_back(idx); i+=2; } - std::cout<<"Prolate : Test : New()"<<std::endl; + // Instantiating object InterpolatorType::Pointer prolate = InterpolatorType::New(); - std::cout<<"Prolate : Test : Reader::New()"<<std::endl; ReaderType::Pointer reader = ReaderType::New(); - std::cout<<"Prolate : Test : Reader::SetFileName : "<<infname<<std::endl; reader->SetFileName(infname); - std::cout<<"Prolate : Test : Update()"<<std::endl; reader->Update(); - std::cout<<"Prolate : Test : prolate::Setinput()"<<std::endl; prolate->SetInputImage(reader->GetOutput()); - std::cout<<"Prolate : Test : prolate::SetRadius() : "<<atoi(argv[3])<<std::endl; prolate->SetRadius(atoi(argv[3])); - std::cout<<"Prolate : Test : output name : "<<outfname<<std::endl; - /************************ - typedef InterpolatorType::VectorType VectType; - VectType tutu=prolate->GetResampledProfil(); - for (unsigned int i=0; i<tutu.size(); i++) - { - std::cout << i+1 << " : " << tutu[i] << std::endl; - } - ************************/ + std::ofstream file; file.open(outfname); - std::cout<<"Prolate : Test : for : "<<outfname<<std::endl; + unsigned mmm = 1; for(std::vector<ContinuousIndexType>::iterator it = indicesList.begin();it!=indicesList.end();++it) - {std::cout<<"Prolate : Test : for indice : "<<mmm<<std::endl; + { file<<(*it)<<" -> "<<prolate->EvaluateAtContinuousIndex((*it))<<std::endl; mmm++; } - std::cout<<"Prolate : Test : close"<<std::endl; file.close(); return EXIT_SUCCESS; -- GitLab