ADD: MuellerToCircularPolarisationImageFilter class (C Tison source)

Code/SARPolarimetry/otbMuellerToCircularPolarisationImageFilter.h

+/*=========================================================================
+
+  Program:   OTB
+  Module:    $otbMuellerToCircularPolarisationImageFilter.h$
+  Language:  C++
+  Date:      $Date: 24/10/2006 $
+  Version:   $Revision: 1.0 $
+
+ 
+
+=========================================================================*/
+#ifndef __MuellerToCircularPolarisationImageFilter_h
+#define __MuellerToCircularPolarisationImageFilter_h
+
+#include "itkImageToImageFilter.h"
+#include "itkImage.h"
+#include "itkNumericTraits.h"
+#include "itkImageRegionIterator.h"
+#include "itkImageRegionConstIterator.h"
+#include "itkImageRegionConstIteratorWithIndex.h"
+#include "itkVector.h"
+
+#include "otbImage.h"
+
+namespace otb
+ {
+/** \class otbMuellerToCircularPolarisationImageFilter
+ * \brief Compute the Circular Polarisation image (3 channels : LL, RR and LR) from the Mueller Matrix image (10 channels)
+ */
+template <class TPixel>
+class  MuellerToCircularPolarisationImageFilter :
+   public itk::ImageToImageFilter< otb::Image<itk::Vector<TPixel,10>,2>, otb::Image<itk::Vector<TPixel,3>,2> >
+{
+public:
+  /** Convenient typedefs for simplifying declarations. */
+  typedef TPixel InputPixelType;
+  typedef otb::Image<itk::Vector<TPixel,10>,2> InputImageType;//L'image de dpart est de type Mueller
+  typedef otb::Image<itk::Vector<TPixel,3>,2> OutputImageType;// L'image d'arrive contient les images LL, RR et LR
+  typedef otb::Image<itk::Vector<TPixel,10>,2> MuellerImageType;
+
+  /** Extract dimension from input and output image. */
+  itkStaticConstMacro(InputImageDimension, unsigned int,
+                      InputImageType::ImageDimension);
+  itkStaticConstMacro(OutputImageDimension, unsigned int,
+                      OutputImageType::ImageDimension);
+
+  
+  /** Standard class typedefs. */
+  typedef MuellerToCircularPolarisationImageFilter Self;
+  typedef itk::ImageToImageFilter< InputImageType, OutputImageType> 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(MuellerToCircularPolarisationImageFilter, ImageToImageFilter);
+  
+  /** Image typedef support. */
+  typedef typename InputImageType::PixelType InputPixelType;
+  typedef typename OutputImageType::PixelType OutputPixelType;
+  
+  typedef typename InputImageType::RegionType InputImageRegionType;
+  typedef typename OutputImageType::RegionType OutputImageRegionType;
+
+  typedef typename InputImageType::SizeType InputSizeType;
+  typedef typename InputImageType::IndexType IndexType ;
+
+  /** MLC type **/
+  typedef itk::Vector<float,9> MLCType;
+
+  /** Coherency type **/
+  typedef itk::Vector<float,9> CoherencyType;
+  /** Mueller type **/
+  typedef itk::Vector<float,10> MuellerType; 
+
+  /**CircularType **/
+  typedef itk::Vector<float, 3> CircularType;
+  
+  
+  /** MuellerToCircularPolarisationImageFilter needs a larger input requested region than
+   * the output requested region.  As such, MuellerToCircularPolarisationImageFilter 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:
+  MuellerToCircularPolarisationImageFilter();
+  virtual ~MuellerToCircularPolarisationImageFilter() {}
+  void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+  /** MuellerToCircularPolarisationImageFilter 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 );
+
+private:
+  MuellerToCircularPolarisationImageFilter(const Self&); //purposely not implemented
+  void operator=(const Self&); //purposely not implemented
+
+ 
+ 
+
+};
+  
+} // end namespace otb
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "otbMuellerToCircularPolarisationImageFilter.txx"
+#endif
+
+#endif

Code/SARPolarimetry/otbMuellerToCircularPolarisationImageFilter.txx

+/*=========================================================================
+
+  Program:   Insight Segmentation & Registration Toolkit
+  Module:    $otbMuellerToCircularPolarisationImageFilter.txx$
+  Language:  C++
+  Date:      $Date: 24/10/2006 $
+  Version:   $Version: 1.0 $
+
+=========================================================================*/
+#ifndef __MuellerToCircularPolarisationImageFilter_txx
+#define __MuellerToCircularPolarisationImageFilter_txx
+
+#include "otbMuellerToCircularPolarisationImageFilter.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 TPixel>
+  MuellerToCircularPolarisationImageFilter<TPixel>::MuellerToCircularPolarisationImageFilter()
+  {
+    // !! TInputImage should be MuellerType and TOutputImage should be CircularType (otherwise fail)
+  
+  }
+
+  template <class TPixel>
+  void 
+  MuellerToCircularPolarisationImageFilter<TPixel>
+  ::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< InputImageType * >( 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 InputImageType::RegionType inputRequestedRegion;
+    inputRequestedRegion = inputPtr->GetRequestedRegion();
+    
+   
+    
+    // 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__);
+	e.SetLocation(ITK_LOCATION);
+	e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
+	e.SetDataObject(inputPtr);
+	throw e;
+      }
+  }
+
+
+  template< class TPixel>
+  void
+  MuellerToCircularPolarisationImageFilter< TPixel>
+  ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
+			 int threadId)
+  {
+
+    
+
+  float M1, M2, M3, M4, M5, M6, M7, M8, M9, M10;
+    
+  itk::ImageRegionConstIterator<InputImageType> inIt (this->GetInput(),outputRegionForThread);
+  itk::ImageRegionIterator<OutputImageType> outIt (this->GetOutput(),outputRegionForThread);
+    CircularType vectorValue;
+    
+    
+    int numPix =0;
+    
+    for ( inIt.GoToBegin(), outIt.GoToBegin(); !inIt.IsAtEnd(); ++inIt, ++outIt)
+    {
+      M1 = inIt.Get()[0];
+      M2 = inIt.Get()[1];
+      M3 = inIt.Get()[2];
+      M4 = inIt.Get()[3];
+      M5 = inIt.Get()[4];
+      M6 = inIt.Get()[5];
+      M7 = inIt.Get()[6];
+      M8 = inIt.Get()[7];
+      M9 = inIt.Get()[8];
+      M10= inIt.Get()[9];
+      
+
+   //  if (numPix>(128*128/2-128)) {std::cout << M1 << M2 << M3 << M4 << M5 << M6 << M7 << M8 << M9 << M10 <<  std::endl;}
+
+      vectorValue[0]= M1+2*M7+M4;   // LL
+      vectorValue[1]= M1-2*M7+M4;   // RR
+      vectorValue[2]= M1-M4;        // LR
+      
+
+      
+      outIt.Set(vectorValue);
+
+    
+      numPix++;
+    }
+  }
+  
+
+  /**
+   * Standard "PrintSelf" method
+   */
+  template <class TPixel>
+  void
+  MuellerToCircularPolarisationImageFilter<TPixel>
+  ::PrintSelf(
+	      std::ostream& os, 
+	      itk::Indent indent) const
+  {
+    Superclass::PrintSelf( os, indent );
+    
+   }
+  
+} // end namespace otb
+
+#endif