From 01cc2e688f8eb7679bbf9f6f194d49c68bd8828e Mon Sep 17 00:00:00 2001
From: Thomas Feuvrier <thomas.feuvrier@c-s.fr>
Date: Tue, 21 Mar 2006 08:01:06 +0000
Subject: [PATCH] nomsg
---
Code/BasicFilters/otbFrostImageFilter.txx | 206 ++++++++++++++++++++++
Code/BasicFilters/otbLeeImageFilter.h | 131 ++++++++++++++
Code/BasicFilters/otbLeeImageFilter.txx | 181 +++++++++++++++++++
3 files changed, 518 insertions(+)
create mode 100644 Code/BasicFilters/otbFrostImageFilter.txx
create mode 100644 Code/BasicFilters/otbLeeImageFilter.h
create mode 100644 Code/BasicFilters/otbLeeImageFilter.txx
diff --git a/Code/BasicFilters/otbFrostImageFilter.txx b/Code/BasicFilters/otbFrostImageFilter.txx
new file mode 100644
index 0000000000..a38ad1de60
--- /dev/null
+++ b/Code/BasicFilters/otbFrostImageFilter.txx
@@ -0,0 +1,206 @@
+/*=========================================================================
+
+ 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
diff --git a/Code/BasicFilters/otbLeeImageFilter.h b/Code/BasicFilters/otbLeeImageFilter.h
new file mode 100644
index 0000000000..252d9043ea
--- /dev/null
+++ b/Code/BasicFilters/otbLeeImageFilter.h
@@ -0,0 +1,131 @@
+/*=========================================================================
+
+ 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
diff --git a/Code/BasicFilters/otbLeeImageFilter.txx b/Code/BasicFilters/otbLeeImageFilter.txx
new file mode 100644
index 0000000000..ad8edc4bde
--- /dev/null
+++ b/Code/BasicFilters/otbLeeImageFilter.txx
@@ -0,0 +1,181 @@
+/*=========================================================================
+
+ 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
--
GitLab