From 2bd206e2c57b3a36f9773e6837fce79f17310215 Mon Sep 17 00:00:00 2001
From: Guillaume Borrut <guillaume.borrut@c-s.fr>
Date: Wed, 13 Aug 2008 08:24:25 +0000
Subject: [PATCH] Ajout des tests de syntheses de polarimetrie
---
...MultiChannelsPolarimetricSynthesisFilter.h | 194 +++++++++
...ltiChannelsPolarimetricSynthesisFilter.txx | 403 ++++++++++++++++++
2 files changed, 597 insertions(+)
create mode 100755 Code/SARPolarimetry/otbMultiChannelsPolarimetricSynthesisFilter.h
create mode 100755 Code/SARPolarimetry/otbMultiChannelsPolarimetricSynthesisFilter.txx
diff --git a/Code/SARPolarimetry/otbMultiChannelsPolarimetricSynthesisFilter.h b/Code/SARPolarimetry/otbMultiChannelsPolarimetricSynthesisFilter.h
new file mode 100755
index 0000000000..a8bda74e37
--- /dev/null
+++ b/Code/SARPolarimetry/otbMultiChannelsPolarimetricSynthesisFilter.h
@@ -0,0 +1,194 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef __otbMultiChannelsPolarimetricSynthesisFilter_h
+#define __otbMultiChannelsPolarimetricSynthesisFilter_h
+
+#include "itkInPlaceImageFilter.h"
+#include "otbPolarimetricFunctor.h"
+#include "itkFixedArray.h"
+#include <complex>
+
+namespace otb
+{
+
+/** \class MultiChannelsPolarimetricSynthesisFilter
+ * \brief Implements
+ *
+ * This class is parameterized over the type of the input image and
+ * the type of the output image. It is also parameterized by the
+ * operation to be applied, using a Functor style.
+ *
+ */
+
+
+template <class TInputImage, class TOutputImage,
+ class TFunction = Functor::PolarimetricFunctor<
+ typename TInputImage::InternalPixelType,
+ typename TInputImage::InternalPixelType,
+ typename TInputImage::InternalPixelType,
+ typename TInputImage::InternalPixelType,
+ typename TOutputImage::PixelType> >
+class ITK_EXPORT MultiChannelsPolarimetricSynthesisFilter : public itk::InPlaceImageFilter<TInputImage,TOutputImage>
+{
+public:
+ /** Standard class typedefs. */
+ typedef MultiChannelsPolarimetricSynthesisFilter Self;
+ typedef itk::InPlaceImageFilter<TInputImage,TOutputImage> Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Method for creation through the object factory. */
+ itkNewMacro(Self);
+
+ /** Run-time type information (and related methods). */
+ itkTypeMacro(MultiChannelsPolarimetricSynthesisFilter, InPlaceImageFilter);
+
+ /** Some typedefs. */
+ typedef TFunction FunctorType;
+ typedef TInputImage InputImageType;
+ typedef typename InputImageType::ConstPointer InputImagePointer;
+ typedef typename InputImageType::RegionType InputImageRegionType;
+ typedef typename InputImageType::PixelType InputImagePixelType;
+ typedef TOutputImage OutputImageType;
+ typedef typename OutputImageType::Pointer OutputImagePointer;
+ typedef typename OutputImageType::RegionType OutputImageRegionType;
+ typedef typename OutputImageType::PixelType OutputImagePixelType;
+ typedef typename std::complex <double> ComplexType;
+ typedef typename itk::FixedArray<ComplexType,2> ComplexArrayType;
+
+ /** Get the functor object. The functor is returned by reference.
+ * (Functors do not have to derive from itk::LightObject, so they do
+ * not necessarily have a reference count. So we cannot return a
+ * SmartPointer.) */
+ FunctorType& GetFunctor() { return m_Functor; };
+ const FunctorType& GetFunctor() const { return m_Functor; };
+
+
+ /** Set the functor object. This replaces the current Functor with a
+ * copy of the specified Functor. This allows the user to specify a
+ * functor that has ivars set differently than the default functor.
+ * This method requires an operator!=() be defined on the functor
+ * (or the compiler's default implementation of operator!=() being
+ * appropriate). */
+ void SetFunctor(const FunctorType& functor)
+ {
+ if (m_Functor != functor)
+ {
+ m_Functor = functor;
+ this->Modified();
+ }
+ }
+
+ void SetEi(ComplexArrayType ei)
+ {
+ m_Ei = ei;
+ this->GetFunctor().SetEi(ei);
+ this->Modified();
+ }
+
+ void SetEr(ComplexArrayType er)
+ {
+ m_Er = er;
+ this->GetFunctor().SetEr(er);
+ this->Modified();
+ }
+
+ /** Set/Get PsiI */
+ itkSetMacro(PsiI,double);
+ itkGetMacro(PsiI,double);
+ /** Set/Get TauI */
+ itkSetMacro(TauI,double);
+ itkGetMacro(TauI,double);
+ /** Set/Get PsiR */
+ itkSetMacro(PsiR,double);
+ itkGetMacro(PsiR,double);
+ /** Set/Get TauR */
+ itkSetMacro(TauR,double);
+ itkGetMacro(TauR,double);
+
+
+protected:
+ MultiChannelsPolarimetricSynthesisFilter();
+ virtual ~MultiChannelsPolarimetricSynthesisFilter(){};
+
+ /** MultiChannelsPolarimetricSynthesisFilter can produce an image which is a different
+ * resolution than its input image. As such, MultiChannelsPolarimetricSynthesisFilter
+ * needs to provide an implementation for
+ * GenerateOutputInformation() in order to inform the pipeline
+ * execution model. The original documentation of this method is
+ * below.
+ *
+ * \sa ProcessObject::GenerateOutputInformaton() */
+ virtual void GenerateOutputInformation();
+
+ virtual void BeforeThreadedGenerateData();
+
+ /** MultiChannelsPolarimetricSynthesisFilter can be implemented as a multithreaded filter.
+ * Therefore, this implementation provides a ThreadedGenerateData() routine
+ * which is called for each processing thread. The output image data is
+ * allocated automatically by the superclass prior to calling
+ * ThreadedGenerateData(). ThreadedGenerateData can only write to the
+ * portion of the output image specified by the parameter
+ * "outputRegionForThread"
+ *
+ * \sa ImageToImageFilter::ThreadedGenerateData(),
+ * ImageToImageFilter::GenerateData() */
+ void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
+ int threadId );
+
+ /** Computation of the electromagnetic fields Ei Er */
+ void ComputeElectromagneticFields();
+
+ /** Verify and force the inputs, if only 2 or 3 channels are present */
+ void VerifyAndForceInputs();
+
+ void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+private :
+ MultiChannelsPolarimetricSynthesisFilter(const Self&); //purposely not implemented
+
+ /** Psi Incident */
+ double m_PsiI;
+ /** Tau Incident */
+ double m_TauI;
+ /** Psi Relechi */
+ double m_PsiR;
+ /** Tau Relechi */
+ double m_TauR;
+
+ /** Champs Electromagnetic Incident */
+ ComplexArrayType m_Ei;
+ /** Champs Electromagnetic Reflechi */
+ ComplexArrayType m_Er;
+
+ /** Conversion coefficient Degre To Radian */
+ static const double m_DTOR=M_PI/180;
+
+ /** Channel indexes 0 1 2 3 on verra */
+
+ FunctorType m_Functor;
+
+};
+
+} // end namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbMultiChannelsPolarimetricSynthesisFilter.txx"
+#endif
+
+#endif
diff --git a/Code/SARPolarimetry/otbMultiChannelsPolarimetricSynthesisFilter.txx b/Code/SARPolarimetry/otbMultiChannelsPolarimetricSynthesisFilter.txx
new file mode 100755
index 0000000000..4bb4242456
--- /dev/null
+++ b/Code/SARPolarimetry/otbMultiChannelsPolarimetricSynthesisFilter.txx
@@ -0,0 +1,403 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef __otbMultiChannelsPolarimetricSynthesisFilter_txx
+#define __otbMultiChannelsPolarimetricSynthesisFilter_txx
+
+#include "otbMultiChannelsPolarimetricSynthesisFilter.h"
+#include "itkImageRegionIterator.h"
+#include "itkImageRegionConstIterator.h"
+#include "itkProgressReporter.h"
+#include <complex>
+
+namespace otb
+{
+
+/**
+ * Constructor
+ */
+template <class TInputImage, class TOutputImage, class TFunction >
+MultiChannelsPolarimetricSynthesisFilter<TInputImage,TOutputImage,TFunction>
+::MultiChannelsPolarimetricSynthesisFilter()
+{
+ this->SetNumberOfRequiredInputs( 1 );
+ //this->InPlaceOff();
+}
+
+ /** PolarimetricSynthesisFilter
+ *
+ *
+ * \sa ProcessObject::GenerateOutputInformaton()
+ */
+template <class TInputImage, class TOutputImage, class TFunction>
+void
+MultiChannelsPolarimetricSynthesisFilter<TInputImage,TOutputImage,TFunction>
+::GenerateOutputInformation()
+{
+ // do not call the superclass' implementation of this method since
+ // this filter allows the input the output to be of different dimensions
+
+ // get pointers to the input and output
+ typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
+ typename Superclass::InputImageConstPointer inputPtr = this->GetInput();
+
+ if ( !outputPtr || !inputPtr)
+ {
+ return;
+ }
+
+ // Set the output image largest possible region. Use a RegionCopier
+ // so that the input and output images can be different dimensions.
+ OutputImageRegionType outputLargestPossibleRegion;
+ this->CallCopyInputRegionToOutputRegion(outputLargestPossibleRegion,
+ inputPtr->GetLargestPossibleRegion());
+ outputPtr->SetLargestPossibleRegion( outputLargestPossibleRegion );
+
+ // Set the output spacing and origin
+ const itk::ImageBase<Superclass::InputImageDimension> *phyData;
+
+ phyData
+ = dynamic_cast<const itk::ImageBase<Superclass::InputImageDimension>*>(this->GetInput());
+
+ if (phyData)
+ {
+ // Copy what we can from the image from spacing and origin of the input
+ // This logic needs to be augmented with logic that select which
+ // dimensions to copy
+ unsigned int i, j;
+ const typename InputImageType::SpacingType&
+ inputSpacing = inputPtr->GetSpacing();
+ const typename InputImageType::PointType&
+ inputOrigin = inputPtr->GetOrigin();
+ const typename InputImageType::DirectionType&
+ inputDirection = inputPtr->GetDirection();
+
+ typename OutputImageType::SpacingType outputSpacing;
+ typename OutputImageType::PointType outputOrigin;
+ typename OutputImageType::DirectionType outputDirection;
+
+ // copy the input to the output and fill the rest of the
+ // output with zeros.
+ for (i=0; i < Superclass::InputImageDimension; ++i)
+ {
+ outputSpacing[i] = inputSpacing[i];
+ outputOrigin[i] = inputOrigin[i];
+ for (j=0; j < Superclass::OutputImageDimension; j++)
+ {
+ if (j < Superclass::InputImageDimension)
+ {
+ outputDirection[j][i] = inputDirection[j][i];
+ }
+ else
+ {
+ outputDirection[j][i] = 0.0;
+ }
+ }
+ }
+ for (; i < Superclass::OutputImageDimension; ++i)
+ {
+ outputSpacing[i] = 1.0;
+ outputOrigin[i] = 0.0;
+ for (j=0; j < Superclass::OutputImageDimension; j++)
+ {
+ if (j == i)
+ {
+ outputDirection[j][i] = 1.0;
+ }
+ else
+ {
+ outputDirection[j][i] = 0.0;
+ }
+ }
+ }
+
+ // set the spacing and origin
+ outputPtr->SetSpacing( outputSpacing );
+ outputPtr->SetOrigin( outputOrigin );
+ outputPtr->SetDirection( outputDirection );
+
+ }
+ else
+ {
+ // pointer could not be cast back down
+ itkExceptionMacro(<< "otb::MultiChannelRAndNIRVegetationIndexImageFilter::GenerateOutputInformation "
+ << "cannot cast input to "
+ << typeid(itk::ImageBase<Superclass::InputImageDimension>*).name() );
+ }
+}
+
+
+/**
+ * ThreadedGenerateData Performs the pixel-wise addition
+ */
+template <class TInputImage, class TOutputImage, class TFunction >
+void
+MultiChannelsPolarimetricSynthesisFilter<TInputImage,TOutputImage,TFunction>
+::ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread,
+ int threadId)
+{
+ InputImagePointer inputPtr = this->GetInput();
+ OutputImagePointer outputPtr = this->GetOutput(0);
+
+ // Define the portion of the input to walk for this thread, using
+ // the CallCopyOutputRegionToInputRegion method allows for the input
+ // and output images to be different dimensions
+ InputImageRegionType inputRegionForThread;
+ this->CallCopyOutputRegionToInputRegion(inputRegionForThread, outputRegionForThread);
+
+ // Define the iterators
+ itk::ImageRegionConstIterator<TInputImage> inputIt(inputPtr, inputRegionForThread);
+ itk::ImageRegionIterator<TOutputImage> outputIt(outputPtr, outputRegionForThread);
+
+ itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
+
+ inputIt.GoToBegin();
+ outputIt.GoToBegin();
+
+ // Computation with 4 channels
+ if ( inputPtr->GetNumberOfComponentsPerPixel()==4 )
+ {
+ while( !inputIt.IsAtEnd() )
+ {
+ outputIt.Set( m_Functor( inputIt.Get()[0], inputIt.Get()[1], inputIt.Get()[2], inputIt.Get()[3] ) );
+ ++inputIt;
+ ++outputIt;
+ progress.CompletedPixel(); // potential exception thrown here
+ }
+ }
+
+}
+
+/**
+ * Computation of the electromagnetic fields Ei Er
+ */
+template <class TInputImage, class TOutputImage, class TFunction >
+void
+MultiChannelsPolarimetricSynthesisFilter<TInputImage,TOutputImage,TFunction>
+::ComputeElectromagneticFields()
+{
+ ComplexType Ei0,Ei1,Er0,Er1;
+ ComplexArrayType AEi, AEr;
+
+ Ei0.real() = vcl_cos(m_DTOR*m_PsiI)*vcl_cos(m_DTOR*m_TauI);
+ Ei0.imag() = -vcl_sin(m_DTOR*m_PsiI)*vcl_sin(m_DTOR*m_TauI);
+
+ Ei1.real() = vcl_sin(m_DTOR*m_PsiI)*vcl_cos(m_DTOR*m_TauI);
+ Ei1.imag() = vcl_cos(m_DTOR*m_PsiI)*vcl_sin(m_DTOR*m_TauI);
+
+ Er0.real() = vcl_cos(m_DTOR*m_PsiR)*vcl_cos(m_DTOR*m_TauR);
+ Er0.imag() = -vcl_sin(m_DTOR*m_PsiR)*vcl_sin(m_DTOR*m_TauR);
+
+ Er1.real() = vcl_sin(m_DTOR*m_PsiR)*vcl_cos(m_DTOR*m_TauR);
+ Er1.imag() = vcl_cos(m_DTOR*m_PsiR)*vcl_sin(m_DTOR*m_TauR);
+
+ AEi[0]=Ei0;
+ AEi[1]=Ei1;
+ AEr[0]=Er0;
+ AEr[1]=Er1;
+
+ this->SetEi(AEi);
+ this->SetEr(AEr);
+
+ std::cout<<"PsiI: "<<m_PsiI<<std::endl;
+ std::cout<<"TauI: "<<m_TauI<<std::endl;
+ std::cout<<"PsiR: "<<m_PsiR<<std::endl;
+ std::cout<<"TauR: "<<m_TauR<<std::endl;
+
+ std::cout<<"Ei0 im: "<<m_Ei[0].imag()<<std::endl;
+ std::cout<<"Ei0 re: "<<m_Ei[0].real()<<std::endl;
+ std::cout<<"Ei1 im: "<<m_Ei[1].imag()<<std::endl;
+ std::cout<<"Ei1 re: "<<m_Ei[1].real()<<std::endl;
+
+ std::cout<<"Er0 im: "<<m_Er[0].imag()<<std::endl;
+ std::cout<<"Er0 re: "<<m_Er[0].real()<<std::endl;
+ std::cout<<"Er1 im: "<<m_Er[1].imag()<<std::endl;
+ std::cout<<"Er1 re: "<<m_Er[1].real()<<std::endl;
+
+ std::cout<<"DTOR: "<<m_DTOR<<std::endl;
+}
+
+/**
+ * Verify and force the inputs, if only 2 or 3 channels are present
+ */
+template <class TInputImage, class TOutputImage, class TFunction >
+void
+MultiChannelsPolarimetricSynthesisFilter<TInputImage,TOutputImage,TFunction>
+::VerifyAndForceInputs()
+{
+/* InputPixelType pix;
+ pix.imag()=0;
+ pix.real()=0;
+
+ // With 3 channels : HH VH VV
+ if ( ( this->GetInput(0)!=0 && this->GetInput(1)==0 )&&
+ ( this->GetInput(2)!=0 && this->GetInput(3)!=0 ) )
+ {
+ typename HVInputImageType::ConstPointer hvImage = dynamic_cast<const HVInputImageType*>((itk::ProcessObject::GetInput(2)));;
+ this->SetInputHV(hvImage);
+ std::cout<<"Case 3 channels !!"<<std::endl;
+ }
+ else
+ // With 3 channels : HH HV VV
+ if ( ( this->GetInput(0)!=0 && this->GetInput(1)!=0 )&&
+ ( this->GetInput(2)==0 && this->GetInput(3)!=0 ) )
+ {
+ //this->SetInputVH(static_cast<typename HVInputImageType::Pointer>(this->GetInput(1)));
+ typename VHInputImageType::ConstPointer vhImage = dynamic_cast<const VHInputImageType*>((itk::ProcessObject::GetInput(1)));;
+ this->SetInputVH(vhImage);
+ std::cout<<"Case 3 channels !!"<<std::endl;
+ }
+ else
+ // Only VH and VV are present
+ if ( ( this->GetInput(0)==0 && this->GetInput(1)==0 ) &&
+ ( this->GetInput(2)!=0 && this->GetInput(3)!=0 ) )
+ {
+
+ std::cout<<"Case VH VV present !!"<<std::endl;
+ // Forcing HH and HV to zero
+ typename HHInputImageType::Pointer inputHH = TInputImageHH::New();
+ typename HVInputImageType::Pointer inputHV = TInputImageHV::New();
+
+ typename VHInputImageType::IndexType start;
+ start[0]=0;
+ start[1]=0;
+
+ typename VHInputImageType::SizeType size;
+ typename VHInputImageType::ConstPointer vhImage = dynamic_cast<const VHInputImageType*>((itk::ProcessObject::GetInput(2)));
+ size = vhImage->GetLargestPossibleRegion().GetSize();
+
+ typename VHInputImageType::RegionType region;
+ region.SetSize(size);
+ region.SetIndex(start);
+ inputHH->SetRegions(region);
+ inputHH->Allocate();
+
+ inputHH->FillBuffer(pix);
+
+ inputHV->SetRegions(region);
+ inputHV->Allocate();
+ inputHV->FillBuffer(pix);
+
+ this->SetInputHH(inputHH);
+ this->SetInputHV(inputHV);
+
+ // Forcing TauI=0 PsiI=90
+ this->SetTauI(0);
+ this->SetPsiI(90);
+ }
+ else
+ // Only HH and HV are present
+ if ( ( this->GetInput(0)!=0 && this->GetInput(1)!=0 ) &&
+ ( this->GetInput(2)==0 && this->GetInput(3)==0 ) )
+ {
+ std::cout<<"Case HH HV present !!"<<std::endl;
+ // Forcing HH and HV to zero
+ typename VVInputImageType::Pointer inputVV = TInputImageVV::New();
+ typename VHInputImageType::Pointer inputVH = TInputImageVH::New();
+ typename VHInputImageType::IndexType start;
+ start[0]=0;
+ start[1]=0;
+ typename VHInputImageType::SizeType size;
+ typename VHInputImageType::ConstPointer vhImage = dynamic_cast<const VHInputImageType*>((itk::ProcessObject::GetInput(1)));
+
+ size = vhImage->GetLargestPossibleRegion().GetSize();
+
+ typename VHInputImageType::RegionType region;
+ region.SetSize(size);
+ region.SetIndex(start);
+ inputVV->SetRegions(region);
+ inputVV->Allocate();
+ inputVV->FillBuffer(pix);
+
+ inputVH->SetRegions(region);
+ inputVH->Allocate();
+ inputVH->FillBuffer(pix);
+
+ this->SetInputVV(inputVV);
+ this->SetInputVH(inputVH);
+
+ // Forcing TauI=0 PsiI=0
+ this->SetTauI(0);
+ this->SetPsiI(0);
+ }
+ else
+ // Only HH and VV are present
+ if ( ( this->GetInput(0)!=0 && this->GetInput(1)==0 )&&
+ ( this->GetInput(2)==0 && this->GetInput(3)!=0 ) )
+ {
+ std::cout<<"Case HH VV present !!"<<std::endl;
+ itkExceptionMacro("Only the HH and VV channels are available : Polarimetric synthesis is impossible !");
+ return;
+ }
+ std::cout<<"Fin VerifyAndForceInputs !!"<<std::endl;
+
+ */
+}
+
+/**
+ *
+ */
+template <class TInputImage, class TOutputImage, class TFunction >
+void
+MultiChannelsPolarimetricSynthesisFilter<TInputImage,TOutputImage,TFunction>
+::BeforeThreadedGenerateData()
+{
+ try{
+ // First Part. Verify and force the inputs
+ VerifyAndForceInputs();
+
+ /* std::cout<<"image 1 "<<this->GetInput(0)<<std::endl;
+ std::cout<<"image 2 "<<this->GetInput(1)<<std::endl;
+ std::cout<<"image 3 "<<this->GetInput(2)<<std::endl;
+ std::cout<<"image 4 "<<this->GetInput(3)<<std::endl;
+ */
+ // Second Part. Estimation of the incident field Ei and the reflected field Er
+ ComputeElectromagneticFields();
+ }
+ catch( itk::ExceptionObject & err )
+ {
+ std::cout << "Exception itk::ExceptionObject levee !" << std::endl;
+ std::cout << "Exception in PolarimetricSynthesisFilter class :"<< err << std::endl;
+ return;
+ }
+ catch (...)
+ {
+ // itkGenericExceptionMacro( <<"Unknown exception in PolarimetricSynthesisFilter class (catch(...)");
+ std::cout << "Exception levee inconnue !" << std::endl;
+ return;
+ }
+
+}
+
+/**
+ * Printself
+ */
+template <class TInputImage, class TOutputImage, class TFunction >
+void
+MultiChannelsPolarimetricSynthesisFilter<TInputImage,TOutputImage,TFunction>
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+ this->Superclass::PrintSelf(os,indent);
+ os << indent << "Psi I: "<<m_PsiI<<std::endl;
+ os << indent << "Tau I: "<<m_TauI<<std::endl;
+ os << indent << "Psi R: "<<m_PsiR<<std::endl;
+ os << indent << "Tau R: "<<m_TauR<<std::endl;
+}
+
+} // end namespace otb
+
+#endif
--
GitLab