diff --git a/Code/BasicFilters/otbMeanShiftImageFilter.txx b/Code/BasicFilters/otbMeanShiftImageFilter.txx
index 9d8614950eadd46135dcd357fe596d61b9ca16f1..1bd3aa742e13f6f475fc51356f41dfe62fdeaef8 100644
--- a/Code/BasicFilters/otbMeanShiftImageFilter.txx
+++ b/Code/BasicFilters/otbMeanShiftImageFilter.txx
@@ -192,7 +192,7 @@ MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter
template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
void
MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>
-::ThreadedGenerateData(const RegionType& outputRegionForThread, int threadId)
+::ThreadedGenerateData(const RegionType& outputRegionForThread, int itkNotUsed(threadId))
{
// Input and output pointers
typename InputImageType::ConstPointer inputPtr = this->GetInput();
diff --git a/Code/BasicFilters/otbStreamingShrinkImageFilter.cxx b/Code/BasicFilters/otbStreamingShrinkImageFilter.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..3d4bd8d9e92b175790fb5236b101caa10bd7bf1e
--- /dev/null
+++ b/Code/BasicFilters/otbStreamingShrinkImageFilter.cxx
@@ -0,0 +1,103 @@
+/*=========================================================================
+
+ 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.
+
+=========================================================================*/
+#include "otbStreamingShrinkImageFilter.h"
+
+namespace otb
+{
+
+unsigned int
+StreamingShrinkImageRegionSplitter
+::GetNumberOfSplits(const RegionType& region, unsigned int requestedNumber)
+{
+ unsigned int theoricalNbPixelPerTile = region.GetNumberOfPixels() / requestedNumber;
+ unsigned int theoricalTileDimension = static_cast<unsigned int> (vcl_sqrt(static_cast<double>(theoricalNbPixelPerTile)) );
+
+ // Take the previous multiple of m_TileSizeAlignment (eventually generate more splits than requested)
+ m_TileDimension = theoricalTileDimension / m_TileSizeAlignment * m_TileSizeAlignment;
+
+ // Minimal tile size is m_TileSizeAlignment * m_TileSizeAlignment
+ if (m_TileDimension < m_TileSizeAlignment)
+ {
+ otbMsgDevMacro(<< "Using the minimal tile size : " << m_TileSizeAlignment << " * " << m_TileSizeAlignment);
+ m_TileDimension = m_TileSizeAlignment;
+ }
+
+ // Use the computed tile size, and generate (m_TileDimension * 1) tiles
+ const SizeType& regionSize = region.GetSize();
+ m_SplitsPerDimension[0] = (regionSize[0] + m_TileDimension - 1) / m_TileDimension;
+ m_SplitsPerDimension[1] = regionSize[1] / m_TileSizeAlignment;
+
+ unsigned int numPieces = 1;
+ for (unsigned int j = 0; j < ImageDimension; ++j)
+ {
+ numPieces *= m_SplitsPerDimension[j];
+ }
+
+ otbMsgDevMacro(<< "Tile dimension : " << m_TileDimension)
+ otbMsgDevMacro(<< "Number of splits per dimension : " << m_SplitsPerDimension[0] << " " << m_SplitsPerDimension[1])
+
+ return numPieces;
+}
+
+StreamingShrinkImageRegionSplitter::RegionType
+StreamingShrinkImageRegionSplitter
+::GetSplit(unsigned int i, unsigned int numberOfPieces, const RegionType& region)
+{
+ RegionType splitRegion;
+ IndexType splitIndex;
+
+ // Compute the actual number of splits
+ unsigned int numPieces = 1;
+ for (unsigned int j = 0; j < ImageDimension; ++j)
+ {
+ numPieces *= m_SplitsPerDimension[j];
+ }
+
+ if (i >= numPieces)
+ {
+ itkExceptionMacro("Requested split number " << i << " but region contains only " << numPieces << " splits");
+ }
+
+ // Compute the split index in the streaming grid
+ splitIndex[1] = i / m_SplitsPerDimension[0];
+ splitIndex[0] = i % m_SplitsPerDimension[0];
+
+ // Transform the split index to the actual coordinates
+ splitRegion.SetIndex(0, region.GetIndex(0) + m_TileDimension * splitIndex[0]);
+ splitRegion.SetIndex(1, region.GetIndex(1) + m_TileSizeAlignment * splitIndex[1]);
+
+ splitRegion.SetSize(0, m_TileDimension);
+ splitRegion.SetSize(1, 1);
+
+ // Handle the borders
+ splitRegion.Crop(region);
+
+ return splitRegion;
+}
+
+void
+StreamingShrinkImageRegionSplitter
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+ Superclass::PrintSelf(os, indent);
+ os << indent << "SplitsPerDimension : " << m_SplitsPerDimension << std::endl;
+ os << indent << "TileDimension : " << m_TileDimension << std::endl;
+ os << indent << "TileSizeAlignment : " << m_TileSizeAlignment << std::endl;
+}
+
+} // End namespace otb
diff --git a/Code/BasicFilters/otbStreamingShrinkImageFilter.h b/Code/BasicFilters/otbStreamingShrinkImageFilter.h
index c077557e021e7e242895479b509ad395c5efa305..b6b3e9eabc91da4b5c32fc56b3e5d7b8f29e407b 100644
--- a/Code/BasicFilters/otbStreamingShrinkImageFilter.h
+++ b/Code/BasicFilters/otbStreamingShrinkImageFilter.h
@@ -25,11 +25,134 @@
#include "otbPersistentImageFilter.h"
#include "otbPersistentFilterStreamingDecorator.h"
-#include "itkTimeProbe.h"
+#include "otbStreamingManager.h"
namespace otb
{
+class ITK_EXPORT StreamingShrinkImageRegionSplitter : public itk::ImageRegionSplitter<2>
+{
+public:
+ /** Standard class typedefs. */
+ typedef StreamingShrinkImageRegionSplitter Self;
+ typedef itk::ImageRegionSplitter<2> 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(StreamingShrinkImageRegionSplitter, itk::Object);
+
+ /** Dimension of the image available at compile time. */
+ itkStaticConstMacro(ImageDimension, unsigned int, 2);
+
+ /** Dimension of the image available at run time. */
+ static unsigned int GetImageDimension()
+ {
+ return ImageDimension;
+ }
+
+ /** Index typedef support. An index is used to access pixel values. */
+ typedef itk::Index<ImageDimension> IndexType;
+ typedef IndexType::IndexValueType IndexValueType;
+
+ /** Size typedef support. A size is used to define region bounds. */
+ typedef itk::Size<ImageDimension> SizeType;
+ typedef SizeType::SizeValueType SizeValueType;
+
+ /** Region typedef support. */
+ typedef itk::ImageRegion<ImageDimension> RegionType;
+
+ /** How many pieces can the specified region be split? A given region
+ * cannot always be divided into the requested number of pieces. For
+ * instance, if the numberOfPieces exceeds the number of pixels along
+ * a certain dimensions, then some splits will not be possible.
+ */
+ virtual unsigned int GetNumberOfSplits(const RegionType& region,
+ unsigned int requestedNumber);
+
+ /** Get a region definition that represents the ith piece a specified region.
+ * The "numberOfPieces" specified should be less than or equal to what
+ * GetNumberOfSplits() returns. */
+ virtual RegionType GetSplit(unsigned int i, unsigned int numberOfPieces,
+ const RegionType& region);
+
+ itkGetMacro(TileSizeAlignment, unsigned int);
+ itkSetMacro(TileSizeAlignment, unsigned int);
+
+ itkGetMacro(TileDimension, unsigned int);
+
+ itkSetMacro(ShrinkFactor, unsigned int);
+ itkGetMacro(ShrinkFactor, unsigned int);
+
+protected:
+ StreamingShrinkImageRegionSplitter() : m_SplitsPerDimension(0U), m_ShrinkFactor(10) {}
+ virtual ~StreamingShrinkImageRegionSplitter() {}
+ void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+private:
+ StreamingShrinkImageRegionSplitter(const StreamingShrinkImageRegionSplitter &); //purposely not implemented
+ void operator =(const StreamingShrinkImageRegionSplitter&); //purposely not implemented
+
+ itk::FixedArray<unsigned int, ImageDimension> m_SplitsPerDimension;
+ unsigned int m_TileDimension;
+ unsigned int m_TileSizeAlignment;
+ unsigned int m_ShrinkFactor;
+};
+
+
+template <class TInputImage>
+class ITK_EXPORT StreamingShrinkStreamingManager : public StreamingManager<TInputImage>
+{
+public:
+ /** Standard class typedefs. */
+ typedef StreamingShrinkStreamingManager Self;
+ typedef StreamingManager<TInputImage> Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Method for creation through the object factory. */
+ itkNewMacro(Self);
+
+ /** Runtime information support. */
+ itkTypeMacro(StreamingShrinkStreamingManager, StreamingManager);
+
+ typedef TInputImage ImageType;
+ typedef typename ImageType::Pointer ImagePointerType;
+ typedef typename ImageType::RegionType RegionType;
+ typedef typename RegionType::IndexType IndexType;
+ typedef typename RegionType::SizeType SizeType;
+ typedef typename ImageType::InternalPixelType PixelType;
+
+ itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension);
+
+ /** Actually computes the stream divisions, according to the specified streaming mode,
+ * eventually using the input parameter to estimate memory consumption */
+ virtual void PrepareStreaming(itk::DataObject * input, const RegionType ®ion);
+
+ void SetShrinkFactor(unsigned int val)
+ {
+ m_ShrinkFactor = val;
+ }
+
+ unsigned int GetShrinkFactor() const
+ {
+ return m_ShrinkFactor;
+ }
+
+protected:
+ StreamingShrinkStreamingManager();
+ virtual ~StreamingShrinkStreamingManager();
+
+private:
+ StreamingShrinkStreamingManager(const StreamingShrinkStreamingManager &); //purposely not implemented
+ void operator =(const StreamingShrinkStreamingManager&); //purposely not implemented
+
+ unsigned int m_ShrinkFactor;
+};
+
/** \class PersistentShrinkImageFilter
* \brief
@@ -116,9 +239,6 @@ private:
/** The shrink factor */
unsigned int m_ShrinkFactor;
-
- itk::TimeProbe m_Chrono;
-
}; // end of class PersistentStatisticsVectorImageFilter
@@ -156,10 +276,14 @@ public:
typedef TOutputImage OutputImageType;
typedef typename Superclass::FilterType PersistentFilterType;
+ typedef StreamingShrinkStreamingManager<InputImageType> StreamingShrinkStreamingManagerType;
+ typedef typename StreamingShrinkStreamingManagerType::Pointer StreamingShrinkStreamingManagerPointerType;
+
void SetInput(InputImageType * input)
{
this->GetFilter()->SetInput(input);
}
+
const InputImageType * GetInput()
{
return this->GetFilter()->GetInput();
@@ -173,15 +297,29 @@ public:
otbSetObjectMemberMacro(Filter, ShrinkFactor, unsigned int);
otbGetObjectMemberMacro(Filter, ShrinkFactor, unsigned int);
+ virtual void Update(void)
+ {
+ m_StreamingManager->SetShrinkFactor( this->GetFilter()->GetShrinkFactor() );
+ Superclass::Update();
+ }
+
protected:
/** Constructor */
- StreamingShrinkImageFilter() {}
+ StreamingShrinkImageFilter()
+ {
+ // Use a specific StreamingManager implementation
+ m_StreamingManager = StreamingShrinkStreamingManagerType::New();
+ this->GetStreamer()->SetStreamingManager( m_StreamingManager );
+ }
+
/** Destructor */
virtual ~StreamingShrinkImageFilter() {}
private:
StreamingShrinkImageFilter(const Self &); //purposely not implemented
void operator =(const Self&); //purposely not implemented
+
+ StreamingShrinkStreamingManagerPointerType m_StreamingManager;
};
} // End namespace otb
diff --git a/Code/BasicFilters/otbStreamingShrinkImageFilter.txx b/Code/BasicFilters/otbStreamingShrinkImageFilter.txx
index 4c7a810c1bf1d2ea5026ed9822ccd697cae73ca0..cd045e31f27cf8c625ed7bb97dda07de7a163956 100644
--- a/Code/BasicFilters/otbStreamingShrinkImageFilter.txx
+++ b/Code/BasicFilters/otbStreamingShrinkImageFilter.txx
@@ -26,6 +26,34 @@
namespace otb
{
+template <class TImage>
+StreamingShrinkStreamingManager<TImage>::StreamingShrinkStreamingManager()
+{
+}
+
+template <class TImage>
+StreamingShrinkStreamingManager<TImage>::~StreamingShrinkStreamingManager()
+{
+}
+
+template <class TImage>
+void
+StreamingShrinkStreamingManager<TImage>::PrepareStreaming( itk::DataObject * input, const RegionType ®ion )
+{
+ typedef otb::StreamingShrinkImageRegionSplitter TileSplitterType;
+ TileSplitterType::Pointer splitter = TileSplitterType::New();
+ splitter->SetTileSizeAlignment(m_ShrinkFactor);
+ this->m_Splitter = splitter;
+
+ unsigned long nbDivisions = EstimateOptimalNumberOfDivisions(input, region);
+ this->m_ComputedNumberOfSplits = this->m_Splitter->GetNumberOfSplits(region, nbDivisions);
+ otbMsgDevMacro(<< "Number of split : " << this->m_ComputedNumberOfSplits)
+
+ // Save the region to generate the splits later
+ this->m_Region = region;
+}
+
+
/** Constructor */
template <class TInputImage, class TOutputImage>
PersistentShrinkImageFilter<TInputImage, TOutputImage>
@@ -79,31 +107,22 @@ PersistentShrinkImageFilter<TInputImage, TOutputImage>
// Nothing that needs to be allocated for the remaining outputs
}
+
template<class TInputImage, class TOutputImage>
void
PersistentShrinkImageFilter<TInputImage, TOutputImage>
::Reset()
{
- // Reinit the chrono
- m_Chrono = itk::TimeProbe();
-
- // get pointers to the input and output
+ // Get pointers to the input and output
InputImageType* inputPtr = const_cast<InputImageType*>(this->GetInput());
inputPtr->UpdateOutputInformation();
m_ShrinkedOutput = OutputImageType::New();
- // we need to compute the output spacing, the output image size, and the
- // output image start index
const typename InputImageType::SpacingType&
inputSpacing = inputPtr->GetSpacing();
const typename InputImageType::SizeType& inputSize
= inputPtr->GetLargestPossibleRegion().GetSize();
- const typename InputImageType::IndexType& inputStartIndex
- = inputPtr->GetLargestPossibleRegion().GetIndex();
-
- otbMsgDebugMacro(<< "Input index " << inputStartIndex);
- otbMsgDebugMacro(<< "Input size: " << inputSize);
typename OutputImageType::SpacingType shrinkedOutputSpacing;
typename OutputImageType::RegionType shrinkedOutputLargestPossibleRegion;
@@ -113,9 +132,7 @@ PersistentShrinkImageFilter<TInputImage, TOutputImage>
for (unsigned int i = 0; i < OutputImageType::ImageDimension; ++i)
{
shrinkedOutputSpacing[i] = inputSpacing[i] * static_cast<double>(m_ShrinkFactor);
- //shrinkedOutputSize[i] = static_cast<int>(static_cast<double>(inputSize[i]) / static_cast<double>(m_ShrinkFactor));
shrinkedOutputSize[i] = inputSize[i] / m_ShrinkFactor;
- //outputStartIndex[i] = inputStartIndex[i];
shrinkedOutputStartIndex[i] = 0;
}
@@ -134,7 +151,6 @@ void
PersistentShrinkImageFilter<TInputImage, TOutputImage>
::Synthetize()
{
- otbMsgDevMacro( "Shrink time : " << m_Chrono.GetTotal() )
}
template<class TInputImage, class TOutputImage>
@@ -142,7 +158,6 @@ void
PersistentShrinkImageFilter<TInputImage, TOutputImage>
::BeforeThreadedGenerateData()
{
- m_Chrono.Start();
}
template<class TInputImage, class TOutputImage>
@@ -175,7 +190,6 @@ void
PersistentShrinkImageFilter<TInputImage, TOutputImage>
::AfterThreadedGenerateData()
{
- m_Chrono.Stop();
}
template <class TImage, class TOutputImage>
diff --git a/Code/Common/otbFilterWatcherBase.cxx b/Code/Common/otbFilterWatcherBase.cxx
index 2e47d89fc3ba3fc19c98f3368af84b29b705ac8d..742eb90d370ac5ece85d2ae495e7a9e3a5bf48ca 100644
--- a/Code/Common/otbFilterWatcherBase.cxx
+++ b/Code/Common/otbFilterWatcherBase.cxx
@@ -46,11 +46,11 @@ FilterWatcherBase
// Assign the callbacks
m_StartFilterCommand->SetCallbackFunction(this,
- &FilterWatcherBase::StartFilter);
+ &FilterWatcherBase::StartFilterCallback);
m_EndFilterCommand->SetCallbackFunction(this,
- &FilterWatcherBase::EndFilter);
+ &FilterWatcherBase::EndFilterCallback);
m_ProgressFilterCommand->SetCallbackFunction(this,
- &FilterWatcherBase::ShowProgress);
+ &FilterWatcherBase::ShowProgressCallback);
// Add the commands as observers
m_StartTag = m_Process->AddObserver(itk::StartEvent(),
@@ -90,11 +90,11 @@ FilterWatcherBase
// Assign the callbacks
m_StartFilterCommand->SetCallbackFunction(this,
- &FilterWatcherBase::StartFilter);
+ &FilterWatcherBase::StartFilterCallback);
m_EndFilterCommand->SetCallbackFunction(this,
- &FilterWatcherBase::EndFilter);
+ &FilterWatcherBase::EndFilterCallback);
m_ProgressFilterCommand->SetCallbackFunction(this,
- &FilterWatcherBase::ShowProgress);
+ &FilterWatcherBase::ShowProgressCallback);
// Add the commands as observers
m_StartTag = m_Process->AddObserver(itk::StartEvent(), m_StartFilterCommand);
@@ -132,11 +132,11 @@ FilterWatcherBase
// Assign the callbacks
m_StartFilterCommand->SetCallbackFunction(this,
- &FilterWatcherBase::StartFilter);
+ &FilterWatcherBase::StartFilterCallback);
m_EndFilterCommand->SetCallbackFunction(this,
- &FilterWatcherBase::EndFilter);
+ &FilterWatcherBase::EndFilterCallback);
m_ProgressFilterCommand->SetCallbackFunction(this,
- &FilterWatcherBase::ShowProgress);
+ &FilterWatcherBase::ShowProgressCallback);
// Add the commands as observers
m_StartTag = m_Process->AddObserver(itk::StartEvent(), m_StartFilterCommand);
diff --git a/Code/Common/otbFilterWatcherBase.h b/Code/Common/otbFilterWatcherBase.h
index 9e64c176323186f99f421e1e8bc31933df6fbf41..ef739e490c713ab76d91a39705275a4fbd4711ae 100644
--- a/Code/Common/otbFilterWatcherBase.h
+++ b/Code/Common/otbFilterWatcherBase.h
@@ -85,6 +85,24 @@ public:
protected:
+ /** Callback method to show the ProgressEvent */
+ virtual void ShowProgressCallback()
+ {
+ this->ShowProgress();
+ }
+
+ /** Callback method to show the StartEvent */
+ virtual void StartFilterCallback()
+ {
+ this->StartFilter();
+ }
+
+ /** Callback method to show the EndEvent */
+ virtual void EndFilterCallback()
+ {
+ this->EndFilter();
+ }
+
/** Callback method to show the ProgressEvent */
virtual void ShowProgress() = 0;
diff --git a/Code/Common/otbImageRegionSquareTileSplitter.h b/Code/Common/otbImageRegionSquareTileSplitter.h
new file mode 100644
index 0000000000000000000000000000000000000000..10dd5b71e45436d6640a2e1f401d2a46d0ca59ed
--- /dev/null
+++ b/Code/Common/otbImageRegionSquareTileSplitter.h
@@ -0,0 +1,146 @@
+/*=========================================================================
+
+ 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 __otbImageRegionSquareTileSplitter_h
+#define __otbImageRegionSquareTileSplitter_h
+
+#include "itkObject.h"
+#include "itkRegion.h"
+#include "itkImageRegion.h"
+#include "itkImageRegionSplitter.h"
+#include "itkObjectFactory.h"
+#include "itkIndex.h"
+#include "itkSize.h"
+
+namespace otb
+{
+
+/** \class ImageRegionSquareTileSplitter
+ * \brief Divide a region into several pieces.
+ *
+ * ImageRegionSquareTileSplitter divides an ImageRegion into smaller regions.
+ * ImageRegionSquareTileSplitter is used by the StreamingImageFilter to divide a
+ * requested output region into a series of smaller requests of the
+ * pipeline. This object has two basic methods: GetNumberOfSplits()
+ * and GetSplit().
+ *
+ * GetNumberOfSplits() is used to determine how may subregions a given
+ * region can be divided. You call GetNumberOfSplits with an argument
+ * that is the number of subregions you want. If the image region can
+ * support that number of subregions, that number is returned.
+ * Otherwise, the maximum number of splits a region can support will
+ * be returned. For example, if a region splitter class only divides
+ * a region into horizontal slabs, then the maximum number of splits
+ * will be the number of rows in the region.
+ *
+ * GetSplit() returns the ith of N subregions (as an ImageRegion object).
+ *
+ * This ImageRegionSquareTileSplitter class divides a region along the outermost
+ * dimension. If the outermost dimension has size 1 (i.e. a volume
+ * with a single slice), the ImageRegionSquareTileSplitter will divide the
+ * region along the next outermost dimension. If that dimension has size 1,
+ * the process continues with the next outermost dimension.
+ *
+ * Regions obtained by the ImageRegionSquareTileSplitter are aligned on a grid
+ * with width of 256. Divisions can occur only at line defined as k*256.
+ *
+ * Other ImageRegionSquareTileSplitter subclasses could divide an image into
+ * more uniform shaped regions instead of slabs.
+ *
+ * \sa ImageRegionMultidimensionalSplitter
+ *
+ * \ingroup ITKSystemObjects
+ * \ingroup DataProcessing
+ */
+
+template <unsigned int VImageDimension>
+class ITK_EXPORT ImageRegionSquareTileSplitter : public itk::ImageRegionSplitter<VImageDimension>
+{
+public:
+ /** Standard class typedefs. */
+ typedef ImageRegionSquareTileSplitter Self;
+ typedef itk::ImageRegionSplitter<VImageDimension> 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(ImageRegionSquareTileSplitter, itk::Object);
+
+ /** Dimension of the image available at compile time. */
+ itkStaticConstMacro(ImageDimension, unsigned int, VImageDimension);
+
+ /** Dimension of the image available at run time. */
+ static unsigned int GetImageDimension()
+ {
+ return VImageDimension;
+ }
+
+ /** Index typedef support. An index is used to access pixel values. */
+ typedef itk::Index<VImageDimension> IndexType;
+ typedef typename IndexType::IndexValueType IndexValueType;
+
+ /** Size typedef support. A size is used to define region bounds. */
+ typedef itk::Size<VImageDimension> SizeType;
+ typedef typename SizeType::SizeValueType SizeValueType;
+
+ /** Region typedef support. */
+ typedef itk::ImageRegion<VImageDimension> RegionType;
+
+ /** How many pieces can the specified region be split? A given region
+ * cannot always be divided into the requested number of pieces. For
+ * instance, if the numberOfPieces exceeds the number of pixels along
+ * a certain dimensions, then some splits will not be possible.
+ */
+ virtual unsigned int GetNumberOfSplits(const RegionType& region,
+ unsigned int requestedNumber);
+
+ /** Get a region definition that represents the ith piece a specified region.
+ * The "numberOfPieces" specified should be less than or equal to what
+ * GetNumberOfSplits() returns. */
+ virtual RegionType GetSplit(unsigned int i, unsigned int numberOfPieces,
+ const RegionType& region);
+
+ itkGetMacro(TileSizeAlignment, unsigned int);
+ itkSetMacro(TileSizeAlignment, unsigned int);
+
+ itkGetMacro(TileDimension, unsigned int);
+
+protected:
+ ImageRegionSquareTileSplitter() : m_SplitsPerDimension(0U), m_TileDimension(0), m_TileSizeAlignment(16) {}
+ virtual ~ImageRegionSquareTileSplitter() {}
+ void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+private:
+ ImageRegionSquareTileSplitter(const ImageRegionSquareTileSplitter &); //purposely not implemented
+ void operator =(const ImageRegionSquareTileSplitter&); //purposely not implemented
+
+ itk::FixedArray<unsigned int, VImageDimension> m_SplitsPerDimension;
+ unsigned int m_TileDimension;
+ unsigned int m_TileSizeAlignment;
+};
+
+} // end namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+# include "otbImageRegionSquareTileSplitter.txx"
+#endif
+
+#endif
diff --git a/Code/Common/otbImageRegionSquareTileSplitter.txx b/Code/Common/otbImageRegionSquareTileSplitter.txx
new file mode 100644
index 0000000000000000000000000000000000000000..6870b8dd034b2cef25ba00ab21b28721a2aa40bd
--- /dev/null
+++ b/Code/Common/otbImageRegionSquareTileSplitter.txx
@@ -0,0 +1,120 @@
+/*=========================================================================
+
+ 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 __otbImageRegionSquareTileSplitter_txx
+#define __otbImageRegionSquareTileSplitter_txx
+
+#include "otbImageRegionSquareTileSplitter.h"
+#include "otbMath.h"
+#include "otbMacro.h"
+
+namespace otb
+{
+
+template <unsigned int VImageDimension>
+unsigned int
+ImageRegionSquareTileSplitter<VImageDimension>
+::GetNumberOfSplits(const RegionType& region, unsigned int requestedNumber)
+{
+ unsigned int theoricalNbPixelPerTile = region.GetNumberOfPixels() / requestedNumber;
+ unsigned int theoricalTileDimension = static_cast<unsigned int> (vcl_sqrt(static_cast<double>(theoricalNbPixelPerTile)) );
+
+ // Take the previous multiple of m_TileSizeAlignment (eventually generate more splits than requested)
+ m_TileDimension = theoricalTileDimension / m_TileSizeAlignment * m_TileSizeAlignment;
+
+ // Minimal tile size is m_TileSizeAlignment * m_TileSizeAlignment
+ if (m_TileDimension < m_TileSizeAlignment)
+ {
+ otbMsgDevMacro(<< "Warning: clamping tile size to " << m_TileSizeAlignment << " * " << m_TileSizeAlignment);
+ m_TileDimension = m_TileSizeAlignment;
+ }
+
+ unsigned int numPieces = 1;
+ const SizeType& regionSize = region.GetSize();
+ for (unsigned int j = 0; j < VImageDimension; ++j)
+ {
+ m_SplitsPerDimension[j] = (regionSize[j] + m_TileDimension - 1) / m_TileDimension;
+ numPieces *= m_SplitsPerDimension[j];
+ }
+
+ otbMsgDevMacro(<< "Tile dimension : " << m_TileDimension)
+ otbMsgDevMacro(<< "Number of splits per dimension : " << m_SplitsPerDimension[0] << " " << m_SplitsPerDimension[1])
+
+ return numPieces;
+}
+
+template <unsigned int VImageDimension>
+itk::ImageRegion<VImageDimension>
+ImageRegionSquareTileSplitter<VImageDimension>
+::GetSplit(unsigned int i, unsigned int itkNotUsed(numberOfPieces), const RegionType& region)
+{
+ RegionType splitRegion;
+ IndexType splitIndex;
+
+ // Compute the actual number of splits
+ unsigned int numPieces = 1;
+ for (unsigned int j = 0; j < VImageDimension; ++j)
+ {
+ numPieces *= m_SplitsPerDimension[j];
+ }
+
+ // Sanity check
+ if (i >= numPieces)
+ {
+ itkExceptionMacro("Asked for split number " << i << " but region contains only " << numPieces << " splits");
+ }
+
+ // Compute the split index in the streaming grid
+ unsigned int remaining = i;
+ for (unsigned int j = VImageDimension - 1; j > 0; --j)
+ {
+ splitIndex[j] = remaining / m_SplitsPerDimension[VImageDimension - 1 - j];
+ remaining = remaining % m_SplitsPerDimension[VImageDimension - 1 - j];
+ }
+ splitIndex[0] = remaining;
+
+ // Transform the split index to the actual coordinates
+ for (unsigned int j = 0; j < VImageDimension; ++j)
+ {
+ splitRegion.SetIndex(j, region.GetIndex(j) + m_TileDimension * splitIndex[j]);
+ splitRegion.SetSize(j, m_TileDimension);
+ }
+
+ // Handle the borders
+ splitRegion.Crop(region);
+
+ return splitRegion;
+}
+
+/**
+ *
+ */
+template <unsigned int VImageDimension>
+void
+ImageRegionSquareTileSplitter<VImageDimension>
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+ Superclass::PrintSelf(os, indent);
+ os << indent << "SplitsPerDimension : " << m_SplitsPerDimension << std::endl;
+ os << indent << "TileDimension : " << m_TileDimension << std::endl;
+ os << indent << "TileSizeAlignment : " << m_TileSizeAlignment << std::endl;
+
+}
+
+} // end namespace itk
+
+#endif
diff --git a/Code/Common/otbPersistentImageFilter.h b/Code/Common/otbPersistentImageFilter.h
index ccea263130c7e5727c5757df6b5802dc8fc1a4b1..87d78f2202c5f16572dae4d689f5bc005c0e4ce6 100644
--- a/Code/Common/otbPersistentImageFilter.h
+++ b/Code/Common/otbPersistentImageFilter.h
@@ -66,7 +66,7 @@ public:
protected:
/** Constructor */
- PersistentImageFilter() {};
+ PersistentImageFilter() {}
/** Destructor */
virtual ~PersistentImageFilter() {}
/**PrintSelf method */
diff --git a/Code/Common/otbPipelineMemoryPrintCalculator.cxx b/Code/Common/otbPipelineMemoryPrintCalculator.cxx
index 397ba8b1126fbb30d860dccc6d5d8967c3cf5e6c..b8dff8c995087bc77f8d936ed71e8a3c18ad58d7 100644
--- a/Code/Common/otbPipelineMemoryPrintCalculator.cxx
+++ b/Code/Common/otbPipelineMemoryPrintCalculator.cxx
@@ -61,7 +61,7 @@ void
PipelineMemoryPrintCalculator
::Compute()
{
- // Clear the visisted process objects set
+ // Clear the visited process objects set
m_VisitedProcessObjects.clear();
// Dry run of pipeline synchronisation
@@ -96,6 +96,7 @@ PipelineMemoryPrintCalculator::MemoryPrintType
PipelineMemoryPrintCalculator
::EvaluateMemoryPrint(ProcessObjectType * process)
{
+ otbMsgDevMacro(<< "EvaluateMemoryPrint for " << process->GetNameOfClass() << " (" << process << ")")
// This variable will store the final print
MemoryPrintType print = 0;
@@ -150,6 +151,7 @@ PipelineMemoryPrintCalculator::MemoryPrintType
PipelineMemoryPrintCalculator
::EvaluateDataObjectPrint(DataObjectType * data) const
{
+ otbMsgDevMacro(<< "EvaluateMemoryPrint for " << data->GetNameOfClass() << " (" << data << ")")
#define OTB_IMAGE_SIZE_BLOCK(type) \
if(dynamic_cast<itk::Image<type, 2> *>(data) != NULL) \
diff --git a/Code/Common/otbStandardFilterWatcher.cxx b/Code/Common/otbStandardFilterWatcher.cxx
index 1fd3e7d4dc90205145af320af5f45e75f5238914..ce47e316a5e22818ab789ac98a6a59c4994f788e 100644
--- a/Code/Common/otbStandardFilterWatcher.cxx
+++ b/Code/Common/otbStandardFilterWatcher.cxx
@@ -29,6 +29,7 @@ StandardFilterWatcher
: FilterWatcherBase(process, comment)
{
m_StarsCount = 50;
+ m_CurrentNbStars = -1;
}
StandardFilterWatcher
@@ -80,9 +81,14 @@ StandardFilterWatcher
progressPercent = 100;
}
- std::string stars(nbStars, '*');
- std::string blanks(nbBlanks, ' ');
- std::cout << "\rProcessing progress: " << progressPercent << "% [" << stars << blanks << "]" << std::flush;
+ if (nbStars > m_CurrentNbStars)
+ {
+ std::string stars(nbStars, '*');
+ std::string blanks(nbBlanks, ' ');
+ std::cout << "\rProcessing progress: " << progressPercent << "% [" << stars << blanks << "]" << std::flush;
+ }
+
+ m_CurrentNbStars = nbStars;
}
}
diff --git a/Code/Common/otbStandardFilterWatcher.h b/Code/Common/otbStandardFilterWatcher.h
index 1f3cbe6e91d41608d74108a307b628b6e6f3198c..5e4fd7f69037855a477e3e2c9b89a0605bf2bc07 100644
--- a/Code/Common/otbStandardFilterWatcher.h
+++ b/Code/Common/otbStandardFilterWatcher.h
@@ -92,6 +92,8 @@ private:
/** Stars coutning */
int m_StarsCount;
+
+ int m_CurrentNbStars;
};
} // end namespace otb
diff --git a/Code/Common/otbStreamingManager.h b/Code/Common/otbStreamingManager.h
new file mode 100644
index 0000000000000000000000000000000000000000..6d5d81b43aca507d7c8b05df9fb7ba93dc099047
--- /dev/null
+++ b/Code/Common/otbStreamingManager.h
@@ -0,0 +1,196 @@
+/*=========================================================================
+
+ 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 __otbStreamingManager_h
+#define __otbStreamingManager_h
+
+#include "otbMacro.h"
+#include "otbConfigure.h"
+
+#include "itkLightObject.h"
+#include "itkImageRegionSplitter.h"
+
+namespace otb
+{
+
+namespace StreamingManagement
+{
+ enum StreamingMode
+ {
+ /** Estimates the memory used by the pipeline
+ * and set the strips size to fit the specified
+ * available RAM
+ */
+ STRIPPED_AVAILABLE_RAM,
+
+ /** Force the use of a specific number of lines per strips
+ */
+ STRIPPED_SET_NUMBEROFLINES,
+
+ /** Estimates the memory used that will be used during the pipeline
+ * execution and set the tile size to fit the specified
+ * available RAM.
+ */
+ TILED_AVAILABLE_RAM,
+
+ /** Force the use of a specific tile dimension
+ * The associated parameter is the size used in each dimension
+ */
+ TILED_SET_TILE_SIZE,
+ };
+}
+
+/** \class StreamingManager
+ * \brief This class handles the streaming process used in the writers implementation
+ *
+ * The streaming mode can be chosen with either SetStrippedRAMStreamingMode, SetStrippedNumberOfLinesStreamingMode,
+ * SetTiledRAMStreamingMode, or SetTiledTileDimensionStreamingMode.
+ *
+ * Then, PrepareStreaming must be called so that the stream type and dimensions are computed
+ * This involves passing the actual DataObject who will be written, since it will be used
+ * during memory estimation for some specific streaming modes.
+ *
+ * After PrepareStreaming has been called, the actual number of splits and streaming mode which will be used
+ * can be retrieved with GetStreamingMode and GetNumberOfSplits.
+ * The different splits can be retrieved with GetSplit
+ *
+ * \sa StreamingImageFileWriter
+ * \sa StreamingImageVirtualFileWriter
+ */
+template<class TImage>
+class ITK_EXPORT StreamingManager : public itk::LightObject
+{
+public:
+ /** Standard class typedefs. */
+ typedef StreamingManager Self;
+ typedef itk::LightObject Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ typedef TImage ImageType;
+ typedef typename ImageType::Pointer ImagePointerType;
+ typedef typename ImageType::RegionType RegionType;
+ typedef typename RegionType::IndexType IndexType;
+ typedef typename RegionType::SizeType SizeType;
+ typedef typename ImageType::InternalPixelType PixelType;
+
+ typedef StreamingManagement::StreamingMode StreamingModeType;
+
+ /** Creation through object factory macro */
+ itkNewMacro(Self);
+
+ /** Type macro */
+ itkTypeMacro(StreamingManager, itk::LightObject);
+
+ /** Dimension of input image. */
+ itkStaticConstMacro(ImageDimension, unsigned int, ImageType::ImageDimension);
+
+ /** Use stripped mode.
+ * The number of lines of the strips are computed by estimating the
+ * memory footprint of the pipeline, and chosen so that no more than
+ * availableRAM MBytes of memory are used.
+ * If no parameter is given, then the available RAM is retrieved from the
+ * OTB configuration file or build configuration option */
+ virtual void SetStrippedRAMStreamingMode( unsigned int availableRAMInMB = 0 );
+
+ /** Use stripped mode.
+ * The number of lines of the strips are explicitely specified */
+ virtual void SetStrippedNumberOfLinesStreamingMode( unsigned int numberOfLines );
+
+ /** Use tiled mode.
+ * The dimensions of the tiles are computed by estimating the
+ * memory footprint of the pipeline, and chosen so that no more than
+ * availableRAM MBytes of memory are used.
+ *
+ * If no parameter is given, then the available RAM is retrieved from the
+ * OTB configuration file or build configuration option.
+ *
+ * The tiles are set to be square, with dimension aligned with multiple of 16 */
+ virtual void SetTiledRAMStreamingMode( unsigned int availableRAMInMB = 0 );
+
+ /** Use tiled mode.
+ * The dimension of the tile are explicitely specified.
+ * The parameter specifies the size of the tile in each dimension. */
+ virtual void SetTiledTileDimensionStreamingMode( unsigned int tileDimension );
+
+ /** Actually computes the stream divisions, accorfing to the specified streaming mode,
+ * eventually using the input parameter to estimate memory consumption */
+ virtual void PrepareStreaming(itk::DataObject * input, const RegionType ®ion);
+
+ /** Returns the actual streaming mode that will be used to process the image.
+ * PrepareStreaming() must have been called before.
+ * This can be different than the required streaming mode. For example, if
+ * the input passed to PrepareStreaming() is fully buffered, then
+ * the STRIPPED_SET_NUMBEROFLINES mode is used with only one strip */
+ virtual StreamingManagement::StreamingMode GetStreamingMode();
+
+ /** Returns the actual number of pieces that will be used to process the image.
+ * PrepareStreaming() must have been called before.
+ * This can be different than the requested number */
+ virtual unsigned int GetNumberOfSplits();
+
+ /** Get a region definition that represents the ith piece a specified region.
+ * The "numberOfPieces" must be equal to what
+ * GetNumberOfSplits() returns. */
+ virtual RegionType GetSplit(unsigned int i);
+
+
+public:
+ StreamingManager();
+ virtual ~StreamingManager();
+
+ virtual unsigned int EstimateOptimalNumberOfDivisions(itk::DataObject * input, const RegionType ®ion);
+
+ /** The desired streaming mode specified by the user */
+ StreamingManagement::StreamingMode m_DesiredMode;
+
+ /** The actual streaming mode which will be used */
+ StreamingManagement::StreamingMode m_ActualMode;
+
+ /** The available RAM set as parameter when using the STRIPPED_AVAILABLE_RAM or TILED_AVAILABLE_RAM mode */
+ unsigned int m_AvailableRAMInMB;
+
+ /** The desired number of lines when using the STRIPPED_SET_NUMBEROFLINES streaming mode */
+ unsigned int m_DesiredNumberOfLines;
+
+ /** The desired tile dimension when using the TILED_SET_TILE_SIZE streaming mode */
+ unsigned int m_DesiredTileDimension;
+
+ /** The computed number of splits after PrepareStreaming has been called */
+ unsigned int m_ComputedNumberOfSplits;
+
+ /** The region to stream */
+ RegionType m_Region;
+
+ /** The splitter used to compute the different strips */
+ typedef itk::ImageRegionSplitter<itkGetStaticConstMacro(ImageDimension)> AbstractSplitterType;
+ typedef typename AbstractSplitterType::Pointer AbstractSplitterPointerType;
+ AbstractSplitterPointerType m_Splitter;
+
+private:
+ StreamingManager(const StreamingManager &); //purposely not implemented
+ void operator =(const StreamingManager&); //purposely not implemented
+
+};
+
+} // End namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbStreamingManager.txx"
+#endif
+
+#endif
diff --git a/Code/Common/otbStreamingManager.txx b/Code/Common/otbStreamingManager.txx
new file mode 100644
index 0000000000000000000000000000000000000000..86d5929ff1cc066ec5ea7dd3514b4e30db9769cc
--- /dev/null
+++ b/Code/Common/otbStreamingManager.txx
@@ -0,0 +1,298 @@
+/*=========================================================================
+
+ 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 __otbStreamingManager_txx
+#define __otbStreamingManager_txx
+
+#include "otbStreamingManager.h"
+#include "otbMacro.h"
+#include "otbConfigure.h"
+#include "otbConfigurationFile.h"
+#include "otbPipelineMemoryPrintCalculator.h"
+#include "itkExtractImageFilter.h"
+
+#include "otbImageRegionSquareTileSplitter.h"
+#include "itkImageRegionSplitter.h"
+
+namespace otb
+{
+
+template <class TImage>
+StreamingManager<TImage>::StreamingManager()
+ : m_DesiredMode(StreamingManagement::TILED_AVAILABLE_RAM),
+ m_ActualMode(StreamingManagement::TILED_AVAILABLE_RAM),
+ m_AvailableRAMInMB(0),
+ m_DesiredNumberOfLines(0),
+ m_DesiredTileDimension(0),
+ m_ComputedNumberOfSplits(0)
+{
+}
+
+template <class TImage>
+StreamingManager<TImage>::~StreamingManager()
+{
+}
+
+template <class TImage>
+void
+StreamingManager<TImage>::SetStrippedRAMStreamingMode( unsigned int availableRAMInMB )
+{
+ m_DesiredMode = StreamingManagement::STRIPPED_AVAILABLE_RAM;
+ m_AvailableRAMInMB = availableRAMInMB;
+}
+
+template <class TImage>
+void
+StreamingManager<TImage>::SetStrippedNumberOfLinesStreamingMode( unsigned int numberOfLines )
+{
+ m_DesiredMode = StreamingManagement::STRIPPED_SET_NUMBEROFLINES;
+ m_DesiredNumberOfLines = numberOfLines;
+}
+
+template <class TImage>
+void
+StreamingManager<TImage>::SetTiledRAMStreamingMode( unsigned int availableRAMInMB )
+{
+ otbMsgDevMacro(<< "StreamingManager::SetTiledRAMStreamingMode " << availableRAMInMB)
+ m_DesiredMode = StreamingManagement::TILED_AVAILABLE_RAM;
+ m_AvailableRAMInMB = availableRAMInMB;
+}
+
+template <class TImage>
+void
+StreamingManager<TImage>::SetTiledTileDimensionStreamingMode( unsigned int tileDimension )
+{
+ m_DesiredMode = StreamingManagement::TILED_SET_TILE_SIZE;
+ m_DesiredTileDimension = tileDimension;
+}
+
+template <class TImage>
+void
+StreamingManager<TImage>::PrepareStreaming( itk::DataObject * input, const RegionType ®ion )
+{
+ switch (m_DesiredMode)
+ {
+ case StreamingManagement::STRIPPED_AVAILABLE_RAM:
+ {
+ otbMsgDevMacro(<< "Activating STRIPPED_AVAILABLE_RAM streaming mode")
+ unsigned long nbDivisions = EstimateOptimalNumberOfDivisions(input, region);
+ m_Splitter = itk::ImageRegionSplitter<itkGetStaticConstMacro(ImageDimension)>::New();
+ m_ComputedNumberOfSplits = m_Splitter->GetNumberOfSplits(region, nbDivisions);
+
+ otbMsgDevMacro(<< "Number of split : " << m_ComputedNumberOfSplits)
+ }
+ break;
+
+ case StreamingManagement::STRIPPED_SET_NUMBEROFLINES:
+ {
+ otbMsgDevMacro(<< "Activating STRIPPED_SET_NUMBEROFLINES streaming mode")
+ if (m_DesiredNumberOfLines < 1)
+ {
+ itkWarningMacro(<< "DesiredNumberOfLines set to 0 : streaming disabled")
+ }
+
+ /* Calculate number of split */
+ unsigned long numberLinesOfRegion = region.GetSize()[1]; // Y dimension
+ unsigned long nbSplit;
+ if (numberLinesOfRegion > m_DesiredNumberOfLines && m_DesiredNumberOfLines > 0)
+ {
+ nbSplit =
+ static_cast<unsigned long>(vcl_ceil(static_cast<double>(numberLinesOfRegion) /
+ static_cast<double>(m_DesiredNumberOfLines)));
+ }
+ else
+ {
+ // Don't stream
+ nbSplit = 1;
+ }
+
+ m_Splitter = itk::ImageRegionSplitter<itkGetStaticConstMacro(ImageDimension)>::New();
+ m_ComputedNumberOfSplits = m_Splitter->GetNumberOfSplits(region, nbSplit);
+ otbMsgDevMacro(<< "Number of split : " << m_ComputedNumberOfSplits)
+ }
+ break;
+
+ case StreamingManagement::TILED_AVAILABLE_RAM:
+ {
+ otbMsgDevMacro(<< "Activating TILED_AVAILABLE_RAM streaming mode")
+ unsigned long nbDivisions = EstimateOptimalNumberOfDivisions(input, region);
+ m_Splitter = otb::ImageRegionSquareTileSplitter<itkGetStaticConstMacro(ImageDimension)>::New();
+ m_ComputedNumberOfSplits = m_Splitter->GetNumberOfSplits(region, nbDivisions);
+ otbMsgDevMacro(<< "Number of split : " << m_ComputedNumberOfSplits)
+ }
+ break;
+
+ case StreamingManagement::TILED_SET_TILE_SIZE:
+ {
+ if (m_DesiredTileDimension == 0)
+ {
+ itkWarningMacro(<< "DesiredTileDimension is 0 : switching to mode STRIPPED_SET_NUMBEROFLINES, disabling streaming")
+ this->SetStrippedNumberOfLinesStreamingMode(0);
+ this->PrepareStreaming(input, region);
+ break;
+ }
+
+ otbMsgDevMacro(<< "Activating TILED_SET_TILE_SIZE streaming mode")
+ if (m_DesiredTileDimension < 16)
+ {
+ itkWarningMacro(<< "DesiredTileDimension inferior to 16 : use 16 as tile dimension")
+ m_DesiredTileDimension = 16;
+ }
+
+ // Calculate number of split
+ m_Splitter = otb::ImageRegionSquareTileSplitter<itkGetStaticConstMacro(ImageDimension)>::New();
+ unsigned int nbDesiredTiles = itk::Math::Ceil<unsigned int>( double(region.GetNumberOfPixels()) / (m_DesiredTileDimension * m_DesiredTileDimension) );
+ m_ComputedNumberOfSplits = m_Splitter->GetNumberOfSplits(region, nbDesiredTiles);
+ otbMsgDevMacro(<< "Number of split : " << m_ComputedNumberOfSplits)
+ }
+ break;
+ }
+
+ // Save the region to generate the splits later
+ m_Region = region;
+}
+
+
+template <class TImage>
+unsigned int
+StreamingManager<TImage>::EstimateOptimalNumberOfDivisions(itk::DataObject * input, const RegionType ®ion)
+{
+ otbMsgDevMacro(<< "m_AvailableRAMInMB " << m_AvailableRAMInMB)
+ unsigned int availableRAMInBytes = m_AvailableRAMInMB * 1024 * 1024;
+
+ if (availableRAMInBytes == 0)
+ {
+ otbMsgDevMacro(<< "Retrieving available RAM size from configuration")
+ // Retrieve it from the configuration
+ try
+ {
+ typedef otb::ConfigurationFile ConfigurationType;
+ ConfigurationType::Pointer conf = ConfigurationType::GetInstance();
+
+ availableRAMInBytes = conf->GetParameter<unsigned int>(
+ "OTB_STREAM_MAX_SIZE_BUFFER_FOR_STREAMING");
+ }
+ catch(...)
+ {
+ // We should never have to go here if the configuration file is
+ // correct and found.
+ // In case it is not fallback on the cmake
+ // defined constants.
+ availableRAMInBytes = OTB_STREAM_MAX_SIZE_BUFFER_FOR_STREAMING;
+ }
+ }
+
+ otbMsgDevMacro("RAM used to estimate memory footprint : " << availableRAMInBytes / 1024 / 1024 << " MB")
+
+ otb::PipelineMemoryPrintCalculator::Pointer memoryPrintCalculator;
+ memoryPrintCalculator = otb::PipelineMemoryPrintCalculator::New();
+
+ memoryPrintCalculator->SetAvailableMemory( availableRAMInBytes );
+
+ // Trick to avoid having the resampler compute the whole
+ // deformation field
+ double regionTrickFactor = 1;
+ ImageType* inputImage = dynamic_cast<ImageType*>(input);
+ //inputImage = 0;
+ if (inputImage)
+ {
+
+ typedef itk::ExtractImageFilter<ImageType, ImageType> ExtractFilterType;
+ typename ExtractFilterType::Pointer extractFilter = ExtractFilterType::New();
+ extractFilter->SetInput(inputImage);
+
+ // Define a small region to run the memory footprint estimation,
+ // around the image center, 100 pixels wide in each dimension
+ SizeType smallSize;
+ smallSize.Fill(100);
+ IndexType index;
+ index[0] = region.GetIndex()[0] + region.GetSize()[0]/2 - 50;
+ index[1] = region.GetIndex()[1] + region.GetSize()[1]/2 - 50;
+
+ RegionType smallRegion;
+ smallRegion.SetSize(smallSize);
+ smallRegion.SetIndex(index);
+
+ // In case the image is smaller than 100 pixels in a direction
+ smallRegion.Crop(region);
+
+ extractFilter->SetExtractionRegion(smallRegion);
+
+ bool smallRegionSuccess = smallRegion.Crop(region);
+
+ if (smallRegionSuccess)
+ {
+ otbMsgDevMacro("Using an extract to estimate memory : " << smallRegion)
+ // the region is well behaved, inside the largest possible region
+ memoryPrintCalculator->SetDataToWrite(extractFilter->GetOutput());
+ regionTrickFactor = static_cast<double>( region.GetNumberOfPixels() )
+ / static_cast<double>(smallRegion.GetNumberOfPixels());
+
+ memoryPrintCalculator->SetBiasCorrectionFactor(regionTrickFactor);
+ }
+ else
+ {
+ otbMsgDevMacro("Using the input region to estimate memory : " << region)
+ // the region is not well behaved
+ // use the full region
+ memoryPrintCalculator->SetDataToWrite(input);
+ memoryPrintCalculator->SetBiasCorrectionFactor(1.0);
+ }
+
+ memoryPrintCalculator->Compute();
+ }
+ else
+ {
+ // Use the original object to estimate memory footprint
+ memoryPrintCalculator->SetDataToWrite(input);
+ memoryPrintCalculator->SetBiasCorrectionFactor(1.0);
+
+ memoryPrintCalculator->Compute();
+ }
+
+ otbMsgDevMacro( "Estimated Memory print for the full image : "
+ << static_cast<unsigned int>(memoryPrintCalculator->GetMemoryPrint() / 1024 / 1024 ) << std::endl)
+ otbMsgDevMacro( "Optimal number of stream divisions: "
+ << memoryPrintCalculator->GetOptimalNumberOfStreamDivisions() << std::endl)
+
+ return memoryPrintCalculator->GetOptimalNumberOfStreamDivisions();
+}
+
+template <class TImage>
+StreamingManagement::StreamingMode
+StreamingManager<TImage>::GetStreamingMode()
+{
+ return m_ActualMode;
+}
+
+template <class TImage>
+unsigned int
+StreamingManager<TImage>::GetNumberOfSplits()
+{
+ return m_ComputedNumberOfSplits;
+}
+
+template <class TImage>
+typename StreamingManager<TImage>::RegionType
+StreamingManager<TImage>::GetSplit(unsigned int i)
+{
+ return m_Splitter->GetSplit(i, m_ComputedNumberOfSplits, m_Region);
+}
+
+} // End namespace otb
+
+#endif
diff --git a/Code/Gui/otbFltkFilterWatcher.cxx b/Code/Gui/otbFltkFilterWatcher.cxx
index 65ef6e5452a735dac780f2733e76410593149bf7..8136428c449e14a816a93b49b6cf025d22edb686 100644
--- a/Code/Gui/otbFltkFilterWatcher.cxx
+++ b/Code/Gui/otbFltkFilterWatcher.cxx
@@ -27,7 +27,8 @@ FltkFilterWatcher
::FltkFilterWatcher(itk::ProcessObject* process,
int x, int y, int w, int h,
const char *comment)
- : FilterWatcherBase(process, comment)
+ : FilterWatcherBase(process, comment),
+ m_CurrentProgress(0)
{
m_Window = new Fl_Window(x, y, w + 10, h + 10);
m_Window->label(m_Comment.c_str());
@@ -46,4 +47,39 @@ FltkFilterWatcher
delete m_Window;
}
+void
+FltkFilterWatcher
+::StartFilter()
+{
+ m_Window->show();
+ m_Progress->value(0);
+ m_Progress->show();
+ Fl::check();
+}
+
+void
+FltkFilterWatcher
+::ShowProgress()
+{
+ if (m_Process)
+ {
+ double progress = m_Process->GetProgress();
+
+ // Update only at each 0.5 percent
+ if (progress - m_CurrentProgress > 0.005)
+ {
+ m_Progress->value(progress);
+ m_CurrentProgress = progress;
+ Fl::check();
+ }
+ }
+}
+
+void
+FltkFilterWatcher
+::EndFilter()
+{
+ m_Window->hide();
+}
+
} // end namespace otb
diff --git a/Code/Gui/otbFltkFilterWatcher.h b/Code/Gui/otbFltkFilterWatcher.h
index f15a47bcf437d41aaf4bc240c279e1c384010451..c8bc02ce51516d685557265c15eee4f8f847fc77 100644
--- a/Code/Gui/otbFltkFilterWatcher.h
+++ b/Code/Gui/otbFltkFilterWatcher.h
@@ -38,9 +38,6 @@ namespace otb
class ITK_EXPORT FltkFilterWatcher : public FilterWatcherBase
{
public:
- /** Classes that need access to filter's private data */
- // friend class XMLFilterWatcher;
-
/** Constructor. Takes a ProcessObject to monitor and an optional
* comment string that is prepended to each event message. */
FltkFilterWatcher(itk::ProcessObject * process,
@@ -50,35 +47,21 @@ public:
/** Destructor. */
virtual ~FltkFilterWatcher();
- /** Callback method to show the EndEvent */
- virtual void EndFilter()
- {
- m_Window->hide();
- }
-
protected:
-
/** Callback method to show the ProgressEvent */
- virtual void ShowProgress()
- {
- if (m_Process)
- {
- m_Progress->value(m_Process->GetProgress());
- Fl::check();
- }
- }
+ virtual void ShowProgress();
/** Callback method to show the StartEvent */
- virtual void StartFilter()
- {
- m_Window->show();
- m_Progress->show();
- }
+ virtual void StartFilter();
-private:
+ /** Callback method to show the EndEvent */
+ virtual void EndFilter();
+private:
Fl_Window * m_Window;
Fl_Progress * m_Progress;
+
+ double m_CurrentProgress;
};
} // end namespace otb
diff --git a/Code/Gui/otbFltkWriterWatcher.cxx b/Code/Gui/otbFltkWriterWatcher.cxx
index a3779c4ba24d427f5a20683485808ebaf9d2de2f..b99700274ff272789ca907bca81041a67e105b1f 100644
--- a/Code/Gui/otbFltkWriterWatcher.cxx
+++ b/Code/Gui/otbFltkWriterWatcher.cxx
@@ -27,7 +27,9 @@ FltkWriterWatcher
::FltkWriterWatcher(itk::ProcessObject* process,
int x, int y, int w, int h,
const char *comment)
- : WriterWatcherBase(process, comment)
+ : WriterWatcherBase(process, comment),
+ m_CurrentFilterProgress(0),
+ m_CurrentWriterProgress(0)
{
this->BuildGUI(x, y, w, h, comment);
}
@@ -37,7 +39,9 @@ FltkWriterWatcher
itk::ProcessObject* source,
int x, int y, int w, int h,
const char *comment)
- : WriterWatcherBase(process, source, comment)
+ : WriterWatcherBase(process, source, comment),
+ m_CurrentFilterProgress(0),
+ m_CurrentWriterProgress(0)
{
this->BuildGUI(x, y, w, h, comment);
}
@@ -78,4 +82,74 @@ FltkWriterWatcher
delete m_Window;
}
+
+void
+FltkWriterWatcher
+::StartFilter()
+{
+ m_Window->show();
+ m_FilterProgress->show();
+ m_WriterProgress->show();
+ Fl::check();
+}
+
+void
+FltkWriterWatcher
+::ShowFilterProgress()
+{
+ if (m_SourceProcess)
+ {
+ double progress = m_SourceProcess->GetProgress();
+
+ // Update only at each 0.5 percent
+ if (progress - m_CurrentFilterProgress > 0.005)
+ {
+ m_FilterProgress->value(progress);
+ m_CurrentFilterProgress = progress;
+ Fl::check();
+ }
+ }
+}
+
+void
+FltkWriterWatcher
+::EndFilter()
+{
+}
+
+void
+FltkWriterWatcher
+::StartWriter()
+{
+ m_Window->show();
+ m_FilterProgress->show();
+ m_WriterProgress->show();
+ Fl::check();
+}
+
+void
+FltkWriterWatcher
+::ShowWriterProgress()
+{
+ if (m_Process)
+ {
+ double progress = m_Process->GetProgress();
+
+ // Update only at each 0.5 percent
+ if (progress - m_CurrentWriterProgress > 0.005)
+ {
+ m_WriterProgress->value(progress);
+ m_CurrentWriterProgress = progress;
+ Fl::check();
+ }
+ }
+}
+
+void
+FltkWriterWatcher
+::EndWriter()
+{
+ m_Window->hide();
+}
+
} // end namespace otb
diff --git a/Code/Gui/otbFltkWriterWatcher.h b/Code/Gui/otbFltkWriterWatcher.h
index e03ecd03f0586060f5618936177461f1535c1814..224d69a6ea55db32fa13f022e5e5fc720ad8c1e1 100644
--- a/Code/Gui/otbFltkWriterWatcher.h
+++ b/Code/Gui/otbFltkWriterWatcher.h
@@ -55,54 +55,25 @@ public:
/** Destructor. */
virtual ~FltkWriterWatcher();
- /** Callback method to show the EndEvent */
- virtual void EndWriter()
- {
- m_Window->hide();
- }
+protected:
- virtual void EndFilter()
- {}
+ /** Callback method to show the ProgressEvent from the writer */
+ virtual void ShowWriterProgress();
-protected:
+ /** Callback method to show the StartEvent from the writer*/
+ virtual void StartWriter();
+
+ /** Callback method to show the EndEvent from the writer*/
+ virtual void EndWriter();
- /** Callback method to show the ProgressEvent */
- virtual void ShowFilterProgress()
- {
- if (m_SourceProcess)
- {
- m_FilterProgress->value(m_SourceProcess->GetProgress());
- Fl::check();
- }
- }
-
- /** Callback method to show the ProgressEvent */
- virtual void ShowWriterProgress()
- {
- if (m_Process)
- {
- m_WriterProgress->value(m_Process->GetProgress());
- Fl::check();
- }
- }
-
- /** Callback method to show the StartEvent */
- virtual void StartWriter()
- {
- m_Window->show();
- m_FilterProgress->show();
- m_WriterProgress->show();
- Fl::check();
- }
-
- /** Callback method to show the StartEvent */
- virtual void StartFilter()
- {
- m_Window->show();
- m_FilterProgress->show();
- m_WriterProgress->show();
- Fl::check();
- }
+ /** Callback method to show the ProgressEvent from the filter */
+ virtual void ShowFilterProgress();
+
+ /** Callback method to show the StartEvent from the filter*/
+ virtual void StartFilter();
+
+ /** Callback method to show the EndEvent from the filter*/
+ virtual void EndFilter();
void BuildGUI(int x, int y, int w, int h, const char * comment);
@@ -111,6 +82,9 @@ private:
Fl_Window * m_Window;
Fl_Progress * m_WriterProgress;
Fl_Progress * m_FilterProgress;
+
+ double m_CurrentFilterProgress;
+ double m_CurrentWriterProgress;
};
} // end namespace otb
diff --git a/Code/IO/CMakeLists.txt b/Code/IO/CMakeLists.txt
index 54cc874302720870751213a7ca50d27a73402078..437b1e4ce1fd0d4ff6a06665a7e925a2cb190c3d 100644
--- a/Code/IO/CMakeLists.txt
+++ b/Code/IO/CMakeLists.txt
@@ -29,7 +29,7 @@ ADD_LIBRARY(OTBIO ${OTBIO_SRCS})
# PROPERTIES
# LINK_INTERFACE_LIBRARIES ""
# )
-TARGET_LINK_LIBRARIES (OTBIO ${TIFF_LIBRARY} ${GEOTIFF_LIBRARY} ${GDAL_LIBRARY} ${OGR_LIBRARY} ${JPEG_LIBRARY} OTBCommon)
+TARGET_LINK_LIBRARIES (OTBIO ${TIFF_LIBRARY} ${GEOTIFF_LIBRARY} ${GDAL_LIBRARY} ${OGR_LIBRARY} ${JPEG_LIBRARY} OTBCommon OTBBasicFilters)
TARGET_LINK_LIBRARIES (OTBIO otbossim otbossimplugins ITKIO ITKCommon otbkml tinyXML)
IF (OTB_USE_LIBLAS)
diff --git a/Code/IO/otbMapFileProductWriter.txx b/Code/IO/otbMapFileProductWriter.txx
index 1463bb05def7feb40c1bb7f8079eb38dcff4f9b0..f74f55fff44dc25647c890a3cf25f03b7468a24e 100644
--- a/Code/IO/otbMapFileProductWriter.txx
+++ b/Code/IO/otbMapFileProductWriter.txx
@@ -237,6 +237,7 @@ MapFileProductWriter<TInputImage>
m_StreamingShrinkImageFilter->SetShrinkFactor(sampleRatioValue);
m_StreamingShrinkImageFilter->SetInput(m_VectorImage);
+ m_StreamingShrinkImageFilter->Update();
m_VectorRescaleIntensityImageFilter = VectorRescaleIntensityImageFilterType::New();
m_VectorRescaleIntensityImageFilter->SetInput(m_StreamingShrinkImageFilter->GetOutput());
diff --git a/Code/IO/otbStreamingImageVirtualWriter.h b/Code/IO/otbStreamingImageVirtualWriter.h
index 9ae59e84295c01dce381482e85a1852d6378d888..99e72c1c06736795f2aedfb9f58b02bfa16a7e82 100644
--- a/Code/IO/otbStreamingImageVirtualWriter.h
+++ b/Code/IO/otbStreamingImageVirtualWriter.h
@@ -20,8 +20,7 @@
#include "itkMacro.h"
#include "itkImageToImageFilter.h"
-#include "itkImageRegionSplitter.h"
-#include "otbStreamingTraits.h"
+#include "otbStreamingManager.h"
namespace otb
{
@@ -45,7 +44,7 @@ namespace otb
* \sa PersistentStatisticsImageFilter
* \sa PersistentImageStreamingDecorator.
*/
-template <class TInputImage>
+template <class TInputImage, class TStreamingManager = StreamingManager<TInputImage> >
class ITK_EXPORT StreamingImageVirtualWriter : public itk::ImageToImageFilter<TInputImage, TInputImage>
{
public:
@@ -68,86 +67,62 @@ public:
typedef typename InputImageType::PixelType InputImagePixelType;
/** Streaming traits helper typedef */
- typedef StreamingTraits<InputImageType> StreamingTraitsType;
+ typedef TStreamingManager StreamingManagerType;
+ typedef typename StreamingManagerType::Pointer StreamingManagerPointerType;
/** Dimension of input image. */
itkStaticConstMacro(InputImageDimension, unsigned int,
InputImageType::ImageDimension);
- /** Set/Get the image input of this writer. */
- void SetInput(const InputImageType *input);
- const InputImageType * GetInput(void);
- const InputImageType * GetInput(unsigned int idx);
+ StreamingManagerType* GetStreamingManager(void)
+ {
+ return m_StreamingManager;
+ }
- void SetNthInput(unsigned int idx, const InputImageType *input);
-
-
- /** SmartPointer to a region splitting object */
- typedef itk::ImageRegionSplitter<itkGetStaticConstMacro(InputImageDimension)> SplitterType;
- typedef typename SplitterType::Pointer RegionSplitterPointer;
-
- /** Set buffer memory size (in bytes) use to calculate the number of stream divisions */
- void SetBufferMemorySize(unsigned long);
-
- /** Set the buffer number of lines use to calculate the number of stream divisions */
- void SetBufferNumberOfLinesDivisions(unsigned long);
-
- /** The number of stream divisions is calculate by using
- * OTB_STREAM_IMAGE_SIZE_TO_ACTIVATE_STREAMING and
- * OTB_STREAM_MAX_SIZE_BUFFER_FOR_STREAMING cmake variables.
- */
- void SetAutomaticNumberOfStreamDivisions(void);
-
- /** Set the tiling automatic mode for streaming division */
- void SetTilingStreamDivisions(void);
- /** Choose number of divisions in tiling streaming division */
- void SetTilingStreamDivisions(unsigned long);
-
- /** Return the string to indicate the method use to calculate number of stream divisions. */
- std::string GetMethodUseToCalculateNumberOfStreamDivisions(void);
-
- /** Set the number of pieces to divide the input. The upstream pipeline
- * will be executed this many times. */
- void SetNumberOfStreamDivisions(unsigned long);
-
- /** Get the number of pieces to divide the input. The upstream pipeline
- * will be executed this many times. */
- unsigned long GetNumberOfStreamDivisions(void);
-
- /** Set the helper class for dividing the input into chunks. */
- itkSetObjectMacro(RegionSplitter, SplitterType);
-
- /** Get the helper class for dividing the input into chunks. */
- itkGetObjectMacro(RegionSplitter, SplitterType);
-
- /** Type use to define number of divisions */
- typedef StreamingMode CalculationDivisionEnumType;
-
- virtual void GenerateInputRequestedRegion(void);
+ void SetStreamingManager(StreamingManagerType* streamingManager)
+ {
+ m_StreamingManager = streamingManager;
+ }
protected:
StreamingImageVirtualWriter();
+
virtual ~StreamingImageVirtualWriter();
+
void PrintSelf(std::ostream& os, itk::Indent indent) const;
virtual void GenerateData(void);
+ virtual void GenerateInputRequestedRegion(void);
+
private:
StreamingImageVirtualWriter(const StreamingImageVirtualWriter &); //purposely not implemented
void operator =(const StreamingImageVirtualWriter&); //purposely not implemented
- /** This method calculate the number of stream divisions, by using the CalculationDivision type */
- unsigned long CalculateNumberOfStreamDivisions(void);
+ void ObserveSourceFilterProgress(itk::Object* object, const itk::EventObject & event )
+ {
+ if (typeid(event) != typeid(itk::ProgressEvent))
+ {
+ return;
+ }
+
+ itk::ProcessObject* processObject = dynamic_cast<itk::ProcessObject*>(object);
+ if (processObject)
+ m_DivisionProgress = processObject->GetProgress();
+
+ this->UpdateFilterProgress();
+ }
- /** Use to define the method used to calculate number of divisions */
- unsigned long m_BufferMemorySize;
- unsigned long m_BufferNumberOfLinesDivisions;
- unsigned long m_NumberOfStreamDivisions;
+ void UpdateFilterProgress()
+ {
+ this->UpdateProgress( (m_DivisionProgress + m_CurrentDivision) / m_NumberOfDivisions );
+ }
- RegionSplitterPointer m_RegionSplitter;
+ unsigned int m_NumberOfDivisions;
+ unsigned int m_CurrentDivision;
+ float m_DivisionProgress;
- /** Use to determine method of calculation number of divisions */
- CalculationDivisionEnumType m_CalculationDivision;
+ StreamingManagerPointerType m_StreamingManager;
};
} // end namespace otb
diff --git a/Code/IO/otbStreamingImageVirtualWriter.txx b/Code/IO/otbStreamingImageVirtualWriter.txx
index ce0b5938974e80d17265524e360b278f9aa013df..73df70a34af41b0802a90c3e89872ec9c46ddea0 100644
--- a/Code/IO/otbStreamingImageVirtualWriter.txx
+++ b/Code/IO/otbStreamingImageVirtualWriter.txx
@@ -19,231 +19,55 @@
#define __otbStreamingImageVirtualWriter_txx
#include "otbStreamingImageVirtualWriter.h"
-#include "itkCommand.h"
-#include "itkImageRegionIterator.h"
-#include "itkObjectFactoryBase.h"
-#include "itkImageFileWriter.h"
-#include "itkImageRegionMultidimensionalSplitter.h"
-
#include "otbMacro.h"
#include "otbConfigure.h"
+#include "itkCommand.h"
namespace otb
{
-/**
- *
- */
-template <class TInputImage>
-StreamingImageVirtualWriter<TInputImage>
+
+template <class TInputImage, class TStreamingManager>
+StreamingImageVirtualWriter<TInputImage, TStreamingManager>
::StreamingImageVirtualWriter()
{
- m_BufferMemorySize = 0;
- m_BufferNumberOfLinesDivisions = 0;
- // default to 10 divisions
- m_NumberOfStreamDivisions = 0;
- // default to AUTOMATIC_NUMBER_OF_DIVISIONS
- m_CalculationDivision = SET_AUTOMATIC_NUMBER_OF_STREAM_DIVISIONS;
-
- // create default region splitter
- m_RegionSplitter = itk::ImageRegionSplitter<InputImageDimension>::New();
+ m_StreamingManager = StreamingManagerType::New();
}
-/**
- *
- */
-template <class TInputImage>
-StreamingImageVirtualWriter<TInputImage>
+template <class TInputImage, class TStreamingManager>
+StreamingImageVirtualWriter<TInputImage, TStreamingManager>
::~StreamingImageVirtualWriter()
{
}
-/**
- *
- */
-template <class TInputImage>
-void
-StreamingImageVirtualWriter<TInputImage>
-::SetBufferMemorySize(unsigned long memory_size_divisions)
-{
- m_BufferMemorySize = memory_size_divisions;
- m_CalculationDivision = SET_BUFFER_MEMORY_SIZE;
- this->Modified();
-}
-
-/**
- *
- */
-template <class TInputImage>
-void
-StreamingImageVirtualWriter<TInputImage>
-::SetBufferNumberOfLinesDivisions(unsigned long nb_lines_divisions)
-{
- m_BufferNumberOfLinesDivisions = nb_lines_divisions;
- m_CalculationDivision = SET_BUFFER_NUMBER_OF_LINES;
- this->Modified();
-}
-
-/**
- *
- */
-template <class TInputImage>
-void
-StreamingImageVirtualWriter<TInputImage>
-::SetNumberOfStreamDivisions(unsigned long nb_divisions)
-{
- m_NumberOfStreamDivisions = nb_divisions;
- m_CalculationDivision = SET_NUMBER_OF_STREAM_DIVISIONS;
- this->Modified();
-}
-
-/**
- *
- */
-template <class TInputImage>
-void
-StreamingImageVirtualWriter<TInputImage>
-::SetAutomaticNumberOfStreamDivisions(void)
-{
- m_CalculationDivision = SET_AUTOMATIC_NUMBER_OF_STREAM_DIVISIONS;
- this->Modified();
-}
-
-/**
- *
- */
-template <class TInputImage>
+template <class TInputImage, class TStreamingManager>
void
-StreamingImageVirtualWriter<TInputImage>
-::SetTilingStreamDivisions(void)
-{
- m_CalculationDivision = SET_TILING_WITH_SET_AUTOMATIC_NUMBER_OF_STREAM_DIVISIONS;
- m_RegionSplitter = itk::ImageRegionMultidimensionalSplitter<InputImageDimension>::New();
- this->Modified();
-}
-
-template <class TInputImage>
-void
-StreamingImageVirtualWriter<TInputImage>
-::SetTilingStreamDivisions(unsigned long nb_divisions)
-{
- m_CalculationDivision = SET_TILING_WITH_SET_NUMBER_OF_STREAM_DIVISIONS;
- m_NumberOfStreamDivisions = nb_divisions;
- m_RegionSplitter = itk::ImageRegionMultidimensionalSplitter<InputImageDimension>::New();
- this->Modified();
-}
-
-template <class TInputImage>
-void
-StreamingImageVirtualWriter<TInputImage>
-::SetInput(const InputImageType *input)
-{
- this->SetNthInput(0, input);
-}
-
-template <class TInputImage>
-void
-StreamingImageVirtualWriter<TInputImage>
-::SetNthInput(unsigned int idx, const InputImageType *input)
-{
- // ProcessObject is not const_correct so this cast is required here.
- this->itk::ProcessObject::SetNthInput(idx,
- const_cast<TInputImage *>(input));
-}
-
-template <class TInputImage>
-const typename StreamingImageVirtualWriter<TInputImage>::InputImageType *
-StreamingImageVirtualWriter<TInputImage>
-::GetInput(void)
-{
- if (this->GetNumberOfInputs() < 1)
- {
- return 0;
- }
-
- return static_cast<TInputImage*>
- (this->itk::ProcessObject::GetInput(0));
-}
-
-template <class TInputImage>
-const typename StreamingImageVirtualWriter<TInputImage>::InputImageType *
-StreamingImageVirtualWriter<TInputImage>
-::GetInput(unsigned int idx)
+StreamingImageVirtualWriter<TInputImage, TStreamingManager>
+::PrintSelf(std::ostream& os, itk::Indent indent) const
{
- return static_cast<TInputImage*> (this->itk::ProcessObject::GetInput(idx));
+ Superclass::PrintSelf(os, indent);
}
-template <class TInputImage>
+template <class TInputImage, class TStreamingManager>
void
-StreamingImageVirtualWriter<TInputImage>
+StreamingImageVirtualWriter<TInputImage, TStreamingManager>
::GenerateInputRequestedRegion(void)
{
InputImagePointer inputPtr = const_cast<InputImageType *>(this->GetInput(0));
+
+ InputImageRegionType region;
typename InputImageRegionType::SizeType size;
typename InputImageRegionType::IndexType index;
- InputImageRegionType region;
+
index.Fill(0);
size.Fill(0);
region.SetSize(size);
region.SetIndex(index);
inputPtr->SetRequestedRegion(region);
}
-/**
- *
- */
-template <class TInputImage>
-unsigned long
-StreamingImageVirtualWriter<TInputImage>
-::GetNumberOfStreamDivisions(void)
-{
- return (CalculateNumberOfStreamDivisions());
-}
-template<class TInputImage>
-unsigned long
-StreamingImageVirtualWriter<TInputImage>
-::CalculateNumberOfStreamDivisions(void)
-{
- return StreamingTraitsType
- ::CalculateNumberOfStreamDivisions(this->GetInput(),
- this->GetInput()->GetLargestPossibleRegion(),
- m_RegionSplitter,
- m_CalculationDivision,
- m_NumberOfStreamDivisions,
- m_BufferMemorySize,
- m_BufferNumberOfLinesDivisions);
-}
-/**
- *
- */
-template <class TInputImage>
-std::string
-StreamingImageVirtualWriter<TInputImage>
-::GetMethodUseToCalculateNumberOfStreamDivisions(void)
-{
- return (StreamingTraitsType::GetMethodUseToCalculateNumberOfStreamDivisions(m_CalculationDivision));
-}
-/**
- *
- */
-template <class TInputImage>
-void
-StreamingImageVirtualWriter<TInputImage>
-::PrintSelf(std::ostream& os, itk::Indent indent) const
-{
- Superclass::PrintSelf(os, indent);
- os << indent << "Number of stream divisions: " << m_NumberOfStreamDivisions
- << std::endl;
- if (m_RegionSplitter)
- {
- os << indent << "Region splitter:" << m_RegionSplitter << std::endl;
- }
- else
- {
- os << indent << "Region splitter: (none)" << std::endl;
- }
-}
-template<class TInputImage>
+
+template<class TInputImage, class TStreamingManager>
void
-StreamingImageVirtualWriter<TInputImage>
+StreamingImageVirtualWriter<TInputImage, TStreamingManager>
::GenerateData(void)
{
/**
@@ -257,6 +81,7 @@ StreamingImageVirtualWriter<TInputImage>
* Tell all Observers that the filter is starting
*/
this->InvokeEvent(itk::StartEvent());
+
/**
* Grab the input
*/
@@ -267,25 +92,47 @@ StreamingImageVirtualWriter<TInputImage>
* minimum of what the user specified via SetNumberOfStreamDivisions()
* and what the Splitter thinks is a reasonable value.
*/
- unsigned int numDivisions;
- numDivisions = static_cast<unsigned int>(CalculateNumberOfStreamDivisions());
+ m_StreamingManager->PrepareStreaming(inputPtr, outputRegion);
+ m_NumberOfDivisions = m_StreamingManager->GetNumberOfSplits();
+
+ /**
+ * Register to the ProgressEvent of the source filter
+ */
+ // Get the source process object
+ itk::ProcessObject* source = inputPtr->GetSource();
+
+ // Check if source exists
+ if(source)
+ {
+ typedef itk::MemberCommand<Self> CommandType;
+ typedef typename CommandType::Pointer CommandPointerType;
+
+ CommandPointerType command = CommandType::New();
+ command->SetCallbackFunction(this, &Self::ObserveSourceFilterProgress);
+
+ source->AddObserver(itk::ProgressEvent(), command);
+ }
+ else
+ {
+ itkWarningMacro(<< "Could not get the source process object. Progress report might be buggy");
+ }
/**
* Loop over the number of pieces, execute the upstream pipeline on each
* piece, and copy the results into the output image.
*/
InputImageRegionType streamRegion;
- unsigned int piece;
- for (piece = 0;
- piece < numDivisions && !this->GetAbortGenerateData();
- piece++)
+ for (m_CurrentDivision = 0;
+ m_CurrentDivision < m_NumberOfDivisions && !this->GetAbortGenerateData();
+ m_CurrentDivision++, m_DivisionProgress = 0, this->UpdateFilterProgress())
{
- streamRegion = m_RegionSplitter->GetSplit(piece, numDivisions, outputRegion);
+ streamRegion = m_StreamingManager->GetSplit(m_CurrentDivision);
+ otbMsgDevMacro(<< "Processing region : " << streamRegion )
inputPtr->ReleaseData();
inputPtr->SetRequestedRegion(streamRegion);
inputPtr->Update();
- this->UpdateProgress((float) piece / numDivisions);
}
+
/**
* If we ended due to aborting, push the progress up to 1.0 (since
* it probably didn't end there)
@@ -308,11 +155,14 @@ StreamingImageVirtualWriter<TInputImage>
this->GetOutput(idx)->DataHasBeenGenerated();
}
}
+
/**
* Release any inputs if marked for release
*/
this->ReleaseInputs();
}
+
+
} // end namespace otb
#endif
diff --git a/Code/Projections/otbUtmMapProjection.h b/Code/Projections/otbUtmMapProjection.h
index f30ba4e544d5556335f0d9dd1a83dba8b073404c..b031c63ba05f9d1af083d8a87db9774d4dd34a51 100644
--- a/Code/Projections/otbUtmMapProjection.h
+++ b/Code/Projections/otbUtmMapProjection.h
@@ -53,7 +53,7 @@ public:
virtual void SetZone(long zone);
virtual void SetHemisphere(char hemisphere);
virtual int GetZone() const;
- virtual const char GetHemisphere() const;
+ virtual char GetHemisphere() const;
virtual void SetZoneAndHemisphereFromGeoPoint(const InputPointType& geoPoint);
virtual int GetZoneFromGeoPoint(const InputPointType& geoPoint) const;
diff --git a/Code/Projections/otbUtmMapProjection.txx b/Code/Projections/otbUtmMapProjection.txx
index 6d1ffd46ea162ca99625b078f2b89c2b653aba1b..450e166099ed39206dfaeda9483066bc5a1e91db 100644
--- a/Code/Projections/otbUtmMapProjection.txx
+++ b/Code/Projections/otbUtmMapProjection.txx
@@ -68,7 +68,7 @@ int UtmMapProjection<TTransform>
///\return the hemisphere
template <TransformDirection::TransformationDirection TTransform>
-const char UtmMapProjection<TTransform>
+char UtmMapProjection<TTransform>
::GetHemisphere() const
{
char hemisphere = this->m_MapProjection->getHemisphere();
diff --git a/Code/SARPolarimetry/otbHermitianEigenAnalysis.h b/Code/SARPolarimetry/otbHermitianEigenAnalysis.h
deleted file mode 100644
index 16001824ae9bba73af0305911c6044885f74d341..0000000000000000000000000000000000000000
--- a/Code/SARPolarimetry/otbHermitianEigenAnalysis.h
+++ /dev/null
@@ -1,342 +0,0 @@
-/*=========================================================================
-
- 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 __otbHermitianEigenAnalysis_h
-#define __otbHermitianEigenAnalysis_h
-
-#include "itkMacro.h"
-#include "itkVector.h"
-#include "itkArray2D.h"
-
-namespace otb
-{
-
-/** \class HermitianEigenAnalysis
- * \brief Find Eigen values of a real 2D symmetric matrix. It
- * serves as a thread safe alternative to the class:
- * vnl_symmetric_eigensystem, which uses netlib routines.
- *
- * The class is templated over the input matrix, (which is expected to provide
- * access to its elements with the [][] operator), matrix to store eigen
- * values, (must provide write operations on its elements with the [] operator),
- * EigenMatrix to store eigen vectors (must provide write access to its elements
- * with the [][] operator).
- *
- * The SetOrderEigenValues() method can be used to order eigen values (and their
- * corresponding eigen vectors if computed) in ascending order. This is the
- * default ordering scheme. Eigen vectors and values can be obtained without
- * ordering by calling SetOrderEigenValues(false)
- *
- * The SetOrderEigenMagnitudes() method can be used to order eigen values (and
- * their corresponding eigen vectors if computed) by magnitude in ascending order.
- *
- * The user of this class is explicitly supposed to set the dimension of the
- * 2D matrix using the SetDimension() method.
- *
- * The class contains routines taken from netlib sources. (www.netlib.org).
- * netlib/tql1.c
- * netlib/tql2.c
- * netlib/tred1.c
- * netlib/tred2.c
- *
- * Reference:
- * num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
- * wilkinson.
- * handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
- *
- * \ingroup SARPolarimetry
- */
-
-template < typename TMatrix, typename TVector, typename TEigenMatrix=TMatrix >
-class HermitianEigenAnalysis
-{
-public:
- typedef enum {
- OrderByValue=1,
- OrderByMagnitude,
- DoNotOrder
- }EigenValueOrderType;
-
- typedef itk::Vector<float, 2> ComplexVectorType;
- typedef itk::Vector<ComplexVectorType, 3> HermitianVectorType;
- typedef itk::Vector<HermitianVectorType, 3> HermitianMatrixType;
-
-
- HermitianEigenAnalysis():
- m_Dimension(0),
- m_Order(0),
- m_OrderEigenValues(OrderByValue)
- {};
-
- HermitianEigenAnalysis( const unsigned int dimension ):
- m_Dimension(dimension),
- m_Order(dimension),
- m_OrderEigenValues(OrderByValue)
- {};
-
- virtual ~HermitianEigenAnalysis() {};
-
- typedef TMatrix MatrixType;
- typedef TEigenMatrix EigenMatrixType;
- typedef TVector VectorType;
-
-
-
-
- /** Compute Eigen values of A
- * A is any type that overloads the [][] operator and contains the
- * symmetric matrix. In practice only the upper triangle of the
- * matrix will be accessed. (Both itk::Matrix and vnl_matrix
- * overload [][] operator.)
- *
- * 'EigenValues' is any type that overloads the [] operator and will contain
- * the eigen values.
- *
- * No size checking is performed. A is expected to be a square matrix of size
- * m_Dimension. 'EigenValues' is expected to be of length m_Dimension.
- * The matrix is not checked to see if it is symmetric.
- */
- unsigned int ComputeEigenValues(
- const TMatrix & A,
- TVector & EigenValues) const;
-
- /** Compute Eigen values and vectors of A
- * A is any type that overloads the [][] operator and contains the
- * symmetric matrix. In practice only the upper triangle of the
- * matrix will be accessed. (Both itk::Matrix and vnl_matrix
- * overload [][] operator.)
- *
- * 'EigenValues' is any type that overloads the [] operator and will contain
- * the eigen values.
- *
- * 'EigenVectors' is any type that provides access to its elements with the
- * [][] operator. It is expected be of size m_Dimension * m_Dimension.
- *
- * No size checking is performed. A is expected to be a square matrix of size
- * m_Dimension. 'EigenValues' is expected to be of length m_Dimension.
- * The matrix is not checked to see if it is symmetric.
- *
- * Each row of the matrix 'EigenVectors' represents an eigen vector. (unlike MATLAB
- * where the columns of the [EigenVectors, EigenValues] = eig(A) contains the
- * eigenvectors).
- */
- unsigned int ComputeEigenValuesAndVectors(
- const TMatrix & A,
- TVector & EigenValues,
- TEigenMatrix & EigenVectors ) const;
-
-
- /** Matrix order. Defaults to matrix dimension if not set **/
- void SetOrder(const unsigned int n)
- {
- m_Order = n;
- }
-
- /** Get the Matrix order. Will be 0 unless explicitly set, or unless a
- * call to SetDimension has been made in which case it will be the
- * matrix dimension. */
- unsigned int GetOrder() const { return m_Order; }
-
- /** Set/Get methods to order the eigen values in ascending order.
- * This is the default. ie lambda_1 < lambda_2 < ....
- */
- void SetOrderEigenValues( const bool b )
- {
- if (b) { m_OrderEigenValues = OrderByValue; }
- else { m_OrderEigenValues = DoNotOrder; }
- }
- bool GetOrderEigenValues() const { return (m_OrderEigenValues == OrderByValue); }
-
- /** Set/Get methods to order the eigen value magnitudes in ascending order.
- * In other words, |lambda_1| < |lambda_2| < .....
- */
- void SetOrderEigenMagnitudes( const bool b )
- {
- if (b) { m_OrderEigenValues = OrderByMagnitude; }
- else { m_OrderEigenValues = DoNotOrder; }
- }
- bool GetOrderEigenMagnitudes() const { return (m_OrderEigenValues == OrderByMagnitude); }
-
- /** Set the dimension of the input matrix A. A is a square matrix of
- * size m_Dimension. */
- void SetDimension( const unsigned int n )
- {
- m_Dimension = n;
- if (m_Order == 0 )
- {
- m_Order = m_Dimension;
- }
- }
-
- /** Get Matrix dimension, Will be 0 unless explicitly set by a
- * call to SetDimension. */
- unsigned int GetDimension() const { return m_Dimension; }
-
-
-private:
- unsigned int m_Dimension;
- unsigned int m_Order;
- EigenValueOrderType m_OrderEigenValues;
-
-
- /** Reduces a real symmetric matrix to a symmetric tridiagonal matrix using
- * orthogonal similarity transformations.
- * 'inputMatrix' contains the real symmetric input matrix. Only the lower
- * triangle of the matrix need be supplied. The upper triangle is unaltered.
- * 'd' contains the diagonal elements of the tridiagonal matrix.
- * 'e' contains the subdiagonal elements of the tridiagonal matrix in its
- * last n-1 positions. e(1) is set to zero.
- * 'e2' contains the squares of the corresponding elements of e.
- * 'e2' may coincide with e if the squares are not needed.
- * questions and comments should be directed to burton s. garbow.
- * mathematics and computer science div, argonne national laboratory
- * this version dated august 1983.
- *
- * Function Adapted from netlib/tred1.c.
- * [Changed: remove static vars, enforce const correctness.
- * Use vnl routines as necessary].
- * Reference:
- * num. math. 11, 181-195(1968) by martin, reinsch, and wilkinson.
- * handbook for auto. comp., vol.ii-linear algebra, 212-226(1971). */
- void ReduceToTridiagonalMatrix(double *inputMatrix, VectorType &d,
- double *e, double *e2) const;
-
-
-
- /** Reduces a real symmetric matrix to a symmetric tridiagonal matrix using
- * and accumulating orthogonal similarity transformations.
- * 'inputMatrix' contains the real symmetric input matrix. Only the lower
- * triangle of the matrix need be supplied. The upper triangle is unaltered.
- * 'diagonalElements' will contains the diagonal elements of the tridiagonal
- * matrix.
- * 'subDiagonalElements' will contain the subdiagonal elements of the tridiagonal
- * matrix in its last n-1 positions. subDiagonalElements(1) is set to zero.
- * 'transformMatrix' contains the orthogonal transformation matrix produced
- * in the reduction.
- *
- * questions and comments should be directed to burton s. garbow.
- * mathematics and computer science div, argonne national laboratory
- * this version dated august 1983.
- *
- * Function Adapted from netlib/tred2.c.
- * [Changed: remove static vars, enforce const correctness.
- * Use vnl routines as necessary].
- * Reference:
- * num. math. 11, 181-195(1968) by martin, reinsch, and wilkinson.
- * handbook for auto. comp., vol.ii-linear algebra, 212-226(1971). */
- void ReduceToTridiagonalMatrixAndGetTransformation(
- double *inputMatrix, VectorType &diagonalElements,
- double *subDiagonalElements, double *transformMatrix) const;
-
-
-
- /* Finds the eigenvalues of a symmetric tridiagonal matrix by the ql method.
- *
- * On input:
- * 'd' contains the diagonal elements of the input matrix.
- * 'e' contains the subdiagonal elements of the input matrix
- * in its last n-1 positions. e(1) is arbitrary.
- * On Output:
- * 'd' contains the eigenvalues.
- * 'e' has been destroyed.
- *
- * Returns:
- * zero for normal return,
- * j if the j-th eigenvalue has not been
- * determined after 30 iterations.
- *
- *
- * Reference
- * this subroutine is a translation of the algol procedure tql1,
- * num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
- * wilkinson.
- * handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
- *
- * Questions and comments should be directed to burton s. garbow,
- * mathematics and computer science div, argonne national laboratory
- * this version dated august 1983.
- *
- * Function Adapted from netlib/tql1.c.
- * [Changed: remove static vars, enforce const correctness.
- * Use vnl routines as necessary] */
- unsigned int ComputeEigenValuesUsingQL(
- VectorType &d, double *e) const;
-
-
- /* Finds the eigenvalues and eigenvectors of a symmetric tridiagonal matrix
- * by the ql method.
- *
- * On input:
- * 'd' contains the diagonal elements of the input matrix.
- * 'e' contains the subdiagonal elements of the input matrix
- * in its last n-1 positions. e(1) is arbitrary.
- * 'z' contains the transformation matrix produced in the reduction by
- * ReduceToTridiagonalMatrixAndGetTransformation(), if performed. If the
- * eigenvectors of the tridiagonal matrix are desired, z must contain
- * the identity matrix.
-
- * On Output:
- * 'd' contains the eigenvalues.
- * 'e' has been destroyed.
- * 'z' contains orthonormal eigenvectors of the symmetric tridiagonal
- * (or full) matrix.
- *
- * Returns:
- * zero for normal return,
- * j if the j-th eigenvalue has not been
- * determined after 1000 iterations.
- *
- * Reference
- * this subroutine is a translation of the algol procedure tql1,
- * num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
- * wilkinson.
- * handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
- *
- * Questions and comments should be directed to burton s. garbow,
- * mathematics and computer science div, argonne national laboratory
- * this version dated august 1983.
- *
- * Function Adapted from netlib/tql2.c.
- * [Changed: remove static vars, enforce const correctness.
- * Use vnl routines as necessary]
- */
- unsigned int ComputeEigenValuesAndVectorsUsingQL(
- VectorType &d, double *e, double *z) const;
-
-};
-
-template< typename TMatrix, typename TVector, typename TEigenMatrix >
-std::ostream & operator<<(std::ostream& os,
- const HermitianEigenAnalysis< TMatrix, TVector, TEigenMatrix > &s)
-{
- os << "[ClassType: HermitianEigenAnalysis]" << std::endl;
- os << " Dimension : " << s.GetDimension() << std::endl;
- os << " Order : " << s.GetOrder() << std::endl;
- os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
- os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
- return os;
-}
-
-} // end namespace otb
-
-
-#ifndef ITK_MANUAL_INSTANTIATION
-#include "otbHermitianEigenAnalysis.txx"
-#endif
-
-#endif
-
diff --git a/Code/SARPolarimetry/otbHermitianEigenAnalysis.txx b/Code/SARPolarimetry/otbHermitianEigenAnalysis.txx
deleted file mode 100644
index 04d32590c9f90d07754167c52d4efa424cef122c..0000000000000000000000000000000000000000
--- a/Code/SARPolarimetry/otbHermitianEigenAnalysis.txx
+++ /dev/null
@@ -1,983 +0,0 @@
-/*=========================================================================
-
- 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 __otbHermitianEigenAnalysis_txx
-#define __otbHermitianEigenAnalysis_txx
-
-#include "otbHermitianEigenAnalysis.h"
-#include "vnl/vnl_math.h"
-
-namespace otb
-{
-
-template< class TMatrix, class TVector, class TEigenMatrix >
-unsigned int
-HermitianEigenAnalysis< TMatrix, TVector, TEigenMatrix >::
-ComputeEigenValues(const TMatrix & A,
- TVector & D) const
-{
- double *workArea1 = new double[ m_Dimension ];
-
- // Copy the input matrix
- double *inputMatrix = new double [ m_Dimension * m_Dimension ];
-
- unsigned int k = 0;
-
- for( unsigned int row=0; row < m_Dimension; row++ )
- {
- for( unsigned int col=0; col < m_Dimension; col++ )
- {
- inputMatrix[k++] = A(row, col);
- }
- }
-
- ReduceToTridiagonalMatrix( inputMatrix, D, workArea1, workArea1 );
- const unsigned int eigenErrIndex =
- ComputeEigenValuesUsingQL( D, workArea1 );
-
- delete[] workArea1;
- delete[] inputMatrix;
-
- return eigenErrIndex; //index of eigen value that could not be computed
-}
-
-
-/*******************************************************************************
-Routine : diagonalisation
-Authors : Eric POTTIER, Laurent FERRO-FAMIL
-Creation : 01/2002
-Update :
-*-------------------------------------------------------------------------------
-Description : Computes the eigenvectors and eigenvalues of a N*N hermitian
-matrix (with N < 10)
-*-------------------------------------------------------------------------------
-Inputs arguments :
-MatrixDim : Dimension of the Hermitian Matrix (N)
-HermitianMatrix : N*N*2 complex hermitian matrix
-Returned values :
-EigenVect : N*N*2 complex eigenvector matrix
-EigenVal : N elements eigenvalue real vector
-*******************************************************************************/
-
-static void HermitianDiagonalisation( float HM[3][3][2], float EigenVect[3][3][2], float EigenVal[3])
-{
-
- double a[10][10][2], v[10][10][2], d[10], b[10], z[10];
- double w[2], s[2], c[2], titi[2], gc[2], hc[2];
- double sm, tresh, x, toto, e, f, g, h, r, d1, d2;
- int n, p, q, ii, i, j, k;
-
- //n = MatrixDim;
- n=3;
-
-
-
- for (i = 1; i < n + 1; i++)
- {
- for (j = 1; j < n + 1; j++)
- {
- a[i][j][0] = HM[i - 1][j - 1][0];
- a[i][j][1] = HM[i - 1][j - 1][1];
- v[i][j][0] = 0.;
- v[i][j][1] = 0.;
- }
- v[i][i][0] = 1.;
- v[i][i][1] = 0.;
- }
-
- for (p = 1; p < n + 1; p++)
- {
- d[p] = a[p][p][0];
- b[p] = d[p];
- z[p] = 0.;
- }
-
- for (ii = 1; ii < 1000 * n * n; ii++)
- {
- sm = 0.;
- for (p = 1; p < n; p++)
- {
- for (q = p + 1; q < n + 1; q++)
- {
- sm = sm + 2. * sqrt(a[p][q][0] * a[p][q][0] + a[p][q][1] * a[p][q][1]);
- }
- }
- sm = sm / (n * (n - 1));
- if (sm < 1.E-16)
- {
- goto Sortie;
- }
- tresh = 1.E-17;
- if (ii < 4)
- {
- tresh = (long) 0.2 *sm / (n * n);
- }
- x = -1.E-15;
- for (i = 1; i < n; i++)
- {
- for (j = i + 1; j < n + 1; j++)
- {
- toto = sqrt(a[i][j][0] * a[i][j][0] + a[i][j][1] * a[i][j][1]);
- if (x < toto)
- {
- x = toto;
- p = i;
- q = j;
- }
- }
- }
- toto = sqrt(a[p][q][0] * a[p][q][0] + a[p][q][1] * a[p][q][1]);
- if (toto > tresh)
- {
- e = d[p] - d[q];
- w[0] = a[p][q][0];
- w[1] = a[p][q][1];
- g = sqrt(w[0] * w[0] + w[1] * w[1]);
- g = g * g;
- f = sqrt(e * e + 4. * g);
- d1 = e + f;
- d2 = e - f;
- if (fabs(d2) > fabs(d1))
- {
- d1 = d2;
- }
- r = fabs(d1) / sqrt(d1 * d1 + 4. * g);
- s[0] = r;
- s[1] = 0.;
- titi[0] = 2. * r / d1;
- titi[1] = 0.;
- c[0] = titi[0] * w[0] - titi[1] * w[1];
- c[1] = titi[0] * w[1] + titi[1] * w[0];
- r = sqrt(s[0] * s[0] + s[1] * s[1]);
- r = r * r;
- h = (d1 / 2. + 2. * g / d1) * r;
- d[p] = d[p] - h;
- z[p] = z[p] - h;
- d[q] = d[q] + h;
- z[q] = z[q] + h;
- a[p][q][0] = 0.;
- a[p][q][1] = 0.;
-
- for (j = 1; j < p; j++)
- {
- gc[0] = a[j][p][0];
- gc[1] = a[j][p][1];
- hc[0] = a[j][q][0];
- hc[1] = a[j][q][1];
- a[j][p][0] = c[0] * gc[0] - c[1] * gc[1] - s[0] * hc[0] - s[1] * hc[1];
- a[j][p][1] = c[0] * gc[1] + c[1] * gc[0] - s[0] * hc[1] + s[1] * hc[0];
- a[j][q][0] = s[0] * gc[0] - s[1] * gc[1] + c[0] * hc[0] + c[1] * hc[1];
- a[j][q][1] = s[0] * gc[1] + s[1] * gc[0] + c[0] * hc[1] - c[1] * hc[0];
- }
- for (j = p + 1; j < q; j++)
- {
- gc[0] = a[p][j][0];
- gc[1] = a[p][j][1];
- hc[0] = a[j][q][0];
- hc[1] = a[j][q][1];
- a[p][j][0] = c[0] * gc[0] + c[1] * gc[1] - s[0] * hc[0] - s[1] * hc[1];
- a[p][j][1] = c[0] * gc[1] - c[1] * gc[0] + s[0] * hc[1] - s[1] * hc[0];
- a[j][q][0] = s[0] * gc[0] + s[1] * gc[1] + c[0] * hc[0] + c[1] * hc[1];
- a[j][q][1] = -s[0] * gc[1] + s[1] * gc[0] + c[0] * hc[1] - c[1] * hc[0];
- }
- for (j = q + 1; j < n + 1; j++)
- {
- gc[0] = a[p][j][0];
- gc[1] = a[p][j][1];
- hc[0] = a[q][j][0];
- hc[1] = a[q][j][1];
- a[p][j][0] = c[0] * gc[0] + c[1] * gc[1] - s[0] * hc[0] + s[1] * hc[1];
- a[p][j][1] = c[0] * gc[1] - c[1] * gc[0] - s[0] * hc[1] - s[1] * hc[0];
- a[q][j][0] = s[0] * gc[0] + s[1] * gc[1] + c[0] * hc[0] - c[1] * hc[1];
- a[q][j][1] = s[0] * gc[1] - s[1] * gc[0] + c[0] * hc[1] + c[1] * hc[0];
- }
- for (j = 1; j < n + 1; j++)
- {
- gc[0] = v[j][p][0];
- gc[1] = v[j][p][1];
- hc[0] = v[j][q][0];
- hc[1] = v[j][q][1];
- v[j][p][0] = c[0] * gc[0] - c[1] * gc[1] - s[0] * hc[0] - s[1] * hc[1];
- v[j][p][1] = c[0] * gc[1] + c[1] * gc[0] - s[0] * hc[1] + s[1] * hc[0];
- v[j][q][0] = s[0] * gc[0] - s[1] * gc[1] + c[0] * hc[0] + c[1] * hc[1];
- v[j][q][1] = s[0] * gc[1] + s[1] * gc[0] + c[0] * hc[1] - c[1] * hc[0];
- }
- }
- }
-
- Sortie:
-
- for (k = 1; k < n + 1; k++)
- {
- d[k] = 0;
- for (i = 1; i < n + 1; i++)
- {
- for (j = 1; j < n + 1; j++)
- {
- d[k] = d[k] + v[i][k][0] * (HM[i - 1][j - 1][0] * v[j][k][0] - HM[i - 1][j - 1][1] * v[j][k][1]);
- d[k] = d[k] + v[i][k][1] * (HM[i - 1][j - 1][0] * v[j][k][1] + HM[i - 1][j - 1][1] * v[j][k][0]);
- }
- }
- }
-
- for (i = 1; i < n + 1; i++)
- {
- for (j = i + 1; j < n + 1; j++)
- {
- if (d[j] > d[i])
- {
- x = d[i];
- d[i] = d[j];
- d[j] = x;
- for (k = 1; k < n + 1; k++)
- {
- c[0] = v[k][i][0];
- c[1] = v[k][i][1];
- v[k][i][0] = v[k][j][0];
- v[k][i][1] = v[k][j][1];
- v[k][j][0] = c[0];
- v[k][j][1] = c[1];
- }
- }
- }
- }
-
- for (i = 0; i < n; i++)
- {
- EigenVal[i] = d[i + 1];
- for (j = 0; j < n; j++)
- {
- EigenVect[i][j][0] = v[i + 1][j + 1][0];
- EigenVect[i][j][1] = v[i + 1][j + 1][1];
- }
- }
-
-}
-
-
-template< class TMatrix, class TVector, class TEigenMatrix >
-unsigned int
-HermitianEigenAnalysis< TMatrix, TVector, TEigenMatrix >::
-ComputeEigenValuesAndVectors(
- const TMatrix & T,
- TVector & EigenValues,
- TEigenMatrix & EigenVectors ) const
-{
- // double *workArea1 = new double[ m_Dimension ];
-// double *workArea2 = new double [ m_Dimension * m_Dimension ];
-
- // Copy the input matrix
-
- unsigned int k = 0;
-
- float HM[3][3][2];
-
-
-// Passage de la matrice T T1....T9 une matrice 3*3*2 compatible avec la methode HermitianDiagonalisation
-
-
- HM[0][0][0]=T[0]; HM[0][0][1]=0.;
- HM[0][1][0]=T[3]; HM[0][1][1]=-T[4];
- HM[0][2][0]=T[5]; HM[0][2][1]=-T[6];
- HM[1][0][0]=T[3]; HM[1][0][1]=T[4];
- HM[1][1][0]=T[1]; HM[1][1][1]=0.;
- HM[1][2][0]=T[7]; HM[1][2][1]=-T[8];
- HM[2][0][0]=T[5]; HM[2][0][1]=T[6];
- HM[2][1][0]=T[7]; HM[2][1][1]=T[8];
- HM[2][2][0]=T[2]; HM[2][2][1]=0.;
-
-
-
-
- //Acrobatie pour convectir EigenVectors en un float *** et EigenValues en float*
-
- float eigenVect[3][3][2];
- float eigenVal[3];
-
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++)
- {
- eigenVect[i][j][0]=0.;
- eigenVect[i][j][1]=0.;
- }
-
- for (int i=0; i<3; i++)
- eigenVal[i]=0.;
-
-
-
- HermitianDiagonalisation( HM, eigenVect, eigenVal);
-
-
- //std::cout << HM[0][0][0] << " " << HM[0][1][0]<<"+j"<<HM[0][1][1] << " " << HM[0][2][0]<<"+j"<<HM[0][2][1]<<std::endl;
- //std::cout << HM[1][0][0]<<"+j"<<HM[1][0][1] << " " << HM[1][1][0] << " " << HM[1][2][0]<<"+j"<<HM[1][2][1]<<std::endl;
- //std::cout << HM[2][0][0]<<"+j"<<HM[2][0][1] << " " << HM[2][1][0]<<"+j"<< HM[2][1][1] << " " <<HM[2][2][0]<<std::endl;
- //std::cout << "Valeurs propres : " << eigenVal[0] << "...." << eigenVal[1] << "...." << eigenVal[2] << std::endl<< std::endl;
-
-
- // Recuperation des sorties
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++)
- {
- EigenVectors[i][2*j]=eigenVect[i][j][0];
- EigenVectors[i][2*j+1]=eigenVect[i][j][1];
- }
-
-
- for (int i=0; i<3; i++)
- EigenValues[i]=eigenVal[i];
-
-
-
-
-
- int eigenErrIndex=0;
- return eigenErrIndex; //index of eigen value that could not be computed
-}
-
-
-template< class TMatrix, class TVector, class TEigenMatrix >
-void
-HermitianEigenAnalysis< TMatrix, TVector, TEigenMatrix >::
-ReduceToTridiagonalMatrix(double * a, VectorType &d,
- double *e, double *e2) const
-{
- double d__1;
-
- /* Local variables */
- double f, g, h;
- int i, j, k, l;
- double scale;
-
- for (i = 0; i < static_cast< int >(m_Order); ++i)
- {
- d[i] = a[m_Order-1 + i * m_Dimension];
- a[m_Order-1 + i * m_Dimension] = a[i + i * m_Dimension];
- }
-
-
- for (i = m_Order-1; i >= 0; --i)
- {
- l = i - 1;
- h = 0.;
- scale = 0.;
-
- /* .......... scale row (algol tol then not needed) .......... */
- for (k = 0; k <= l; ++k)
- {
- scale += vnl_math_abs(d[k]);
- }
- if (scale == 0.)
- {
- for (j = 0; j <= l; ++j)
- {
- d[j] = a[l + j * m_Dimension];
- a[l + j * m_Dimension] = a[i + j * m_Dimension];
- a[i + j * m_Dimension] = 0.;
- }
- e[i] = 0.;
- e2[i] = 0.;
- continue;
- }
- for (k = 0; k <= l; ++k)
- {
- d[k] /= scale;
- h += d[k] * d[k];
- }
-
- e2[i] = scale * scale * h;
- f = d[l];
- d__1 = vcl_sqrt(h);
- g = (-1.0) * vnl_math_sgn0(f) * vnl_math_abs(d__1);
- e[i] = scale * g;
- h -= f * g;
- d[l] = f - g;
- if (l != 0)
- {
-
- /* .......... form a*u .......... */
- for (j = 0; j <= l; ++j)
- {
- e[j] = 0.;
- }
-
- for (j = 0; j <= l; ++j)
- {
- f = d[j];
- g = e[j] + a[j + j * m_Dimension] * f;
-
- for (k = j+1; k <= l; ++k)
- {
- g += a[k + j * m_Dimension] * d[k];
- e[k] += a[k + j * m_Dimension] * f;
- }
- e[j] = g;
- }
-
- /* .......... form p .......... */
- f = 0.;
-
- for (j = 0; j <= l; ++j)
- {
- e[j] /= h;
- f += e[j] * d[j];
- }
-
- h = f / (h + h);
- /* .......... form q .......... */
- for (j = 0; j <= l; ++j)
- {
- e[j] -= h * d[j];
- }
-
- /* .......... form reduced a .......... */
- for (j = 0; j <= l; ++j)
- {
- f = d[j];
- g = e[j];
-
- for (k = j; k <= l; ++k)
- {
- a[k + j * m_Dimension] = a[k + j * m_Dimension] - f * e[k] - g * d[k];
- }
- }
- }
-
- for (j = 0; j <= l; ++j)
- {
- f = d[j];
- d[j] = a[l + j * m_Dimension];
- a[l + j * m_Dimension] = a[i + j * m_Dimension];
- a[i + j * m_Dimension] = f * scale;
- }
- }
-}
-
-
-template< class TMatrix, class TVector, class TEigenMatrix >
-void
-HermitianEigenAnalysis< TMatrix, TVector, TEigenMatrix >::
-ReduceToTridiagonalMatrixAndGetTransformation(double * a, VectorType &d,
- double *e, double *z) const
-{
- double d__1;
-
- /* Local variables */
- double f, g, h;
- unsigned int i, j, k, l;
- double scale, hh;
-
- for (i = 0; i < m_Order; ++i)
- {
- for (j = i; j < m_Order; ++j)
- {
- z[j + i * m_Dimension] = a[j + i * m_Dimension];
- }
- d[i] = a[m_Order -1 + i * m_Dimension];
- }
-
- for (i = m_Order-1; i > 0; --i)
- {
- l = i - 1;
- h = 0.0;
- scale = 0.0;
-
- /* .......... scale row (algol tol then not needed) .......... */
- for (k = 0; k <= l; ++k)
- {
- scale += vnl_math_abs(d[k]);
- }
- if (scale == 0.0)
- {
- e[i] = d[l];
-
- for (j = 0; j <= l; ++j)
- {
- d[j] = z[l + j * m_Dimension];
- z[i + j * m_Dimension] = 0.0;
- z[j + i * m_Dimension] = 0.0;
- }
- }
- else
- {
- for (k = 0; k <= l; ++k)
- {
- d[k] /= scale;
- h += d[k] * d[k];
- }
-
- f = d[l];
- d__1 = vcl_sqrt(h);
- g = (-1.0) * vnl_math_sgn0(f) * vnl_math_abs(d__1);
- e[i] = scale * g;
- h -= f * g;
- d[l] = f - g;
-
- /* .......... form a*u .......... */
- for (j = 0; j <= l; ++j)
- {
- e[j] = 0.0;
- }
-
- for (j = 0; j <= l; ++j)
- {
- f = d[j];
- z[j + i * m_Dimension] = f;
- g = e[j] + z[j + j * m_Dimension] * f;
-
- for (k = j+1; k <= l; ++k)
- {
- g += z[k + j * m_Dimension] * d[k];
- e[k] += z[k + j * m_Dimension] * f;
- }
-
- e[j] = g;
- }
-
- /* .......... form p .......... */
- f = 0.0;
-
- for (j = 0; j <= l; ++j)
- {
- e[j] /= h;
- f += e[j] * d[j];
- }
-
- hh = f / (h + h);
-
- /* .......... form q .......... */
- for (j = 0; j <= l; ++j)
- {
- e[j] -= hh * d[j];
- }
-
- /* .......... form reduced a .......... */
- for (j = 0; j <= l; ++j)
- {
- f = d[j];
- g = e[j];
-
- for (k = j; k <= l; ++k)
- {
- z[k + j * m_Dimension] = z[k + j * m_Dimension] - f * e[k] - g * d[k];
- }
-
- d[j] = z[l + j * m_Dimension];
- z[i + j * m_Dimension] = 0.0;
- }
- }
-
- d[i] = h;
- }
-
- /* .......... accumulation of transformation matrices .......... */
- for (i = 1; i < m_Order; ++i)
- {
- l = i - 1;
- z[m_Order-1 + l * m_Dimension] = z[l + l * m_Dimension];
- z[l + l * m_Dimension] = 1.0;
- h = d[i];
- if (h != 0.0)
- {
- for (k = 0; k <= l; ++k)
- {
- d[k] = z[k + i * m_Dimension] / h;
- }
-
- for (j = 0; j <= l; ++j)
- {
- g = 0.0;
-
- for (k = 0; k <= l; ++k)
- {
- g += z[k + i * m_Dimension] * z[k + j * m_Dimension];
- }
-
- for (k = 0; k <= l; ++k)
- {
- z[k + j * m_Dimension] -= g * d[k];
- }
- }
- }
-
- for (k = 0; k <= l; ++k)
- {
- z[k + i * m_Dimension] = 0.0;
- }
-
- }
-
- for (i = 0; i < m_Order; ++i)
- {
- d[i] = z[m_Order-1 + i * m_Dimension];
- z[m_Order-1 + i * m_Dimension] = 0.0;
- }
-
- z[m_Order-1 + (m_Order-1) * m_Dimension] = 1.0;
- e[0] = 0.0;
-
-}
-
-
-template< class TMatrix, class TVector, class TEigenMatrix >
-unsigned int
-HermitianEigenAnalysis< TMatrix, TVector, TEigenMatrix >::
-ComputeEigenValuesUsingQL(VectorType &d, double *e) const
-{
-
- const double c_b10 = 1.0;
-
- /* Local variables */
- double c, f, g, h;
- unsigned int i, j, l, m;
- double p, r, s, c2, c3=0.0;
- double s2=0.0;
- double dl1, el1;
- double tst1, tst2;
-
- unsigned int ierr = 0;
- if (m_Order == 1) {
- return 1;
- }
-
- for (i = 1; i < m_Order; ++i) {
- e[i-1] = e[i];
- }
-
- f = 0.;
- tst1 = 0.;
- e[m_Order-1] = 0.;
-
- for (l = 0; l < m_Order; ++l)
- {
- j = 0;
- h = vnl_math_abs(d[l]) + vnl_math_abs(e[l]);
- if (tst1 < h)
- {
- tst1 = h;
- }
- /* .......... look for small sub-diagonal element .......... */
- for (m = l; m < m_Order; ++m)
- {
- tst2 = tst1 + vnl_math_abs(e[m]);
- if (tst2 == tst1)
- {
- break;
- }
- /* .......... e(n) is always zero, so there is no exit */
- /* through the bottom of the loop .......... */
- }
-
- if (m != l)
- {
-
- do
- {
- if (j == 30)
- {
- /* .......... set error -- no convergence to an */
- /* eigenvalue after 30 iterations .......... */
- ierr = l+1;
- return ierr;
- }
- ++j;
- /* .......... form shift .......... */
- g = d[l];
- p = (d[l+1] - g) / (e[l] * 2.);
- r = vnl_math_hypot(p, c_b10);
- d[l] = e[l] / (p + vnl_math_sgn0(p) * vnl_math_abs(r));
- d[l+1] = e[l] * (p + vnl_math_sgn0(p) * vnl_math_abs(r));
- dl1 = d[l+1];
- h = g - d[l];
-
- for (i = l+2; i < m_Order; ++i)
- {
- d[i] -= h;
- }
-
- f += h;
- /* .......... ql transformation .......... */
- p = d[m];
- c = 1.;
- c2 = c;
- el1 = e[l+1];
- s = 0.;
- for (i = m-1; i >= l; --i)
- {
- c3 = c2;
- c2 = c;
- s2 = s;
- g = c * e[i];
- h = c * p;
- r = vnl_math_hypot(p, e[i]);
- e[i+1] = s * r;
- s = e[i] / r;
- c = p / r;
- p = c * d[i] - s * g;
- d[i+1] = h + s * (c * g + s * d[i]);
- if( i == l )
- {
- break;
- }
- }
-
- p = -s * s2 * c3 * el1 * e[l] / dl1;
- e[l] = s * p;
- d[l] = c * p;
- tst2 = tst1 + vnl_math_abs(e[l]);
- } while (tst2 > tst1);
- }
-
- p = d[l] + f;
-
- if( m_OrderEigenValues == OrderByValue )
- {
- // Order by value
- for (i = l; i > 0; --i)
- {
- if (p >= d[i-1])
- break;
- d[i] = d[i-1];
- }
- d[i] = p;
- }
- else if( m_OrderEigenValues == OrderByMagnitude )
- {
- // Order by magnitude.. make eigen values positive
- for (i = l; i > 0; --i)
- {
- if (vnl_math_abs(p) >= vnl_math_abs(d[i-1]))
- break;
- d[i] = vnl_math_abs(d[i-1]);
- }
- d[i] = vnl_math_abs(p);
- }
- else
- {
- d[l] = p;
- }
- }
-
- return ierr; //ierr'th eigen value that couldn't be computed
-
-}
-
-
-template< class TMatrix, class TVector, class TEigenMatrix >
-unsigned int
-HermitianEigenAnalysis< TMatrix, TVector, TEigenMatrix >::
-ComputeEigenValuesAndVectorsUsingQL(VectorType &d, double *e, double *z) const
-{
-
- const double c_b10 = 1.0;
-
- /* Local variables */
- double c, f, g, h;
- unsigned int i, j, k, l, m;
- double p, r, s, c2, c3=0.0;
- double s2=0.0;
- double dl1, el1;
- double tst1, tst2;
-
- unsigned int ierr = 0;
- if (m_Order == 1)
- {
- return 1;
- }
-
- for (i = 1; i < m_Order; ++i)
- {
- e[i - 1] = e[i];
- }
-
- f = 0.0;
- tst1 = 0.;
- e[m_Order-1] = 0.;
-
- for (l = 0; l < m_Order; ++l)
- {
- j = 0;
- h = vnl_math_abs(d[l]) + vnl_math_abs(e[l]);
- if (tst1 < h)
- {
- tst1 = h;
- }
-
- /* .......... look for small sub-diagonal element .......... */
- for (m = l; m < m_Order; ++m)
- {
- tst2 = tst1 + vnl_math_abs(e[m]);
- if (tst2 == tst1)
- {
- break;
- }
-
- /* .......... e(n) is always zero, so there is no exit */
- /* through the bottom of the loop .......... */
- }
-
- if (m != l)
- {
- do
- {
-
- if (j == 1000)
- {
- /* .......... set error -- no convergence to an */
- /* eigenvalue after 1000 iterations .......... */
- ierr = l+1;
- return ierr;
- }
- ++j;
- /* .......... form shift .......... */
- g = d[l];
- p = (d[l+1] - g) / (e[l] * 2.);
- r = vnl_math_hypot(p, c_b10);
- d[l] = e[l] / (p + vnl_math_sgn0(p) * vnl_math_abs(r));
- d[l+1] = e[l] * (p + vnl_math_sgn0(p) * vnl_math_abs(r));
- dl1 = d[l+1];
- h = g - d[l];
-
- for (i = l+2; i < m_Order; ++i)
- {
- d[i] -= h;
- }
-
- f += h;
- /* .......... ql transformation .......... */
- p = d[m];
- c = 1.0;
- c2 = c;
- el1 = e[l+1];
- s = 0.;
-
- for (i = m-1; i >= l; --i)
- {
- c3 = c2;
- c2 = c;
- s2 = s;
- g = c * e[i];
- h = c * p;
- r = vnl_math_hypot(p, e[i]);
- e[i + 1] = s * r;
- s = e[i] / r;
- c = p / r;
- p = c * d[i] - s * g;
- d[i + 1] = h + s * (c * g + s * d[i]);
-
- /* .......... form vector .......... */
- for (k = 0; k < m_Order; ++k)
- {
- h = z[k + (i + 1) * m_Dimension];
- z[k + (i + 1) * m_Dimension] = s * z[k + i * m_Dimension] + c * h;
- z[k + i * m_Dimension] = c * z[k + i * m_Dimension] - s * h;
- }
- if( i == l )
- {
- break;
- }
- }
-
- p = -s * s2 * c3 * el1 * e[l] / dl1;
- e[l] = s * p;
- d[l] = c * p;
- tst2 = tst1 + vnl_math_abs(e[l]);
- } while (tst2 > tst1);
-
- }
-
- d[l] += f;
- }
-
- /* .......... order eigenvalues and eigenvectors .......... */
- if( m_OrderEigenValues == OrderByValue )
- {
- // Order by value
- for (i = 0; i < m_Order-1; ++i)
- {
- k = i;
- p = d[i];
-
- for (j = i+1; j < m_Order; ++j)
- {
- if (d[j] >= p)
- {
- continue;
- }
- k = j;
- p = d[j];
- }
-
- if (k == i)
- {
- continue;
- }
- d[k] = d[i];
- d[i] = p;
-
- for (j = 0; j < m_Order; ++j)
- {
- p = z[j + i * m_Dimension];
- z[j + i * m_Dimension] = z[j + k * m_Dimension];
- z[j + k * m_Dimension] = p;
- }
- }
- }
- else if( m_OrderEigenValues == OrderByMagnitude )
- {
- // Order by magnitude
- for (i = 0; i < m_Order-1; ++i)
- {
- k = i;
- p = d[i];
-
- for (j = i+1; j < m_Order; ++j)
- {
- if (vnl_math_abs(d[j]) >= vnl_math_abs(p))
- {
- continue;
- }
- k = j;
- p = d[j];
- }
-
- if (k == i)
- {
- continue;
- }
- d[k] = vnl_math_abs(d[i]);
- d[i] = vnl_math_abs(p);
-
- for (j = 0; j < m_Order; ++j)
- {
- p = z[j + i * m_Dimension];
- z[j + i * m_Dimension] = z[j + k * m_Dimension];
- z[j + k * m_Dimension] = p;
- }
- }
- }
-
-
- return ierr;
-}
-
-
-} // end namespace otb
-
-#endif
-
diff --git a/Code/SARPolarimetry/otbReciprocalHAlphaImageFilter.h b/Code/SARPolarimetry/otbReciprocalHAlphaImageFilter.h
index 18292eca3bce5cb5dbd14873378fb639edb9bd31..04863bc98651081ceb15499e3bdbfb982eecbcdc 100644
--- a/Code/SARPolarimetry/otbReciprocalHAlphaImageFilter.h
+++ b/Code/SARPolarimetry/otbReciprocalHAlphaImageFilter.h
@@ -20,9 +20,9 @@
#define __ReciprocalHAlphaImageFilter_h
#include "otbUnaryFunctorImageFilter.h"
-#include "otbHermitianEigenAnalysis.h"
-#include "itkVector.h"
#include "otbMath.h"
+#include "vnl/algo/vnl_complex_eigensystem.h"
+#include <algorithm>
namespace otb
{
@@ -44,60 +44,68 @@ template< class TInput, class TOutput>
class ReciprocalHAlphaFunctor
{
public:
- typedef typename std::complex <double> ComplexType;
- typedef typename TOutput::ValueType OutputValueType;
-
- /** CoherencyMatrix type **/
- typedef itk::Vector<double, 9> CoherencyMatrixType;
-
- /** Vector type used to store eigenvalues. */
- typedef itk::Vector<double, 3> EigenvalueType;
-
- /** Matrix type used to store eigenvectors. */
- typedef itk::Vector<float, 2> VectorType;
- typedef itk::Vector<VectorType, 3> EigenVectorFirstComposantType;
- typedef itk::Vector<VectorType, 3> EigenVectorType; // eigenvector type (real part, imaginary part)
- typedef itk::Vector<itk::Vector<float, 6>, 3> EigenMatrixType;
- typedef itk::Image<EigenVectorType, 2> EigenVectorImageType;
- typedef itk::Image<double, 2> EigenValueImageType;
-
- typedef itk::Vector<double, 3> OutputVectorType;
-
- typedef itk::Vector<float, 2> ComplexVectorType;
- typedef itk::Vector<ComplexVectorType, 3> HermitianVectorType;
- typedef itk::Vector<HermitianVectorType, 3> HermitianMatrixType;
- typedef otb::HermitianEigenAnalysis<CoherencyMatrixType, EigenvalueType, EigenMatrixType> HermitianAnalysisType;
-
-
+ typedef typename std::complex<double> ComplexType;
+ typedef vnl_matrix<ComplexType> VNLMatrixType;
+ typedef vnl_vector<ComplexType> VNLVectorType;
+ typedef vnl_vector<double> VNLDoubleVectorType;
+ typedef std::vector<double> VectorType;
+ typedef typename TOutput::ValueType OutputValueType;
+
+
inline TOutput operator()( const TInput & Coherency ) const
{
TOutput result;
result.SetSize(m_NumberOfComponentsPerPixel);
- CoherencyMatrixType T;
- EigenvalueType eigenValues;
- EigenMatrixType eigenVectors;
- HermitianAnalysisType HermitianAnalysis;
-
- T[0] = static_cast<double>(Coherency[0].real());
- T[1] = static_cast<double>(Coherency[3].real());
- T[2] = static_cast<double>(Coherency[5].real());
- T[3] = static_cast<double>(Coherency[1].real());
- T[4] = static_cast<double>(Coherency[1].imag());
- T[5] = static_cast<double>(Coherency[2].real());
- T[6] = static_cast<double>(Coherency[2].imag());
- T[7] = static_cast<double>(Coherency[4].real());
- T[8] = static_cast<double>(Coherency[4].imag());
- HermitianAnalysis.ComputeEigenValuesAndVectors(T, eigenValues, eigenVectors);
-
+ const double T0 = static_cast<double>(Coherency[0].real());
+ const double T1 = static_cast<double>(Coherency[3].real());
+ const double T2 = static_cast<double>(Coherency[5].real());
+
+ VNLMatrixType vnlMat(3, 3, 0.);
+ vnlMat[0][0] = ComplexType(T0, 0.);
+ vnlMat[0][1] = std::conj(ComplexType(Coherency[1]));
+ vnlMat[0][2] = std::conj(ComplexType(Coherency[2]));
+ vnlMat[1][0] = ComplexType(Coherency[1]);
+ vnlMat[1][1] = ComplexType(T1, 0.);
+ vnlMat[1][2] = std::conj(ComplexType(Coherency[4]));
+ vnlMat[2][0] = ComplexType(Coherency[2]);
+ vnlMat[2][1] = ComplexType(Coherency[4]);
+ vnlMat[2][2] = ComplexType(T2, 0.);
+
+ // Only compute the left symetry to respect the previous Hermitian Analisys code
+ vnl_complex_eigensystem syst(vnlMat, false, true);
+ const VNLMatrixType eigenVectors( syst.L );
+ const VNLVectorType eigenValues(syst.W);
+
// Entropy estimation
double totalEigenValues(0.0);
double p[3];
double entropy;
double alpha;
double anisotropy;
-
- totalEigenValues = static_cast<double>( eigenValues[0] + eigenValues[1] + eigenValues[2]);
+
+ // Sort eigen values in decreasing order
+ VectorType sortedRealEigenValues(3, eigenValues[0].real());
+ sortedRealEigenValues[1] = eigenValues[1].real();
+ sortedRealEigenValues[2] = eigenValues[2].real();
+ std::sort( sortedRealEigenValues.begin(), sortedRealEigenValues.end() );
+ std::reverse( sortedRealEigenValues.begin(), sortedRealEigenValues.end() );
+
+ // Extract the first component of each the eigen vector sorted by eigen value decrease order
+ VNLVectorType sortedGreaterEigenVector(3, eigenVectors[0][0]);
+ for(unsigned int i=0; i<3; i++)
+ {
+ if( vcl_abs( eigenValues[1].real()-sortedRealEigenValues[i] ) < m_Epsilon )
+ {
+ sortedGreaterEigenVector[i] = eigenVectors[1][0];
+ }
+ else if( vcl_abs( eigenValues[2].real()-sortedRealEigenValues[i] ) < m_Epsilon )
+ {
+ sortedGreaterEigenVector[i] = eigenVectors[2][0];
+ }
+ }
+
+ totalEigenValues = sortedRealEigenValues[0] + sortedRealEigenValues[1] + sortedRealEigenValues[2];
if (totalEigenValues <m_Epsilon)
{
totalEigenValues = m_Epsilon;
@@ -105,12 +113,7 @@ public:
for (unsigned int k = 0; k < 3; k++)
{
- if (eigenValues[k] <0.)
- {
- eigenValues[k] = 0.;
- }
-
- p[k] = static_cast<double>(eigenValues[k]) / totalEigenValues;
+ p[k] = std::max(sortedRealEigenValues[k], 0.) / totalEigenValues;
}
if ( (p[0] < m_Epsilon) || (p[1] < m_Epsilon) || (p[2] < m_Epsilon) )
@@ -129,7 +132,7 @@ public:
for(unsigned int k = 0; k < 3; k++)
{
- p[k] = eigenValues[k] / totalEigenValues;
+ p[k] = sortedRealEigenValues[k] / totalEigenValues;
if (p[k] < 0.)
{
@@ -142,13 +145,13 @@ public:
}
}
- val0=sqrt(static_cast<double>(eigenVectors[0][0]*eigenVectors[0][0]) + static_cast<double>(eigenVectors[0][1]*eigenVectors[0][1]));
+ val0 = std::abs(sortedGreaterEigenVector[0]);
a0=acos(vcl_abs(val0)) * CONST_180_PI;
- val1=sqrt(static_cast<double>(eigenVectors[0][2]*eigenVectors[0][2]) + static_cast<double>(eigenVectors[0][3]*eigenVectors[0][3]));
+ val1 = std::abs(sortedGreaterEigenVector[1]);
a1=acos(vcl_abs(val1)) * CONST_180_PI;
- val2=sqrt(static_cast<double>(eigenVectors[0][4]*eigenVectors[0][4]) + static_cast<double>(eigenVectors[0][5]*eigenVectors[0][5]));
+ val2= std::abs(sortedGreaterEigenVector[2]);
a2=acos(vcl_abs(val2)) * CONST_180_PI;
alpha=p[0]*a0 + p[1]*a1 + p[2]*a2;
@@ -159,11 +162,11 @@ public:
}
// Anisotropy estimation
- anisotropy=(eigenValues[1] - eigenValues[2])/(eigenValues[1] + eigenValues[2]+m_Epsilon);
+ anisotropy=(sortedRealEigenValues[1] - sortedRealEigenValues[2])/(sortedRealEigenValues[1] + sortedRealEigenValues[2] + m_Epsilon);
- result[0] = entropy;
- result[1] = alpha;
- result[2] = anisotropy;
+ result[0] = static_cast<OutputValueType>(entropy);
+ result[1] = static_cast<OutputValueType>(alpha);
+ result[2] = static_cast<OutputValueType>(anisotropy);
return result;
}
diff --git a/Testing/Code/BasicFilters/otbMatrixImageFilterTest.cxx b/Testing/Code/BasicFilters/otbMatrixImageFilterTest.cxx
index 60c2881d52cccbf547e65e3942088de115aacdb8..483d6afb84f6a5749726343a7fdcb0da0601e64f 100644
--- a/Testing/Code/BasicFilters/otbMatrixImageFilterTest.cxx
+++ b/Testing/Code/BasicFilters/otbMatrixImageFilterTest.cxx
@@ -28,9 +28,6 @@
int otbMatrixImageFilterNew(int argc, char * argv[])
{
- const char * inputFilename = argv[1];
- const char * outputFilename = argv[2];
-
typedef std::complex<double> PixelType;
typedef otb::VectorImage<PixelType> ImageType;
diff --git a/Testing/Code/BasicFilters/otbMatrixTransposeMatrixImageFilter.cxx b/Testing/Code/BasicFilters/otbMatrixTransposeMatrixImageFilter.cxx
index beebc71c1231462c0535364c697f48eab7e090c4..91f728949ad2065b9788f7ee83f2da0b14ca72d0 100644
--- a/Testing/Code/BasicFilters/otbMatrixTransposeMatrixImageFilter.cxx
+++ b/Testing/Code/BasicFilters/otbMatrixTransposeMatrixImageFilter.cxx
@@ -49,8 +49,7 @@ int otbMatrixTransposeMatrixImageFilter(int argc, char * argv[])
reader1->SetFileName(infname1);
reader2->SetFileName(infname2);
- // filter->SetStreamingMode(otb::SET_NUMBER_OF_STREAM_DIVISIONS);
- filter->GetStreamer()->SetNumberOfStreamDivisions(200);
+ filter->GetStreamer()->GetStreamingManager()->SetStrippedNumberOfLinesStreamingMode(200);
filter->SetFirstInput(reader1->GetOutput());
filter->SetSecondInput(reader2->GetOutput());
filter->SetUsePadFirstInput(true);
diff --git a/Testing/Code/BasicFilters/otbStreamingInnerProductVectorImageFilter.cxx b/Testing/Code/BasicFilters/otbStreamingInnerProductVectorImageFilter.cxx
index 2987a0f63604e0d6c49050e4f7a50c7fef9babd6..842fde27ecac9a100baad9058c6ea2471932672c 100644
--- a/Testing/Code/BasicFilters/otbStreamingInnerProductVectorImageFilter.cxx
+++ b/Testing/Code/BasicFilters/otbStreamingInnerProductVectorImageFilter.cxx
@@ -39,8 +39,7 @@ int otbStreamingInnerProductVectorImageFilter(int argc, char* argv[])
// Instantiation object
FilterType::Pointer filter = FilterType::New();
-
- filter->GetStreamer()->SetNumberOfStreamDivisions(10);
+ filter->GetStreamer()->GetStreamingManager()->SetStrippedNumberOfLinesStreamingMode(10);
filter->SetCenterData(centerdata);
filter->SetInput(reader->GetOutput());
filter->Update();
diff --git a/Testing/Code/BasicFilters/otbStreamingMinMaxImageFilter.cxx b/Testing/Code/BasicFilters/otbStreamingMinMaxImageFilter.cxx
index 0a4493c8306a063c2fb731fcc211a1684372a1a9..8cb67d368012f70209198445bf2ddf1dc34319c2 100644
--- a/Testing/Code/BasicFilters/otbStreamingMinMaxImageFilter.cxx
+++ b/Testing/Code/BasicFilters/otbStreamingMinMaxImageFilter.cxx
@@ -42,8 +42,7 @@ int otbStreamingMinMaxImageFilter(int argc, char * argv[])
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName(infname);
- //filter->SetStreamingMode(otb::SET_NUMBER_OF_STREAM_DIVISIONS);
- filter->GetStreamer()->SetNumberOfStreamDivisions(200);
+ filter->GetStreamer()->GetStreamingManager()->SetStrippedNumberOfLinesStreamingMode(200);
filter->SetInput(reader->GetOutput());
otb::StandardFilterWatcher watcher(filter, "Min Max Computation");
filter->Update();
diff --git a/Testing/Code/BasicFilters/otbStreamingMinMaxVectorImageFilter.cxx b/Testing/Code/BasicFilters/otbStreamingMinMaxVectorImageFilter.cxx
index 778df22c3a34f2eddc7ae02ed10113fb4bdd8985..d114e5d7a0d89ecd86b1d70033843bdbd0850126 100644
--- a/Testing/Code/BasicFilters/otbStreamingMinMaxVectorImageFilter.cxx
+++ b/Testing/Code/BasicFilters/otbStreamingMinMaxVectorImageFilter.cxx
@@ -42,8 +42,7 @@ int otbStreamingMinMaxVectorImageFilter(int argc, char * argv[])
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName(infname);
- //filter->SetStreamingMode(otb::SET_NUMBER_OF_STREAM_DIVISIONS);
- filter->GetStreamer()->SetNumberOfStreamDivisions(200);
+ filter->GetStreamer()->GetStreamingManager()->SetStrippedNumberOfLinesStreamingMode(200);
filter->SetInput(reader->GetOutput());
otb::StandardFilterWatcher watcher(filter, "Min Max Computation");
filter->Update();
diff --git a/Testing/Code/BasicFilters/otbStreamingStatisticsImageFilter.cxx b/Testing/Code/BasicFilters/otbStreamingStatisticsImageFilter.cxx
index adb2fc6c71399e1afc7f9dd2c4d40aa702d292a5..474097031e76b937fca9336ec1b2084ea68e249a 100644
--- a/Testing/Code/BasicFilters/otbStreamingStatisticsImageFilter.cxx
+++ b/Testing/Code/BasicFilters/otbStreamingStatisticsImageFilter.cxx
@@ -42,8 +42,7 @@ int otbStreamingStatisticsImageFilter(int argc, char * argv[])
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName(infname);
- //filter->GetStreamer()->SetStreamingMode(otb::SET_NUMBER_OF_STREAM_DIVISIONS);
- filter->GetStreamer()->SetNumberOfStreamDivisions(200);
+ filter->GetStreamer()->GetStreamingManager()->SetStrippedNumberOfLinesStreamingMode(200);
filter->SetInput(reader->GetOutput());
filter->Update();
diff --git a/Testing/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.cxx b/Testing/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.cxx
index e0c312c80be531fb50907cce675f0ec7d0bbd90a..6ff493c558e7fa82c7b61fafb3d4675028808e0a 100644
--- a/Testing/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.cxx
+++ b/Testing/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.cxx
@@ -41,8 +41,7 @@ int otbStreamingStatisticsVectorImageFilter(int argc, char * argv[])
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName(infname);
- //filter->SetStreamingMode(otb::SET_NUMBER_OF_STREAM_DIVISIONS);
- filter->GetStreamer()->SetNumberOfStreamDivisions(200);
+ filter->GetStreamer()->GetStreamingManager()->SetStrippedNumberOfLinesStreamingMode(200);
filter->SetInput(reader->GetOutput());
filter->Update();
diff --git a/Testing/Code/Common/CMakeLists.txt b/Testing/Code/Common/CMakeLists.txt
index a598f4ed523f02bf5972055c2e6eaaef5588f577..ef1cf514ef8833f72f36486c0297351017271db5 100644
--- a/Testing/Code/Common/CMakeLists.txt
+++ b/Testing/Code/Common/CMakeLists.txt
@@ -736,6 +736,19 @@ ADD_TEST(coTvImageRegionTileMapSplitter ${COMMON_TESTS8}
${TEMP}/coImageRegionTileMapSplitter.txt
)
+# ------------- otb::ImageRegionSquareTileSplitter ----------------------------
+ADD_TEST(coTuImageRegionSquareTileSplitterNew ${COMMON_TESTS8}
+ otbImageRegionSquareTileSplitterNew
+)
+
+ADD_TEST(coTvImageRegionSquareTileSplitter ${COMMON_TESTS8}
+ --compare-ascii ${NOTOL}
+ ${BASELINE_FILES}/coImageRegionSquareTileSplitter.txt
+ ${TEMP}/coImageRegionSquareTileSplitter.txt
+ otbImageRegionSquareTileSplitter
+ ${TEMP}/coImageRegionSquareTileSplitter.txt
+)
+
# ------------- otb::ImageOfVectorsToMonoChannelExtractROI ----------------------------
ADD_TEST(coTuImageOfVectorsToMonoChannelExtractROINew ${COMMON_TESTS8}
otbImageOfVectorsToMonoChannelExtractROINew
@@ -1103,6 +1116,7 @@ otbQuickLookImageGeneratorNew.cxx
otbQuickLookImageGenerator.cxx
otbImageRegionTileMapSplitterNew.cxx
otbImageRegionTileMapSplitter.cxx
+otbImageRegionSquareTileSplitter.cxx
otbImageOfVectorsToMonoChannelExtractROINew.cxx
otbImageOfVectorsToMonoChannelExtractROI.cxx
otbImageRegionNonUniformMultidimensionalSplitterNew.cxx
diff --git a/Testing/Code/Common/otbCommonTests8.cxx b/Testing/Code/Common/otbCommonTests8.cxx
index 63e946336f54ae8d97a12ce41cfe3478225f0214..c948600503a217f1a6f89e33618659d644aab84a 100644
--- a/Testing/Code/Common/otbCommonTests8.cxx
+++ b/Testing/Code/Common/otbCommonTests8.cxx
@@ -30,6 +30,8 @@ void RegisterTests()
REGISTER_TEST(otbQuickLookImageGenerator);
REGISTER_TEST(otbImageRegionTileMapSplitterNew);
REGISTER_TEST(otbImageRegionTileMapSplitter);
+ REGISTER_TEST(otbImageRegionSquareTileSplitterNew);
+ REGISTER_TEST(otbImageRegionSquareTileSplitter);
REGISTER_TEST(otbImageOfVectorsToMonoChannelExtractROINew);
REGISTER_TEST(otbImageOfVectorsToMonoChannelExtractROI);
REGISTER_TEST(otbImageRegionNonUniformMultidimensionalSplitterNew);
diff --git a/Testing/Code/Common/otbImageRegionSquareTileSplitter.cxx b/Testing/Code/Common/otbImageRegionSquareTileSplitter.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..89c790ebca52e0ceb456886d77d74b2a98bd1460
--- /dev/null
+++ b/Testing/Code/Common/otbImageRegionSquareTileSplitter.cxx
@@ -0,0 +1,146 @@
+/*=========================================================================
+
+ 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.
+
+=========================================================================*/
+#include "otbImageRegionSquareTileSplitter.h"
+#include <fstream>
+
+const int Dimension = 2;
+typedef otb::ImageRegionSquareTileSplitter<Dimension> SquareTileSplitterType;
+typedef SquareTileSplitterType::IndexType IndexType;
+typedef SquareTileSplitterType::SizeType SizeType;
+typedef SquareTileSplitterType::RegionType RegionType;
+
+
+int otbImageRegionSquareTileSplitterNew(int argc, char * argv[])
+{
+ SquareTileSplitterType::Pointer splitter = SquareTileSplitterType::New();
+
+ std::cout << splitter << std::endl;
+
+ return EXIT_SUCCESS;
+}
+
+int TestSplitter(const RegionType& region, unsigned int PixelSize, unsigned int MaxTileSize, std::ostream& os)
+{
+ os << "----------------------------------" << std::endl;
+ os << "Region : " << region << std::endl;
+ os << "PixelSize : " << PixelSize << std::endl;
+ os << "MaxTileSize : " << MaxTileSize << std::endl;
+
+ SquareTileSplitterType::Pointer splitter;
+ splitter = SquareTileSplitterType::New();
+
+ unsigned int requestedNbSplits = region.GetNumberOfPixels() * PixelSize / MaxTileSize;
+ if (requestedNbSplits == 0)
+ requestedNbSplits = 1;
+ os << "Requested Number of splits : " << requestedNbSplits << std::endl;
+
+ const unsigned int nbSplits = splitter->GetNumberOfSplits(region, requestedNbSplits);
+ os << "Actual Number of splits : " << nbSplits << std::endl;
+
+
+ RegionType split;
+
+ // First split :
+ split = splitter->GetSplit(0, nbSplits, region);
+ os << "First Split : " << split
+ << "(" << split.GetNumberOfPixels() * PixelSize << " bytes)" << std::endl;
+
+ if (nbSplits > 1)
+ {
+ // Second split :
+ split = splitter->GetSplit(1, nbSplits, region);
+ os << "Second Split : " << split
+ << "(" << split.GetNumberOfPixels() * PixelSize << " bytes)" << std::endl;
+ }
+
+ if (nbSplits > 2)
+ {
+ // Last split :
+ split = splitter->GetSplit(nbSplits - 1, nbSplits, region);
+ os << "Last Split : " << split
+ << "(" << split.GetNumberOfPixels() * PixelSize << " bytes)" << std::endl;
+ }
+
+ return EXIT_SUCCESS;
+}
+
+
+int otbImageRegionSquareTileSplitter(int argc, char * argv[])
+{
+ std::ofstream outfile(argv[1]);
+ RegionType region;
+
+ // Test with a 0-based indexed region
+ region.SetIndex(0, 0);
+ region.SetIndex(1, 0);
+ region.SetSize(0, 1024);
+ region.SetSize(1, 1024);
+ TestSplitter(region, 1, 128, outfile);
+ TestSplitter(region, 1, 512*512, outfile);
+ TestSplitter(region, 2, 512*512, outfile);
+ TestSplitter(region, 4, 512*512, outfile);
+ TestSplitter(region, 8, 512*512, outfile);
+
+ // Test with a shifted region
+ region.SetIndex(0, 42);
+ region.SetIndex(1, 42);
+ region.SetSize(0, 1000);
+ region.SetSize(1, 1000);
+ TestSplitter(region, 1, 128, outfile);
+ TestSplitter(region, 1, 512*512, outfile);
+ TestSplitter(region, 2, 512*512, outfile);
+ TestSplitter(region, 4, 512*512, outfile);
+ TestSplitter(region, 8, 512*512, outfile);
+
+ // Test with a negative shift
+ region.SetIndex(0, -42);
+ region.SetIndex(1, -42);
+ region.SetSize(0, 1000);
+ region.SetSize(1, 1000);
+ TestSplitter(region, 1, 128, outfile);
+ TestSplitter(region, 1, 512*512, outfile);
+ TestSplitter(region, 2, 512*512, outfile);
+ TestSplitter(region, 4, 512*512, outfile);
+ TestSplitter(region, 8, 512*512, outfile);
+
+ // Test with a reduced size
+ region.SetIndex(0, 0);
+ region.SetIndex(1, 0);
+ region.SetSize(0, 1);
+ region.SetSize(1, 1);
+ TestSplitter(region, 1, 128, outfile);
+ TestSplitter(region, 1, 512*512, outfile);
+ TestSplitter(region, 2, 512*512, outfile);
+ TestSplitter(region, 4, 512*512, outfile);
+ TestSplitter(region, 8, 512*512, outfile);
+
+ // Test with a reduced size, shifted
+ region.SetIndex(0, 42);
+ region.SetIndex(1, 42);
+ region.SetSize(0, 1);
+ region.SetSize(1, 1);
+ TestSplitter(region, 1, 128, outfile);
+ TestSplitter(region, 1, 512*512, outfile);
+ TestSplitter(region, 2, 512*512, outfile);
+ TestSplitter(region, 4, 512*512, outfile);
+ TestSplitter(region, 8, 512*512, outfile);
+
+ outfile.close();
+
+ return EXIT_SUCCESS;
+}
diff --git a/Testing/Code/Common/otbQuickLookImageGenerator.cxx b/Testing/Code/Common/otbQuickLookImageGenerator.cxx
index 071f7fd5c642a1bccebd3e308f9f994181a34261..8e86c247de19a0801e2741ae603c9c1805a68608 100644
--- a/Testing/Code/Common/otbQuickLookImageGenerator.cxx
+++ b/Testing/Code/Common/otbQuickLookImageGenerator.cxx
@@ -49,6 +49,8 @@ int otbQuickLookImageGenerator(int argc, char* argv[])
filter->SetMaximumKernelWidth(atoi(argv[6]));
filter->UseImageSpacing(atoi(argv[7]));
+ filter->Update();
+
writer->SetInput(filter->GetOutput());
writer->SetFileName(outputFileName);
diff --git a/Testing/Code/IO/otbImageFileWriterTestWithoutInput.cxx b/Testing/Code/IO/otbImageFileWriterTestWithoutInput.cxx
index d6576ab6e6be1de99c0fdb18cfc6ae756c807d08..6c9733dadd957382625aba51c521c9e133aae49b 100644
--- a/Testing/Code/IO/otbImageFileWriterTestWithoutInput.cxx
+++ b/Testing/Code/IO/otbImageFileWriterTestWithoutInput.cxx
@@ -70,9 +70,6 @@ int otbImageScalarFileWriterTestWithoutInputGeneric(int argc, char* argv[])
IteratorType it(image, image->GetLargestPossibleRegion());
it.GoToBegin();
- double val = 0.;
- ImagePixelType pixVal;
-
while (!it.IsAtEnd())
{
it.Set(static_cast<PixelType> (size[0] * it.GetIndex()[1] + it.GetIndex()[0]));
diff --git a/Testing/Code/IO/otbImageStreamingFileWriterTestWithoutInput.cxx b/Testing/Code/IO/otbImageStreamingFileWriterTestWithoutInput.cxx
index 372218fe4c2692489c0d638f1ff19d39c30f8507..cff8ed3b03c25fbcfde335d591029c4d29b541c4 100644
--- a/Testing/Code/IO/otbImageStreamingFileWriterTestWithoutInput.cxx
+++ b/Testing/Code/IO/otbImageStreamingFileWriterTestWithoutInput.cxx
@@ -71,9 +71,6 @@ int otbImageScalarStreamingFileWriterTestWithoutInputGeneric(int argc, char* arg
IteratorType it(image, image->GetLargestPossibleRegion());
it.GoToBegin();
- double val = 0.;
- ImagePixelType pixVal;
-
while (!it.IsAtEnd())
{
it.Set(static_cast<PixelType> (size[0] * it.GetIndex()[1] + it.GetIndex()[0]));
diff --git a/Testing/Code/IO/otbVectorImageFileWriterTestWithoutInput.cxx b/Testing/Code/IO/otbVectorImageFileWriterTestWithoutInput.cxx
index fbb4b31b2992fd0911f2ff9ab669c934721448e4..00ce431586e39b497718303d0fcc25d1c691d5b6 100644
--- a/Testing/Code/IO/otbVectorImageFileWriterTestWithoutInput.cxx
+++ b/Testing/Code/IO/otbVectorImageFileWriterTestWithoutInput.cxx
@@ -84,7 +84,6 @@ int otbVectorImageFileWriterScalarTestWithoutInputGeneric(int argc, char* argv[]
IteratorType it(image, image->GetLargestPossibleRegion());
it.GoToBegin();
- double val = 0.;
ImagePixelType pixVal;
pixVal.SetSize(atoi(argv[3]));
@@ -179,7 +178,6 @@ int otbVectorImageFileWriterComplexTestWithoutInputGeneric(int argc, char* argv[
IteratorType it(image, image->GetLargestPossibleRegion());
it.GoToBegin();
- double val = 0.;
ImagePixelType pixVal;
pixVal.SetSize(atoi(argv[3]));
diff --git a/Testing/Code/IO/otbVectorImageStreamingFileWriterTestWithoutInput.cxx b/Testing/Code/IO/otbVectorImageStreamingFileWriterTestWithoutInput.cxx
index b5768ec8717bca133d584a2605f5cb3da1750502..023d3c6adb4941ae625c5d63d45370e3f3bc1ba5 100644
--- a/Testing/Code/IO/otbVectorImageStreamingFileWriterTestWithoutInput.cxx
+++ b/Testing/Code/IO/otbVectorImageStreamingFileWriterTestWithoutInput.cxx
@@ -84,7 +84,6 @@ int otbVectorImageStreamingFileWriterScalarTestWithoutInputGeneric(int argc, cha
IteratorType it(image, image->GetLargestPossibleRegion());
it.GoToBegin();
- double val = 0.;
ImagePixelType pixVal;
pixVal.SetSize(atoi(argv[3]));
@@ -105,7 +104,6 @@ int otbVectorImageStreamingFileWriterScalarTestWithoutInputGeneric(int argc, cha
writer->SetInput(image);
writer->Update();
-
return EXIT_SUCCESS;
}
@@ -180,7 +178,6 @@ int otbVectorImageStreamingFileWriterComplexTestWithoutInputGeneric(int argc, ch
IteratorType it(image, image->GetLargestPossibleRegion());
it.GoToBegin();
- double val = 0.;
ImagePixelType pixVal;
pixVal.SetSize(atoi(argv[3]));
diff --git a/Testing/Code/ObjectDetection/otbDescriptorsListSampleGenerator.cxx b/Testing/Code/ObjectDetection/otbDescriptorsListSampleGenerator.cxx
index e52c587fc6da2f835a318a49af06f379b3915c7e..c4a7b9af956898cf608886499a3d2040dc6ccd92 100644
--- a/Testing/Code/ObjectDetection/otbDescriptorsListSampleGenerator.cxx
+++ b/Testing/Code/ObjectDetection/otbDescriptorsListSampleGenerator.cxx
@@ -157,15 +157,7 @@ int otbDescriptorsListSampleGenerator(int argc, char* argv[])
descriptorsGenerator->SetDescriptorsFunction(descriptorsFunction.GetPointer());
descriptorsGenerator->SetNeighborhoodRadius(neighborhood);
- if (streaming == 0)
- {
- descriptorsGenerator->GetStreamer()->SetNumberOfStreamDivisions(1);
- }
- else
- {
- descriptorsGenerator->GetStreamer()->SetNumberOfStreamDivisions(streaming);
- }
-
+ descriptorsGenerator->GetStreamer()->GetStreamingManager()->SetStrippedNumberOfLinesStreamingMode(streaming);
descriptorsGenerator->Update();
@@ -240,16 +232,7 @@ int otbDescriptorsSVMModelCreation(int argc, char* argv[])
descriptorsGenerator->SetSamplesLocations(vectorDataReader->GetOutput());
descriptorsGenerator->SetDescriptorsFunction(descriptorsFunction.GetPointer());
descriptorsGenerator->SetNeighborhoodRadius(5);
-
- if (streaming == 0)
- {
- descriptorsGenerator->GetStreamer()->SetNumberOfStreamDivisions(1);
- }
- else
- {
- descriptorsGenerator->GetStreamer()->SetNumberOfStreamDivisions(streaming);
- }
-
+ descriptorsGenerator->GetStreamer()->GetStreamingManager()->SetStrippedNumberOfLinesStreamingMode(streaming);
descriptorsGenerator->Update();
// Normalize the samples
diff --git a/Testing/Code/ObjectDetection/otbObjectDetectionClassifier.cxx b/Testing/Code/ObjectDetection/otbObjectDetectionClassifier.cxx
index fe3846df055eec4a684eaf7f6594ad4f6f2b7004..e394c2ea2b8c1114090f54b02b368c02a62f4cb9 100644
--- a/Testing/Code/ObjectDetection/otbObjectDetectionClassifier.cxx
+++ b/Testing/Code/ObjectDetection/otbObjectDetectionClassifier.cxx
@@ -135,15 +135,8 @@ int otbObjectDetectionClassifier(int argc, char* argv[])
classifier->SetShifts(statisticsReader->GetStatisticVectorByName("mean"));
classifier->SetScales(statisticsReader->GetStatisticVectorByName("stddev"));
- if (streaming == 0)
- {
- classifier->GetStreamer()->SetNumberOfStreamDivisions(1);
- }
- else
- {
- classifier->GetStreamer()->SetNumberOfStreamDivisions(streaming);
- }
-
+ // 0 means no streaming
+ classifier->GetStreamer()->GetStreamingManager()->SetStrippedNumberOfLinesStreamingMode(streaming);
classifier->Update();
std::vector<ObjectDetectionClassifierType::PointType> points;
diff --git a/Testing/Code/Projections/CMakeLists.txt b/Testing/Code/Projections/CMakeLists.txt
index 7b94e96cfdcadda6f49a1f88b1498dfcc33336ac..c2d6ebdaf1d3bdcf2c5daf7f0fe4e28ce4fb6c24 100644
--- a/Testing/Code/Projections/CMakeLists.txt
+++ b/Testing/Code/Projections/CMakeLists.txt
@@ -842,6 +842,20 @@ ADD_TEST(prTvGenericRSTransformFromImage ${PROJECTIONS_TESTS3}
otbGenericRSTransformFromImage
${INPUTDATA}/WithoutProjRefWithKeywordlist.tif)
+IF(OTB_DATA_USE_LARGEINPUT)
+ADD_TEST(prTvGenericRSTransformQuickbirdToulouseGeodesicPointChecking ${PROJECTIONS_TESTS3}
+otbGenericRSTransformImageAndMNTToWGS84ConversionChecking
+${LARGEINPUT}/QUICKBIRD/TOULOUSE/000000128955_01_P001_PAN/02APR01105228-P1BS-000000128955_01_P001.TIF
+${INPUTDATA}/DEM/srtm_directory/
+16271 15647 # IGN reference point in front of st sernin church pointed in image
+1.48353 43.55968 # IGN reference point lon/lat
+192.205 # IGN reference point elevation (prec < 10 cm)
+5
+0.00001
+)
+
+ENDIF(OTB_DATA_USE_LARGEINPUT)
+
ADD_TEST(prTuVectorDataProjectionFilterNew ${PROJECTIONS_TESTS3} otbVectorDataProjectionFilterNew )
#test points
diff --git a/Testing/Code/Projections/otbGenericRSTransformFromImage.cxx b/Testing/Code/Projections/otbGenericRSTransformFromImage.cxx
index 415d40aefd1fec2b73dfabd8e497cdbacd9c4ba5..2386cea022af4756664f8cff3b93cd4eb1db524c 100644
--- a/Testing/Code/Projections/otbGenericRSTransformFromImage.cxx
+++ b/Testing/Code/Projections/otbGenericRSTransformFromImage.cxx
@@ -24,11 +24,14 @@
#include "otbGenericRSTransform.h"
#include <ogr_spatialref.h>
#include <fstream>
+#include "itkEuclideanDistance.h"
+
typedef otb::VectorImage<double, 2> ImageType;
typedef otb::ImageFileReader<ImageType> ReaderType;
typedef otb::GenericRSTransform<> TransformType;
typedef TransformType::InputPointType PointType;
+typedef itk::Statistics::EuclideanDistance<ImageType::PointType> DistanceType;
int otbGenericRSTransformFromImage(int argc, char* argv[])
{
@@ -64,3 +67,125 @@ int otbGenericRSTransformFromImage(int argc, char* argv[])
return EXIT_SUCCESS;
}
+
+int otbGenericRSTransformImageAndMNTToWGS84ConversionChecking(int argc, char* argv[])
+{
+ std::string infname = argv[1];
+ std::string demdir = argv[2];
+
+ ImageType::PointType refImgPt, refGeoPt, estimatedImgPt, estimatedGeoPt;
+ refImgPt[0] = atof(argv[3]);
+ refImgPt[1] = atof(argv[4]);
+ refGeoPt[0] = atof(argv[5]);
+ refGeoPt[1] = atof(argv[6]);
+ double refHeight = atof(argv[7]);
+ double imgTol = atof(argv[8]);
+ double geoTol = atof(argv[9]);
+
+ bool pass = true;
+
+// Reader
+ ReaderType::Pointer reader = ReaderType::New();
+ reader->SetFileName(argv[1]);
+ reader->UpdateOutputInformation();
+
+ // Build wgs ref
+ OGRSpatialReference oSRS;
+ oSRS.SetWellKnownGeogCS("WGS84");
+ char * wgsRef = NULL;
+ oSRS.exportToWkt(&wgsRef);
+
+ DistanceType::Pointer distance = DistanceType::New();
+
+ // Instanciate WGS->Image transform
+ TransformType::Pointer wgs2img = TransformType::New();
+ wgs2img->SetInputProjectionRef(wgsRef);
+ wgs2img->SetOutputProjectionRef(reader->GetOutput()->GetProjectionRef());
+ wgs2img->SetOutputKeywordList(reader->GetOutput()->GetImageKeywordlist());
+ wgs2img->SetDEMDirectory(argv[2]);
+ wgs2img->InstanciateTransform();
+
+ // Instanciate Image->WGS transform
+ TransformType::Pointer img2wgs = TransformType::New();
+ img2wgs->SetInputProjectionRef(reader->GetOutput()->GetProjectionRef());
+ img2wgs->SetInputKeywordList(reader->GetOutput()->GetImageKeywordlist());
+ img2wgs->SetOutputProjectionRef(wgsRef);
+ img2wgs->SetDEMDirectory(argv[2]);
+ img2wgs->InstanciateTransform();
+
+ std::cout<< std::setprecision(8) << std::endl;
+
+ std::cout<<"References: "<<std::endl;
+ std::cout<<"Geo (wgs84): "<<refGeoPt<<std::endl;
+ std::cout<<"Image: "<<refImgPt<<std::endl;
+ std::cout<<"Height: "<<refHeight<<std::endl;
+
+ std::cout<<"Using elevation from "<<demdir<<std::endl;
+
+ estimatedImgPt = wgs2img->TransformPoint(refGeoPt);
+ estimatedGeoPt = img2wgs->TransformPoint(refImgPt);
+
+ std::cout<<"Forward(refImg): "<<refImgPt<<" -> " << estimatedGeoPt<<std::endl;
+ std::cout<<"Inverse(refGeo): "<<refGeoPt<<" -> "<< estimatedImgPt <<std::endl;
+
+ double imgRes = distance->Evaluate(refImgPt, estimatedImgPt);
+ double geoRes = distance->Evaluate(refGeoPt, estimatedGeoPt);
+
+ if(imgRes > imgTol)
+ {
+ pass = false;
+ std::cerr<<"Image residual is too high: "<<imgRes<<std::endl;
+ }
+
+ if(geoRes > geoTol)
+ {
+ pass = false;
+ std::cerr<<"Geographic (WGS84) residual is too high: "<<geoRes<<std::endl;
+ }
+
+ std::cout<<"Using reference elevation: "<<std::endl;
+ wgs2img = TransformType::New();
+ wgs2img->SetInputProjectionRef(wgsRef);
+ wgs2img->SetOutputProjectionRef(reader->GetOutput()->GetProjectionRef());
+ wgs2img->SetOutputKeywordList(reader->GetOutput()->GetImageKeywordlist());
+ wgs2img->SetAverageElevation(refHeight);
+ wgs2img->InstanciateTransform();
+
+ img2wgs = TransformType::New();
+ img2wgs->SetInputProjectionRef(reader->GetOutput()->GetProjectionRef());
+ img2wgs->SetInputKeywordList(reader->GetOutput()->GetImageKeywordlist());
+ img2wgs->SetOutputProjectionRef(wgsRef);
+ img2wgs->SetAverageElevation(refHeight);
+ img2wgs->InstanciateTransform();
+
+ estimatedImgPt = wgs2img->TransformPoint(refGeoPt);
+ estimatedGeoPt = img2wgs->TransformPoint(refImgPt);
+
+ std::cout<<"Forward(refImg): "<<refImgPt<<" -> " << estimatedGeoPt<<std::endl;
+ std::cout<<"Inverse(refGeo): "<<refGeoPt<<" -> "<< estimatedImgPt <<std::endl;
+
+ imgRes = distance->Evaluate(refImgPt, estimatedImgPt);
+ geoRes = distance->Evaluate(refGeoPt, estimatedGeoPt);
+
+ if(imgRes > imgTol)
+ {
+ pass = false;
+ std::cerr<<"Image residual is too high: "<<imgRes<<std::endl;
+ }
+
+ if(geoRes > geoTol)
+ {
+ pass = false;
+ std::cerr<<"Geographic (WGS84) residual is too high: "<<geoRes<<std::endl;
+ }
+
+ if(pass)
+ {
+ return EXIT_SUCCESS;
+ }
+ else
+ {
+ std::cerr<<"There were imprecise results."<<std::endl;
+ return EXIT_FAILURE;
+ }
+}
diff --git a/Testing/Code/Projections/otbProjectionsTests3.cxx b/Testing/Code/Projections/otbProjectionsTests3.cxx
index cfac112ba3e9ca38fab26036109da7b343ddca08..418a62ad10633f6427417c0bdad453eb32c3b770 100644
--- a/Testing/Code/Projections/otbProjectionsTests3.cxx
+++ b/Testing/Code/Projections/otbProjectionsTests3.cxx
@@ -34,6 +34,7 @@ void RegisterTests()
REGISTER_TEST(otbGenericRSTransform);
REGISTER_TEST(otbGenericRSTransformWithSRID);
REGISTER_TEST(otbGenericRSTransformFromImage);
+ REGISTER_TEST(otbGenericRSTransformImageAndMNTToWGS84ConversionChecking);
REGISTER_TEST(otbVectorDataProjectionFilterNew);
REGISTER_TEST(otbVectorDataProjectionFilter);
REGISTER_TEST(otbVectorDataProjectionFilterFromMapToSensor);
diff --git a/Testing/Code/SARPolarimetry/CMakeLists.txt b/Testing/Code/SARPolarimetry/CMakeLists.txt
index 654b91338eab7045ea798de4dbf7a08b944173ee..c291daf45d2e66c98f55facf982901265abd93a8 100644
--- a/Testing/Code/SARPolarimetry/CMakeLists.txt
+++ b/Testing/Code/SARPolarimetry/CMakeLists.txt
@@ -373,10 +373,6 @@ ADD_TEST(saTvMuellerToCovarianceImageFilter ${SARPOLARIMETRY_TESTS2}
${BASELINE}/saTvSinclairImageFilter_SinclairToMueller.tif
${TEMP}/saTvMuellerToMLCImageFilter.tif
)
-# Hermitian eigen analysis class
-ADD_TEST(saTvHermitianEigenAnalysisTest ${SARPOLARIMETRY_TESTS2}
- otbHermitianEigenAnalysisTest
-)
# Polarimetric Data
ADD_TEST(saTuPolarimetricDataNew ${SARPOLARIMETRY_TESTS2}
@@ -430,7 +426,6 @@ otbMuellerToPolarisationDegreeAndPowerImageFilterNew.cxx
otbMuellerToPolarisationDegreeAndPowerImageFilter.cxx
otbMuellerToCovarianceImageFilterNew.cxx
otbMuellerToCovarianceImageFilter.cxx
-otbHermitianEigenAnalysisTest.cxx
otbPolarimetricData.cxx
)
diff --git a/Testing/Code/SARPolarimetry/otbHermitianEigenAnalysisTest.cxx b/Testing/Code/SARPolarimetry/otbHermitianEigenAnalysisTest.cxx
index d21ca4c3ddde6269fadc9742e7d8606c7d56ee1b..6ef54da7429f59bcfa0037dba639394e1f16b906 100644
--- a/Testing/Code/SARPolarimetry/otbHermitianEigenAnalysisTest.cxx
+++ b/Testing/Code/SARPolarimetry/otbHermitianEigenAnalysisTest.cxx
@@ -26,6 +26,8 @@
#include "otbVectorImage.h"
#include "otbHermitianEigenAnalysis.h"
#include "itkVariableLengthVector.h"
+#include "vnl/algo/vnl_complex_eigensystem.h"
+#include <complex>
int otbHermitianEigenAnalysisTest(int argc, char * argv[])
@@ -37,6 +39,7 @@ int otbHermitianEigenAnalysisTest(int argc, char * argv[])
typedef otb::HermitianEigenAnalysis<MatrixType, EigenvalueType, EigenMatrixType> FilterType;
EigenMatrixType resEigVal;
+
itk::Vector<float, 6> temp;
temp[0] = 0.162793;
temp[1] = -0.432753;
@@ -115,5 +118,43 @@ int otbHermitianEigenAnalysisTest(int argc, char * argv[])
return EXIT_FAILURE;
}
+ std::cout<<"Values: "<<vect<<std::endl;
+ std::cout<<"Vectors: "<<eigVal<<std::endl;
+
+
+ typedef std::complex<double> ComplexType;
+ vnl_matrix<ComplexType> vnlMat(3, 3, 0);
+ vnlMat[0][0] = ComplexType(0., 0.);
+ vnlMat[0][1] = ComplexType(1.5, 2.);
+ vnlMat[0][2] = ComplexType(2.5, 3.);
+ vnlMat[1][0] = ComplexType(1.5, -2.);
+ vnlMat[1][1] = ComplexType(0.5, 0.);
+ vnlMat[1][2] = ComplexType(3.5, 4.);
+ vnlMat[2][0] = ComplexType(2.5, -3.);
+ vnlMat[2][1] = ComplexType(3.5, -4.);
+ vnlMat[2][2] = ComplexType(1., 0.);
+
+ std::cout<<"Matrix:: "<<std::endl;
+ vnlMat.print(std::cout);
+
+ vnl_complex_eigensystem syst(vnlMat, true, true);
+
+ vnl_matrix< ComplexType > pm = syst.L;
+ std::cout<<"Left:: "<<std::endl;
+ pm.print(std::cout);
+ std::cout<<"Right:: "<<std::endl;
+
+
+ pm = syst.R;
+ pm.print(std::cout);
+ std::cout<<"W:: "<<std::endl;
+ vnl_vector< ComplexType > lol = syst.W;
+
+ for(unsigned i=0; i<lol.size(); i++)
+ {
+ std::cout<<" "<<lol[i];
+ }
+ std::cout<<std::endl;
+
return EXIT_SUCCESS;
}
diff --git a/Testing/Code/SARPolarimetry/otbReciprocalHAlphaImageFilter.cxx b/Testing/Code/SARPolarimetry/otbReciprocalHAlphaImageFilter.cxx
index ed970d31910f298a5b56d58a236f42131a5effd2..79c312c5cbeaf8059a450909ebb6e24aa2c14a17 100644
--- a/Testing/Code/SARPolarimetry/otbReciprocalHAlphaImageFilter.cxx
+++ b/Testing/Code/SARPolarimetry/otbReciprocalHAlphaImageFilter.cxx
@@ -89,9 +89,10 @@ int otbReciprocalHAlphaImageFilter(int argc, char * argv[])
PerBandMeanFilterType::Pointer perBandMeanFilter = PerBandMeanFilterType::New();
perBandMeanFilter->SetInput(sinclairToCoherencyFilter->GetOutput());
-
+
FilterType::Pointer filter = FilterType::New();
filter->SetInput(perBandMeanFilter->GetOutput());
+ filter->SetNumberOfThreads(1);
ExtractType::Pointer extract = ExtractType::New();
extract->SetInput(filter->GetOutput());
@@ -102,6 +103,8 @@ int otbReciprocalHAlphaImageFilter(int argc, char * argv[])
writer->SetFileName(outputFilename);
writer->SetInput(extract->GetOutput());
+ writer->SetNumberOfThreads(1);
+writer->SetNumberOfStreamDivisions(1);
writer->Update();
return EXIT_SUCCESS;
diff --git a/Testing/Code/SARPolarimetry/otbSARPolarimetryTests2.cxx b/Testing/Code/SARPolarimetry/otbSARPolarimetryTests2.cxx
index bdda5c32405263bb9f8f4d1bf94de291a20a6d3d..97a622061338743fcbf26f51075f6c820ca6d614 100644
--- a/Testing/Code/SARPolarimetry/otbSARPolarimetryTests2.cxx
+++ b/Testing/Code/SARPolarimetry/otbSARPolarimetryTests2.cxx
@@ -42,7 +42,6 @@ void RegisterTests()
REGISTER_TEST(otbMuellerToPolarisationDegreeAndPowerImageFilter);
REGISTER_TEST(otbMuellerToCovarianceImageFilterNew);
REGISTER_TEST(otbMuellerToCovarianceImageFilter);
- REGISTER_TEST(otbHermitianEigenAnalysisTest);
REGISTER_TEST(otbPolarimetricDataNew);
REGISTER_TEST(otbPolarimetricDataTest);
}
diff --git a/Testing/Code/Visualization/otbImageLayerScalar.cxx b/Testing/Code/Visualization/otbImageLayerScalar.cxx
index d14ac1518b77547d8c0632cd6bc211dbe82fe467..515e5e8dd2209a42cf008fc5bc35693202113bc7 100644
--- a/Testing/Code/Visualization/otbImageLayerScalar.cxx
+++ b/Testing/Code/Visualization/otbImageLayerScalar.cxx
@@ -53,6 +53,7 @@ int otbImageLayerScalar(int argc, char * argv[])
// Quicklook
shrinker->SetInput(reader->GetOutput());
shrinker->SetShrinkFactor(ssrate);
+ shrinker->Update();
// new layer
LayerType::Pointer layer = LayerType::New();
diff --git a/Testing/Code/Visualization/otbImageLayerVector.cxx b/Testing/Code/Visualization/otbImageLayerVector.cxx
index 04ab0b7ca10fd739ef33b9336cdf88fb10a9ad8d..94e46a11af24b5373a0fc062b236dc7152119f13 100644
--- a/Testing/Code/Visualization/otbImageLayerVector.cxx
+++ b/Testing/Code/Visualization/otbImageLayerVector.cxx
@@ -65,6 +65,7 @@ int otbImageLayerVector(int argc, char * argv[])
// Quicklook
shrinker->SetInput(reader->GetOutput());
shrinker->SetShrinkFactor(ssrate);
+ shrinker->Update();
// new layer
LayerType::Pointer layer = LayerType::New();
diff --git a/otbConfigure.h.in b/otbConfigure.h.in
index ab1f7d425524181d86d0022f76f02bbb8822b9be..29e21a97c10853e99033fb7e2de510ead5f3b417 100644
--- a/otbConfigure.h.in
+++ b/otbConfigure.h.in
@@ -1,5 +1,5 @@
/*
- * here is where system comupted values get stored these values should only
+ * here is where system computed values get stored these values should only
* change when the target compile platform changes
*/