nomsg

Code/BasicFilters/otbFrostImageFilter.txx

+/*=========================================================================
+
+  Programme :   OTB (ORFEO ToolBox)
+  Auteurs   :   CS - P.Imbo
+  Language  :   C++
+  Date      :   24 janvier 2006
+  Role      :   Filtre de débruitage de FROST sur une image 
+  $Id$ 
+
+=========================================================================*/
+#ifndef __otbFrostImageFilter_txx
+#define __otbFrostImageFilter_txx
+
+#include "otbFrostImageFilter.h"
+
+#include "itkDataObject.h"
+#include "itkExceptionObject.h"
+#include "itkConstNeighborhoodIterator.h"
+#include "itkNeighborhoodInnerProduct.h"
+#include "itkImageRegionIterator.h"
+#include "itkNeighborhoodAlgorithm.h"
+#include "itkZeroFluxNeumannBoundaryCondition.h"
+#include "itkOffset.h"
+#include "itkProgressReporter.h"
+#include <math.h>
+
+namespace otb
+{
+
+/**
+ *
+ */
+template <class TInputImage, class TOutputImage>
+FrostImageFilter<TInputImage, TOutputImage>::FrostImageFilter()
+{
+  m_Radius.Fill(1);
+}
+
+template <class TInputImage, class TOutputImage>
+void FrostImageFilter<TInputImage, TOutputImage>::GenerateInputRequestedRegion() throw (itk::InvalidRequestedRegionError)
+{
+  // call the superclass' implementation of this method
+  Superclass::GenerateInputRequestedRegion();
+  
+  // get pointers to the input and output
+  typename Superclass::InputImagePointer inputPtr   =  const_cast< TInputImage * >( this->GetInput() );
+  typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
+  
+  if ( !inputPtr || !outputPtr )
+    {
+    return;
+    }
+
+  // get a copy of the input requested region (should equal the output
+  // requested region)
+  typename TInputImage::RegionType inputRequestedRegion;
+  inputRequestedRegion = inputPtr->GetRequestedRegion();
+
+  // pad the input requested region by the operator radius
+  inputRequestedRegion.PadByRadius( m_Radius );
+
+  // crop the input requested region at the input's largest possible region
+  if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) )
+    {
+    inputPtr->SetRequestedRegion( inputRequestedRegion );
+    return;
+    }
+  else
+    {
+    // Couldn't crop the region (requested region is outside the largest
+    // possible region).  Throw an exception.
+
+    // store what we tried to request (prior to trying to crop)
+    inputPtr->SetRequestedRegion( inputRequestedRegion );
+    
+    // build an exception
+    itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
+    itk::OStringStream msg;
+    msg << static_cast<const char *>(this->GetNameOfClass())
+        << "::GenerateInputRequestedRegion()";
+    e.SetLocation(msg.str().c_str());
+    e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
+    e.SetDataObject(inputPtr);
+    throw e;
+    }
+}
+
+
+template< class TInputImage, class TOutputImage>
+void FrostImageFilter< TInputImage, TOutputImage>::ThreadedGenerateData(	
+			const 	OutputImageRegionType& 		outputRegionForThread,
+                       	int 	threadId
+				)
+{
+  unsigned int i;
+  itk::ZeroFluxNeumannBoundaryCondition<InputImageType> 	nbc;
+  itk::ConstNeighborhoodIterator<InputImageType> 		bit;
+  typename itk::ConstNeighborhoodIterator<InputImageType>::OffsetType	off;
+  itk::ImageRegionIterator<OutputImageType> 			it;
+  
+  
+  // Allocate output
+  typename OutputImageType::Pointer     output = this->GetOutput();
+  typename InputImageType::ConstPointer input  = this->GetInput();
+  
+  // Find the data-set boundary "faces"
+  typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>::FaceListType 		faceList;
+  typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>::FaceListType::iterator 	fit;
+
+  itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType> bC;
+  faceList = bC(input, outputRegionForThread, m_Radius);
+
+
+  // support progress methods/callbacks
+  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
+  
+//  InputRealType pixel;
+  InputRealType sum;
+  InputRealType sum2;
+
+  double 	Moyenne, Variance; 
+  double	Alpha;   	 
+  double    	NormFiltre;
+  double	FrostFiltre;
+  double	CoefFiltre;
+  double	dPixel;
+  
+    
+  // Process each of the boundary faces.  These are N-d regions which border
+  // the edge of the buffer.
+  for (fit=faceList.begin(); fit != faceList.end(); ++fit)
+    { 
+    bit = itk::ConstNeighborhoodIterator<InputImageType>(m_Radius, input, *fit);
+    unsigned int neighborhoodSize = bit.Size();
+    unsigned int CenterPos        = bit.GetCenterNeighborhoodIndex(); 
+    it = itk::ImageRegionIterator<OutputImageType>(output, *fit);
+    bit.OverrideBoundaryCondition(&nbc);
+    bit.GoToBegin();
+
+    while ( ! bit.IsAtEnd() )
+      {
+      sum  = itk::NumericTraits<InputRealType>::Zero;
+      sum2 = itk::NumericTraits<InputRealType>::Zero;
+      for (i = 0; i < neighborhoodSize; ++i)
+        {
+        	dPixel = static_cast<double>( bit.GetPixel(i) );
+		sum += dPixel;
+		sum2 += dPixel * dPixel;
+        }
+      Moyenne   = sum  / double(neighborhoodSize);
+      Variance  = sum2 / double(neighborhoodSize) - Moyenne * Moyenne ;
+     
+      Alpha = m_Deramp * Variance / (Moyenne * Moyenne) ;
+      
+      NormFiltre  = 0.0;
+      FrostFiltre = 0.0;
+ 
+      const double rad_x = static_cast<double>(m_Radius[0]);
+      const double rad_y = static_cast<double>(m_Radius[1]);
+                 
+      for (double x = -rad_x; x<= rad_x; x++)      
+      	{	
+       	for (double y = -rad_y; y <= rad_y; y++)
+		{
+		double Dist = double(sqrt(x*x+y*y));
+		off[0]= static_cast<int>(x);
+		off[1]= static_cast<int>(y);
+//		i = (unsigned int)((y+rad_y)*(2*rad_y+1)+(x+rad_x));
+		dPixel= static_cast<double>( bit.GetPixel(off));
+//		dPixel= static_cast<double>( bit.GetPixel(i));
+		CoefFiltre = Alpha * exp(-Alpha *Dist);
+		NormFiltre  = NormFiltre  + CoefFiltre;
+		FrostFiltre = FrostFiltre + (CoefFiltre * dPixel);
+		}
+        }	      
+      dPixel = (FrostFiltre/NormFiltre);      
+      if (finite(dPixel)==0){
+      	dPixel = 0.;
+      }
+                  
+      it.Set( static_cast<OutputPixelType>( dPixel ) );
+      
+      ++bit;
+      ++it;
+      progress.CompletedPixel();
+      
+      }
+    }
+}
+
+/**
+ * Standard "PrintSelf" method
+ */
+template <class TInputImage, class TOutput>
+void 
+FrostImageFilter<TInputImage, TOutput>::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+  Superclass::PrintSelf( os, indent );
+  os << indent << "Radius: " << m_Radius << std::endl;
+}
+
+
+} // end namespace otb
+
+
+#endif

Code/BasicFilters/otbLeeImageFilter.h

+/*=========================================================================
+
+  Programme :   OTB (ORFEO ToolBox)
+  Auteurs   :   CS - P.Imbo
+  Language  :   C++
+  Date      :   23 janvier 2006
+  Role      :  Filtre de débruitage de LEE sur une image 
+  $Id$ 
+
+=========================================================================*/
+#ifndef __otbLeeImageFilter_h
+#define __otbLeeImageFilter_h
+
+#include "itkImageToImageFilter.h"
+#include "itkImage.h"
+#include "itkNumericTraits.h"
+
+namespace otb
+{
+
+/** \class LeeImageFilter
+ * \brief Applique un filtre de débruitage sur une image.
+ *
+ * Ce filtre correspond au filtre Anti-Speckle de LEE :
+ *
+ *	R = E[I] + b(I-E[I]) avec 	b  = C²r / ( C²r + C²v )
+ *					Cv = 1 / sqrt(L)  avec L le nombre de vues
+ *					Cr = sqrt(Var(I)) / E[I] avec Var(I) = E[I²] - E[I]²
+ *
+ */
+
+template <class TInputImage, class TOutputImage>
+class LeeImageFilter :  public itk::ImageToImageFilter< TInputImage, TOutputImage >
+{
+public:
+  /** 	Extrait les dimensions aussi bien des images 
+  	d'entrée (Input) que de sortie (Output). */
+  itkStaticConstMacro(		InputImageDimension,
+  				unsigned int,
+                      		TInputImage::ImageDimension);
+  itkStaticConstMacro(		OutputImageDimension, 
+  				unsigned int,
+                      		TOutputImage::ImageDimension);
+
+  /** "typedef" pour simplifier la définition et la déclaration de variables. */
+  typedef TInputImage InputImageType;
+  /** "typedef" pour simplifier la définition et la déclaration de variables. */
+  typedef TOutputImage OutputImageType;
+
+  /** "typedef" pour les classes standards. */
+  typedef LeeImageFilter Self;
+  typedef itk::ImageToImageFilter< InputImageType, OutputImageType> Superclass;
+  typedef itk::SmartPointer<Self> Pointer;
+  typedef itk::SmartPointer<const Self>  ConstPointer;
+
+  /** Methode pour la gestion "object factory". */
+  itkNewMacro(Self);
+
+  /** Retourne le nom de la classe. */
+  itkTypeMacro(LeeImageFilter, ImageToImageFilter);
+  
+  /** Définition des images supportées. */
+  typedef typename InputImageType::PixelType InputPixelType;
+  typedef typename OutputImageType::PixelType OutputPixelType;
+  /** "typedef" définissant un réel dans cette classe. */
+  typedef typename itk::NumericTraits<InputPixelType>::RealType InputRealType;
+  
+  typedef typename InputImageType::RegionType InputImageRegionType;
+  typedef typename OutputImageType::RegionType OutputImageRegionType;
+
+  /** "typedef" définissant la taille d'une image. */
+  typedef typename InputImageType::SizeType SizeType;
+
+  /** Positionne le rayon définissant le voisinage utilisé pour le calcul du filtre. */
+  itkSetMacro(Radius, SizeType);
+
+  /** Récupère le rayon définissant le voisinage utilisé pour le calcul du filtre. */
+  itkGetConstReferenceMacro(Radius, SizeType);
+  
+  /** Positionne le nombre de vues utilisé pour le calcul du filtre. */
+  itkSetMacro(NbVues, double);
+  /** Récupère le nombre de vues (référencé constant) utilisé pour le calcul du filtre. */
+  itkGetConstReferenceMacro(NbVues, double);
+
+  /** LeeImageFilter a besoin d'une zone de traitement plus large en entrée qu'en sortie 
+   * pemettant une utilisation needs a larger input requested region than
+   * the output requested region.  As such, LeeImageFilter needs
+   * to provide an implementation for GenerateInputRequestedRegion()
+   * in order to inform the pipeline execution model.
+   *
+   * \sa ImageToImageFilter::GenerateInputRequestedRegion() */
+  virtual void GenerateInputRequestedRegion() throw(itk::InvalidRequestedRegionError);
+
+protected:
+  LeeImageFilter();
+  virtual ~LeeImageFilter() {};
+  void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+  /** LeeImageFilter peut etre implémentée pour un traitement de filtre multithreaded.
+   * Ainsi, cette implémentation fournit la méthode ThreadedGenerateData()
+   * qui est appelée pour chaque thread du process. Les données image sont allouées automatiquement 
+   * par la classe "mère" en appelant la méthode ThreadedGenerateData(). ThreadedGenerateData peut seulement 
+   * écrire la portion de l'image spécifiée par le paramètre "outputRegionForThread"
+   *
+   *	Filtre de LEE : 
+   *	R = E[I] + b(I-E[I]) avec 	b  = C²r / ( C²r + C²v )
+   *					Cv = 1 / sqrt(L)  avec L le nombre de vues
+   *					Cr = sqrt(Var(I)) / E[I] avec Var(I) = E[I²] - E[I]²
+   *
+   * \sa ImageToImageFilter::ThreadedGenerateData(),
+   *     ImageToImageFilter::GenerateData() */
+  void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
+                            int threadId );
+
+private:
+  LeeImageFilter(const Self&); //purposely not implemented
+  void operator=(const Self&); //purposely not implemented
+
+  /** Déclaration du rayon */
+  SizeType m_Radius;
+  /** Déclaration du nombre de vues du filtre */
+  double m_NbVues;
+};
+} // end namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbLeeImageFilter.txx"
+#endif
+
+  
+#endif

Code/BasicFilters/otbLeeImageFilter.txx

+/*=========================================================================
+
+  Programme :   OTB (ORFEO ToolBox)
+  Auteurs   :   CS - P.Imbo
+  Language  :   C++
+  Date      :   23 janvier 2006
+  Role      :  Filtre de débruitage de LEE sur une image 
+  $Id$ 
+
+=========================================================================*/
+#ifndef __otbLeeImageFilter_txx
+#define __otbLeeImageFilter_txx
+
+#include "otbLeeImageFilter.h"
+
+#include "itkDataObject.h"
+#include "itkExceptionObject.h"
+#include "itkConstNeighborhoodIterator.h"
+#include "itkNeighborhoodInnerProduct.h"
+#include "itkImageRegionIterator.h"
+#include "itkNeighborhoodAlgorithm.h"
+#include "itkZeroFluxNeumannBoundaryCondition.h"
+#include "itkOffset.h"
+#include "itkProgressReporter.h"
+
+namespace otb
+{
+
+/**
+ *
+ */
+template <class TInputImage, class TOutputImage>
+LeeImageFilter<TInputImage, TOutputImage>::LeeImageFilter()
+{
+  m_Radius.Fill(1);
+}
+
+template <class TInputImage, class TOutputImage>
+void LeeImageFilter<TInputImage, TOutputImage>::GenerateInputRequestedRegion() throw (itk::InvalidRequestedRegionError)
+{
+  // call the superclass' implementation of this method
+  Superclass::GenerateInputRequestedRegion();
+  
+  // get pointers to the input and output
+  typename Superclass::InputImagePointer inputPtr   =  const_cast< TInputImage * >( this->GetInput() );
+  typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
+  
+  if ( !inputPtr || !outputPtr )
+    {
+    return;
+    }
+
+  // get a copy of the input requested region (should equal the output
+  // requested region)
+  typename TInputImage::RegionType inputRequestedRegion;
+  inputRequestedRegion = inputPtr->GetRequestedRegion();
+
+  // pad the input requested region by the operator radius
+  inputRequestedRegion.PadByRadius( m_Radius );
+
+  // crop the input requested region at the input's largest possible region
+  if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) )
+    {
+    inputPtr->SetRequestedRegion( inputRequestedRegion );
+    return;
+    }
+  else
+    {
+    // Couldn't crop the region (requested region is outside the largest
+    // possible region).  Throw an exception.
+
+    // store what we tried to request (prior to trying to crop)
+    inputPtr->SetRequestedRegion( inputRequestedRegion );
+    
+    // build an exception
+    itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
+    itk::OStringStream msg;
+    msg << static_cast<const char *>(this->GetNameOfClass())
+        << "::GenerateInputRequestedRegion()";
+    e.SetLocation(msg.str().c_str());
+    e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
+    e.SetDataObject(inputPtr);
+    throw e;
+    }
+}
+
+
+template< class TInputImage, class TOutputImage>
+void LeeImageFilter< TInputImage, TOutputImage>::ThreadedGenerateData(	
+			const 	OutputImageRegionType& 		outputRegionForThread,
+                       	int 	threadId
+				)
+{
+  unsigned int i;
+  itk::ZeroFluxNeumannBoundaryCondition<InputImageType> nbc;
+
+  itk::ConstNeighborhoodIterator<InputImageType> bit;
+  itk::ImageRegionIterator<OutputImageType> it;
+  
+  // Allocate output
+  typename OutputImageType::Pointer     output = this->GetOutput();
+  typename InputImageType::ConstPointer input  = this->GetInput();
+  
+  // Find the data-set boundary "faces"
+  typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>::FaceListType 		faceList;
+  typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>::FaceListType::iterator 	fit;
+
+  itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType> bC;
+  faceList = bC(input, outputRegionForThread, m_Radius);
+
+
+  // support progress methods/callbacks
+  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
+  
+//  InputRealType pixel;
+  InputRealType sum;
+  InputRealType sum2;
+
+  double Cr2, Cv2, E_I, I, Var_I, dPixel;
+  
+  
+//  dPixel = this->getNbVues();
+  dPixel = m_NbVues;
+  //Calcul du rapport
+  Cv2 = 1./(sqrt(dPixel));
+  Cv2 = Cv2*Cv2;
+
+  // Process each of the boundary faces.  These are N-d regions which border
+  // the edge of the buffer.
+  for (fit=faceList.begin(); fit != faceList.end(); ++fit)
+    { 
+    bit = itk::ConstNeighborhoodIterator<InputImageType>(m_Radius, input, *fit);
+    unsigned int neighborhoodSize = bit.Size();
+    it = itk::ImageRegionIterator<OutputImageType>(output, *fit);
+    bit.OverrideBoundaryCondition(&nbc);
+    bit.GoToBegin();
+
+    while ( ! bit.IsAtEnd() )
+      {
+      sum = itk::NumericTraits<InputRealType>::Zero;
+      sum2 = itk::NumericTraits<InputRealType>::Zero;
+      //Parcours du voisinage
+      for (i = 0; i < neighborhoodSize; ++i)
+        {
+        	dPixel = static_cast<double>( bit.GetPixel(i) );
+		sum += dPixel;
+		sum2 += dPixel * dPixel;
+        }
+      
+      E_I   = sum / double(neighborhoodSize);
+      Var_I = sum2 / double(neighborhoodSize) - E_I*E_I ;
+      I = static_cast<double>( bit.GetCenterPixel() );
+      Cr2    = Var_I / (E_I*E_I);
+      
+      dPixel = E_I + ((I-E_I)*(Cr2))/(Cr2 + Cv2);
+      // get the mean value
+      it.Set( static_cast<OutputPixelType>( dPixel ) );
+      
+      ++bit;
+      ++it;
+      progress.CompletedPixel();
+      }
+    }
+}
+
+/**
+ * Standard "PrintSelf" method
+ */
+template <class TInputImage, class TOutput>
+void 
+LeeImageFilter<TInputImage, TOutput>::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+  Superclass::PrintSelf( os, indent );
+  os << indent << "Radius: " << m_Radius << std::endl;
+}
+
+
+} // end namespace otb
+
+
+#endif