From 043290a42f129bfadff693a59dc6b6ea5a78d89b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ga=C3=ABlle=20USSEGLIO?= <gaelle.usseglio@cnes.fr>
Date: Sat, 11 Apr 2020 14:59:11 +0000
Subject: [PATCH] BUG : First correction for CorrelationRough (streaming)
---
app/otbSARCorrelationRough.cxx | 47 +-
...otbSARStreamingMaximumMinimumImageFilter.h | 326 ++++++++++++++
...bSARStreamingMaximumMinimumImageFilter.txx | 415 ++++++++++++++++++
3 files changed, 779 insertions(+), 9 deletions(-)
create mode 100644 include/otbSARStreamingMaximumMinimumImageFilter.h
create mode 100644 include/otbSARStreamingMaximumMinimumImageFilter.txx
diff --git a/app/otbSARCorrelationRough.cxx b/app/otbSARCorrelationRough.cxx
index 15e767f..5d9f4fa 100644
--- a/app/otbSARCorrelationRough.cxx
+++ b/app/otbSARCorrelationRough.cxx
@@ -24,6 +24,8 @@
#include "itkFFTNormalizedCorrelationImageFilter.h"
#include "itkMinimumMaximumImageCalculator.h"
+#include "otbSARStreamingMaximumMinimumImageFilter.h"
+
#include <iostream>
#include <string>
#include <fstream>
@@ -42,9 +44,10 @@ public:
itkNewMacro(Self);
itkTypeMacro(SARCorrelationRough, otb::Wrapper::Application);
- // Filters
+// Filters
typedef itk::FFTNormalizedCorrelationImageFilter<FloatImageType, FloatImageType> CorFilterType;
typedef itk::MinimumMaximumImageCalculator< FloatImageType> MinMaxCalculatorType;
+typedef otb::SARStreamingMaximumMinimumImageFilter< FloatImageType> MinMaxFilterType;
private:
void DoInit() override
@@ -124,6 +127,8 @@ void DoExecute() override
otbAppLogINFO(<<"ML Factor on azimut :"<<factorML_azi);
otbAppLogINFO(<<"ML Factor on range : "<<factorML_ran);
+ otbAppLogINFO(<<"RAM : "<< GetParameterInt("ram"));
+
// Get master and slave image
FloatImageType::Pointer MasterPtr = GetParameterFloatImage("inmaster");
FloatImageType::Pointer SlavePtr = GetParameterFloatImage("inslave");
@@ -134,22 +139,46 @@ void DoExecute() override
m_Ref.push_back(correlationFilter.GetPointer());
correlationFilter->SetFixedImage(MasterPtr);
correlationFilter->SetMovingImage(SlavePtr);
- correlationFilter->Update();
+ //correlationFilter->Update();
// Min/Max calculator
- MinMaxCalculatorType::Pointer minMaxFilter = MinMaxCalculatorType::New();
- minMaxFilter->SetImage(correlationFilter->GetOutput());
- minMaxFilter->ComputeMaximum();
+ // MinMaxCalculatorType::Pointer minMaxFilter = MinMaxCalculatorType::New();
+ // minMaxFilter->SetImage(correlationFilter->GetOutput());
+ // minMaxFilter->ComputeMaximum();
+
+
+
+ MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New();
+ minMaxFilter->SetInput(correlationFilter->GetOutput());
+ // Adapt streaming with ram parameter (default 256 MB)
+ minMaxFilter->GetStreamer()->SetAutomaticStrippedStreaming(GetParameterInt("ram"));
+ minMaxFilter->Update();
+
+
+ std::cout << "minMaxFilter->GetIndexOfMax() : " << minMaxFilter->GetIndexOfMax() << std::endl;
+ std::cout << "minMaxFilter->GetMax() : " << minMaxFilter->GetMax() << std::endl;
+
+
+ // float shiftSLC_range = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/2.)-
+ // minMaxFilter->GetIndexOfMaximum()[0]) * static_cast<float>(factorML_ran);
+ // float shiftSLC_azimut = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/2.)
+ // -minMaxFilter->GetIndexOfMaximum()[1]) * static_cast<float>(factorML_azi);
+
+ // float shiftML_range = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/2.)-
+ // minMaxFilter->GetIndexOfMaximum()[0]);
+ // float shiftML_azimut = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/2.)
+ // -minMaxFilter->GetIndexOfMaximum()[1]);
+
float shiftSLC_range = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/2.)-
- minMaxFilter->GetIndexOfMaximum()[0]) * static_cast<float>(factorML_ran);
+ minMaxFilter->GetIndexOfMax()[0]) * static_cast<float>(factorML_ran);
float shiftSLC_azimut = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/2.)
- -minMaxFilter->GetIndexOfMaximum()[1]) * static_cast<float>(factorML_azi);
+ -minMaxFilter->GetIndexOfMax()[1]) * static_cast<float>(factorML_azi);
float shiftML_range = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/2.)-
- minMaxFilter->GetIndexOfMaximum()[0]);
+ minMaxFilter->GetIndexOfMax()[0]);
float shiftML_azimut = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/2.)
- -minMaxFilter->GetIndexOfMaximum()[1]);
+ -minMaxFilter->GetIndexOfMax()[1]);
// Assigne Outputs
SetParameterFloat("shiftranslc", shiftSLC_range);
diff --git a/include/otbSARStreamingMaximumMinimumImageFilter.h b/include/otbSARStreamingMaximumMinimumImageFilter.h
new file mode 100644
index 0000000..13861e6
--- /dev/null
+++ b/include/otbSARStreamingMaximumMinimumImageFilter.h
@@ -0,0 +1,326 @@
+/*
+ * Copyright (C) 2005-2018 Centre National d'Etudes Spatiales (CNES)
+ *
+ * This file is part of Orfeo Toolbox
+ *
+ * https://www.orfeo-toolbox.org/
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef otbSARStreamingMaximumMinimumImageFilter_h
+#define otbSARStreamingMaximumMinimumImageFilter_h
+
+#include "otbPersistentImageFilter.h"
+#include "itkNumericTraits.h"
+#include "itkArray.h"
+#include "itkSimpleDataObjectDecorator.h"
+#include "otbPersistentFilterStreamingDecorator.h"
+
+#include <complex>
+#include <cmath>
+
+namespace otb
+{
+
+/** \class PersistentMaximumMinimumImageFilter
+ * \brief Estimate some values for an input image (min/max and corresponding indexes).
+ *
+ * This filter persists its temporary data. It means that if you Update it n times on n different
+ * requested regions, the output statistics will be the statitics of the whole set of n regions.
+ *
+ * To reset the temporary data, one should call the Reset() function.
+ *
+ * To get the wanted values once the regions have been processed via the pipeline, use the Synthetize() method.
+ *
+ * \ingroup DiapOTBModule
+ */
+template<class TInputImage>
+class ITK_EXPORT PersistentMaximumMinimumImageFilter :
+ public PersistentImageFilter<TInputImage, TInputImage>
+{
+public:
+ /** Standard Self typedef */
+ typedef PersistentMaximumMinimumImageFilter Self;
+ typedef PersistentImageFilter<TInputImage, 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(PersistentMaximumMinimumImageFilter, PersistentImageFilter);
+
+ /** Image related typedefs. For Image Input*/
+ typedef TInputImage ImageType;
+ typedef typename TInputImage::Pointer InputImagePointer;
+
+ typedef typename TInputImage::RegionType RegionType;
+ typedef typename TInputImage::SizeType SizeType;
+ typedef typename TInputImage::IndexType IndexType;
+ typedef typename TInputImage::PixelType PixelType;
+ typedef double ValueType;
+ typedef typename ImageType::IndexValueType IndexValueType;
+ typedef typename ImageType::SizeValueType SizeValueType;
+
+ itkStaticConstMacro(InputImageDimension, unsigned int,
+ TInputImage::ImageDimension);
+
+ /** Image related typedefs. */
+ itkStaticConstMacro(ImageDimension, unsigned int,
+ TInputImage::ImageDimension);
+
+ /** Type to use for computations. */
+ // typedef typename itk::NumericTraits<PixelType>::RealType RealType;
+ typedef typename itk::NumericTraits<double>::RealType RealType;
+
+ /** Smart Pointer type to a DataObject. */
+ typedef typename itk::DataObject::Pointer DataObjectPointer;
+ typedef itk::ProcessObject::DataObjectPointerArraySizeType DataObjectPointerArraySizeType;
+
+ /** Type of DataObjects used for scalar outputs */
+ typedef itk::SimpleDataObjectDecorator<RealType> RealObjectType;
+ typedef itk::SimpleDataObjectDecorator<long> LongObjectType;
+ typedef itk::SimpleDataObjectDecorator<unsigned int> UintObjectType;
+ typedef itk::SimpleDataObjectDecorator<PixelType> PixelObjectType;
+ typedef itk::SimpleDataObjectDecorator<ValueType> ValueObjectType;
+ typedef itk::SimpleDataObjectDecorator<IndexType> IndexObjectType;
+
+ /** Return the number of valid points contained into the input grid. */
+ IndexType GetIndexOfMin() const
+ {
+ return this->GetIndexOfMinOutput()->Get();
+ }
+ IndexObjectType* GetIndexOfMinOutput();
+ const IndexObjectType* GetIndexOfMinOutput() const;
+
+ IndexType GetIndexOfMax() const
+ {
+ return this->GetIndexOfMaxOutput()->Get();
+ }
+ IndexObjectType* GetIndexOfMaxOutput();
+ const IndexObjectType* GetIndexOfMaxOutput() const;
+
+ /** Return the Minimum. */
+ ValueType GetMin() const
+ {
+ return this->GetMinOutput()->Get();
+ }
+ ValueObjectType* GetMinOutput();
+ const ValueObjectType* GetMinOutput() const;
+
+ /** Return the Maximum. */
+ ValueType GetMax() const
+ {
+ return this->GetMaxOutput()->Get();
+ }
+ ValueObjectType* GetMaxOutput();
+ const ValueObjectType* GetMaxOutput() const;
+
+ /** Make a DataObject of the correct type to be used as the specified
+ * output. */
+ DataObjectPointer MakeOutput(DataObjectPointerArraySizeType idx) ITK_OVERRIDE;
+ using Superclass::MakeOutput;
+
+ /**
+ * Synthetize and Reset function called by our PersistentFilterStreamingDecorator
+ */
+ void Synthetize(void) ITK_OVERRIDE;
+ void Reset(void) ITK_OVERRIDE;
+
+ itkSetMacro(IgnoreInfiniteValues, bool);
+ itkGetMacro(IgnoreInfiniteValues, bool);
+
+ itkSetMacro(IgnoreUserDefinedValue, bool);
+ itkGetMacro(IgnoreUserDefinedValue, bool);
+
+ itkSetMacro(UserIgnoredValue, RealType);
+ itkGetMacro(UserIgnoredValue, RealType);
+
+protected:
+ PersistentMaximumMinimumImageFilter();
+ ~PersistentMaximumMinimumImageFilter() ITK_OVERRIDE {}
+ void PrintSelf(std::ostream& os, itk::Indent indent) const ITK_OVERRIDE;
+
+
+ /** PersistentMaximumMinimumImageFilter needs a larger input requested region than the output
+ * requested region. As such, SARQuadraticAveragingImageFilter needs to provide an
+ * implementation for GenerateInputRequestedRegion() in order to inform the
+ * pipeline execution model.
+ * \sa ProcessObject::GenerateInputRequestedRegion() */
+ virtual void GenerateInputRequestedRegion() ITK_OVERRIDE;
+
+ /**
+ * OutputRegionToInputRegion returns the input region
+ */
+ RegionType OutputRegionToInputRegion(const RegionType& outputRegion) const;
+
+ /** Multi-thread version GenerateData. */
+ void ThreadedGenerateData(const RegionType&
+ outputRegionForThread,
+ itk::ThreadIdType threadId) ITK_OVERRIDE;
+
+ /** Pass the input through unmodified. Do this by Grafting in the
+ * AllocateOutputs method.
+ */
+ void AllocateOutputs() ITK_OVERRIDE;
+ void GenerateOutputInformation() ITK_OVERRIDE;
+
+private:
+ PersistentMaximumMinimumImageFilter(const Self &); //purposely not implemented
+ void operator =(const Self&); //purposely not implemented
+
+ itk::Array<unsigned int> m_ThreadIndexMaxC;
+ itk::Array<unsigned int> m_ThreadIndexMaxL;
+ itk::Array<unsigned int> m_ThreadIndexMinC;
+ itk::Array<unsigned int> m_ThreadIndexMinL;
+ itk::Array<ValueType> m_ThreadMax;
+ itk::Array<ValueType> m_ThreadMin;
+;
+
+ /* Ignored values */
+ bool m_IgnoreInfiniteValues;
+ bool m_IgnoreUserDefinedValue;
+ RealType m_UserIgnoredValue;
+ std::vector<unsigned long> m_IgnoredInfinitePixelCount;
+ std::vector<unsigned int> m_IgnoredUserPixelCount;
+
+
+}; // end of class PersistentMaximumMinimumImageFilter
+
+/*===========================================================================*/
+
+/** \class SARStreamingMaximumMinimumImageFilter
+ * \brief This class streams the whole input image through the PersistentMaximumMinimumImageFilter.
+ *
+ * This way, it allows computing the first order global statistics of this image. It calls the
+ * Reset() method of the PersistentMaximumMinimumImageFilter before streaming the image and the
+ * Synthetize() method of the PersistentMaximumMinimumImageFilter after having streamed the image
+ * to compute the statistics. The accessor on the results are wrapping the accessors of the
+ * internal PersistentMaximumMinimumImageFilter.
+ * By default infinite values are ignored, use IgnoreInfiniteValues accessor to consider
+ * infinite values in the computation.
+ *
+ * This filter can be used as:
+ * \code
+ * typedef otb::StreamingMaximumMinimumImageFilter<ImageType> StatisticsType;
+ * StatisticsType::Pointer statistics = StatisticsType::New();
+ * statistics->SetInput(reader->GetOutput());
+ * statistics->Update();
+ * std::cout << statistics-> GetIndexOfMax() << std::endl;
+ * std::cout << statistics-> GetMax() << std::endl;
+ * \endcode
+ *
+ * \sa PersistentMaximumMinimumImageFilter
+ * \sa PersistentImageFilter
+ * \sa PersistentFilterStreamingDecorator
+ * \sa StreamingImageVirtualWriter
+ * \ingroup Streamed
+ * \ingroup Multithreaded
+ * \ingroup MathematicalStatisticsImageFilters
+ *
+ * \ingroup DiapOTBModule
+ */
+
+template<class TInputImage>
+class ITK_EXPORT SARStreamingMaximumMinimumImageFilter :
+ public PersistentFilterStreamingDecorator<PersistentMaximumMinimumImageFilter<TInputImage> >
+{
+public:
+ /** Standard Self typedef */
+ typedef SARStreamingMaximumMinimumImageFilter Self;
+ typedef PersistentFilterStreamingDecorator
+ <PersistentMaximumMinimumImageFilter<TInputImage> > Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Type macro */
+ itkNewMacro(Self);
+
+ /** Creation through object factory macro */
+ itkTypeMacro(StreamingMaximumMinimumImageFilter, PersistentFilterStreamingDecorator);
+
+ typedef typename Superclass::FilterType StatFilterType;
+ typedef typename StatFilterType::PixelType PixelType;
+ typedef double ValueType;
+ typedef typename StatFilterType::RealType RealType;
+ typedef TInputImage InputImageType;
+ typedef typename TInputImage::IndexType IndexType;
+
+ /** Type of DataObjects used for scalar outputs */
+ typedef itk::SimpleDataObjectDecorator<RealType> RealObjectType;
+ typedef itk::SimpleDataObjectDecorator<long> LongObjectType;
+ typedef itk::SimpleDataObjectDecorator<PixelType> PixelObjectType;
+ typedef itk::SimpleDataObjectDecorator<ValueType> ValueObjectType;
+ typedef itk::SimpleDataObjectDecorator<IndexType> IndexObjectType;
+
+ using Superclass::SetInput;
+ void SetInput(InputImageType * input)
+ {
+ this->GetFilter()->SetInput(input);
+ }
+ const InputImageType * GetInput()
+ {
+ return this->GetFilter()->GetInput();
+ }
+
+ IndexType GetIndexOfMax()
+ {
+ return this->GetFilter()->GetIndexOfMaxOutput()->Get();
+ }
+
+ IndexType GetIndexOfMin()
+ {
+ return this->GetFilter()->GetIndexOfMinOutput()->Get();
+ }
+
+
+ ValueType GetMax()
+ {
+ return this->GetFilter()->GetMaxOutput()->Get();
+ }
+
+ ValueType GetMin()
+ {
+ return this->GetFilter()->GetMinOutput()->Get();
+ }
+
+
+ otbSetObjectMemberMacro(Filter, IgnoreInfiniteValues, bool);
+ otbGetObjectMemberMacro(Filter, IgnoreInfiniteValues, bool);
+
+ otbSetObjectMemberMacro(Filter, IgnoreUserDefinedValue, bool);
+ otbGetObjectMemberMacro(Filter, IgnoreUserDefinedValue, bool);
+
+ otbSetObjectMemberMacro(Filter, UserIgnoredValue, RealType);
+ otbGetObjectMemberMacro(Filter, UserIgnoredValue, RealType);
+
+protected:
+ /** Constructor */
+ SARStreamingMaximumMinimumImageFilter() {};
+ /** Destructor */
+ ~SARStreamingMaximumMinimumImageFilter() ITK_OVERRIDE {}
+
+private:
+ SARStreamingMaximumMinimumImageFilter(const Self &); //purposely not implemented
+ void operator =(const Self&); //purposely not implemented
+};
+
+} // end namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbSARStreamingMaximumMinimumImageFilter.txx"
+#endif
+
+#endif
diff --git a/include/otbSARStreamingMaximumMinimumImageFilter.txx b/include/otbSARStreamingMaximumMinimumImageFilter.txx
new file mode 100644
index 0000000..1a5840e
--- /dev/null
+++ b/include/otbSARStreamingMaximumMinimumImageFilter.txx
@@ -0,0 +1,415 @@
+/*
+ * Copyright (C) 2005-2018 Centre National d'Etudes Spatiales (CNES)
+ *
+ * This file is part of Orfeo Toolbox
+ *
+ * https://www.orfeo-toolbox.org/
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef otbSARStreamingMaximumMinimumImageFilter_txx
+#define otbSARStreamingMaximumMinimumImageFilter_txx
+#include "otbSARStreamingMaximumMinimumImageFilter.h"
+
+#include "itkImageScanlineConstIterator.h"
+#include "itkImageRegionIterator.h"
+#include "itkProgressReporter.h"
+#include "otbMacro.h"
+
+#include <complex>
+
+namespace otb
+{
+
+template<class TInputImage>
+PersistentMaximumMinimumImageFilter<TInputImage>
+::PersistentMaximumMinimumImageFilter()
+ :m_ThreadMin(1),
+ m_ThreadMax(1),
+ m_ThreadIndexMaxC(0),
+ m_ThreadIndexMaxL(0),
+ m_ThreadIndexMinC(0),
+ m_ThreadIndexMinL(0),
+ m_IgnoreInfiniteValues(true),
+ m_IgnoreUserDefinedValue(false)
+{
+ //this->itk::ProcessObject::SetNumberOfRequiredOutputs(3);
+
+ // first output is a copy of the image, DataObject created by
+ // superclass
+ //
+ // allocate the data objects for the outputs which are
+ // just decorators around unsigned int types (for index)
+ typename IndexObjectType::Pointer output1
+ = static_cast<IndexObjectType*>(this->MakeOutput(3).GetPointer());
+ this->itk::ProcessObject::SetNthOutput(1, output1.GetPointer());
+ typename IndexObjectType::Pointer output2
+ = static_cast<IndexObjectType*>(this->MakeOutput(3).GetPointer());
+ this->itk::ProcessObject::SetNthOutput(2, output2.GetPointer());
+
+ // For Min/Max
+ typename ValueObjectType::Pointer output3
+ = static_cast<ValueObjectType*>(this->MakeOutput(1).GetPointer());
+ this->itk::ProcessObject::SetNthOutput(3, output3.GetPointer());
+ typename ValueObjectType::Pointer output4
+ = static_cast<ValueObjectType*>(this->MakeOutput(1).GetPointer());
+ this->itk::ProcessObject::SetNthOutput(4, output4.GetPointer());
+
+ // Initialisation
+ IndexType indexInit;
+ indexInit[0] = 0;
+ indexInit[1] = 0;
+ this->GetIndexOfMinOutput()->Set(indexInit);
+ this->GetIndexOfMaxOutput()->Set(indexInit);
+ this->GetMinOutput()->Set(10E20);
+ this->GetMaxOutput()->Set(-10E20);
+
+ // Initiate the infinite ignored pixel counters
+ m_IgnoredInfinitePixelCount= std::vector<unsigned long>(this->GetNumberOfThreads(), 0);
+ m_IgnoredUserPixelCount= std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
+
+ this->Reset();
+}
+
+template<class TInputImage>
+typename itk::DataObject::Pointer
+PersistentMaximumMinimumImageFilter<TInputImage>
+::MakeOutput(DataObjectPointerArraySizeType output)
+{
+ switch (output)
+ {
+ case 0:
+ return static_cast<itk::DataObject*>(LongObjectType::New().GetPointer());
+ break;
+ case 1:
+ return static_cast<itk::DataObject*>(ValueObjectType::New().GetPointer());
+ break;
+ case 2:
+ return static_cast<itk::DataObject*>(PixelObjectType::New().GetPointer());
+ break;
+ case 3:
+ return static_cast<itk::DataObject*>(IndexObjectType::New().GetPointer());
+ break;
+ default:
+ // might as well make an image
+ return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
+ break;
+ }
+}
+
+
+template<class TInputImage>
+typename PersistentMaximumMinimumImageFilter<TInputImage>::IndexObjectType*
+PersistentMaximumMinimumImageFilter<TInputImage>
+::GetIndexOfMinOutput()
+{
+ return static_cast<IndexObjectType*>(this->itk::ProcessObject::GetOutput(1));
+}
+
+template<class TInputImage>
+const typename PersistentMaximumMinimumImageFilter<TInputImage>::IndexObjectType*
+PersistentMaximumMinimumImageFilter<TInputImage>
+::GetIndexOfMinOutput() const
+{
+ return static_cast<const IndexObjectType*>(this->itk::ProcessObject::GetOutput(1));
+}
+
+template<class TInputImage>
+typename PersistentMaximumMinimumImageFilter<TInputImage>::IndexObjectType*
+PersistentMaximumMinimumImageFilter<TInputImage>
+::GetIndexOfMaxOutput()
+{
+ return static_cast<IndexObjectType*>(this->itk::ProcessObject::GetOutput(2));
+}
+
+template<class TInputImage>
+const typename PersistentMaximumMinimumImageFilter<TInputImage>::IndexObjectType*
+PersistentMaximumMinimumImageFilter<TInputImage>
+::GetIndexOfMaxOutput() const
+{
+ return static_cast<const IndexObjectType*>(this->itk::ProcessObject::GetOutput(2));
+}
+
+template<class TInputImage>
+typename PersistentMaximumMinimumImageFilter<TInputImage>::ValueObjectType*
+PersistentMaximumMinimumImageFilter<TInputImage>
+::GetMinOutput()
+{
+ return static_cast<ValueObjectType*>(this->itk::ProcessObject::GetOutput(3));
+}
+
+template<class TInputImage>
+const typename PersistentMaximumMinimumImageFilter<TInputImage>::ValueObjectType*
+PersistentMaximumMinimumImageFilter<TInputImage>
+::GetMinOutput() const
+{
+ return static_cast<const ValueObjectType*>(this->itk::ProcessObject::GetOutput(3));
+}
+
+
+template<class TInputImage>
+typename PersistentMaximumMinimumImageFilter<TInputImage>::ValueObjectType*
+PersistentMaximumMinimumImageFilter<TInputImage>
+::GetMaxOutput()
+{
+ return static_cast<ValueObjectType*>(this->itk::ProcessObject::GetOutput(4));
+}
+
+template<class TInputImage>
+const typename PersistentMaximumMinimumImageFilter<TInputImage>::ValueObjectType*
+PersistentMaximumMinimumImageFilter<TInputImage>
+::GetMaxOutput() const
+{
+ return static_cast<const ValueObjectType*>(this->itk::ProcessObject::GetOutput(4));
+}
+
+
+template<class TInputImage>
+void
+PersistentMaximumMinimumImageFilter<TInputImage>
+::GenerateOutputInformation()
+{
+ Superclass::GenerateOutputInformation();
+ if (this->GetInput())
+ {
+ this->GetOutput()->CopyInformation(this->GetInput());
+ this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion());
+
+ if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
+ {
+ this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
+ }
+
+ }
+}
+template<class TInputImage>
+void
+PersistentMaximumMinimumImageFilter<TInputImage>
+::AllocateOutputs()
+{
+ // This is commented to prevent the streaming of the whole image for the first stream strip
+ // It shall not cause any problem because the output image of this filter is not intended to be used.
+ //InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
+ //this->GraftOutput( image );
+ // Nothing that needs to be allocated for the remaining outputs
+}
+
+template<class TInputImage>
+void
+PersistentMaximumMinimumImageFilter<TInputImage>
+::Synthetize()
+{
+ int numberOfThreads = this->GetNumberOfThreads();
+
+ // Init indexes and values
+ IndexType minIndex, maxIndex;
+ minIndex = this->GetIndexOfMin();
+ maxIndex = this->GetIndexOfMax();;
+
+ ValueType minValue, maxValue;
+ minValue = this->GetMin();
+ maxValue = this->GetMax();
+
+ // For each threads, get information into our Thread arrays
+ for (int i = 0; i < numberOfThreads; ++i)
+ {
+ // Conditions for min/max
+ if (minValue > m_ThreadMin[i])
+ {
+ minValue = m_ThreadMin[i];
+ minIndex[0] = m_ThreadIndexMinC[i];
+ minIndex[1] = m_ThreadIndexMinL[i];
+ }
+ if (maxValue < m_ThreadMax[i])
+ {
+ maxValue = m_ThreadMax[i];
+ maxIndex[0] = m_ThreadIndexMaxC[i];
+ maxIndex[1] = m_ThreadIndexMaxL[i];
+ }
+ }
+
+
+ // Set the outputs
+ this->GetIndexOfMinOutput()->Set(minIndex);
+ this->GetIndexOfMaxOutput()->Set(maxIndex);
+ this->GetMinOutput()->Set(minValue);
+ this->GetMaxOutput()->Set(maxValue);
+
+ std::cout << "maxValue : " << maxValue << std::endl;
+ std::cout << "this->GetMax() : " << this->GetMax() << std::endl;
+
+
+}
+
+template<class TInputImage>
+void
+PersistentMaximumMinimumImageFilter<TInputImage>
+::Reset()
+{
+ int numberOfThreads = this->GetNumberOfThreads();
+
+ // Resize the thread temporaries
+ m_ThreadIndexMaxC.SetSize(numberOfThreads);
+ m_ThreadIndexMaxL.SetSize(numberOfThreads);
+ m_ThreadIndexMinC.SetSize(numberOfThreads);
+ m_ThreadIndexMinL.SetSize(numberOfThreads);
+ m_ThreadMax.SetSize(numberOfThreads);
+ m_ThreadMin.SetSize(numberOfThreads);
+ // Initialize the temporaries
+ m_ThreadIndexMaxC.Fill(0);
+ m_ThreadIndexMaxL.Fill(0);
+ m_ThreadIndexMinC.Fill(0);
+ m_ThreadIndexMinL.Fill(0);
+ m_ThreadMin.Fill(10E20);
+ m_ThreadMax.Fill(-10E20);
+
+ if (m_IgnoreInfiniteValues)
+ {
+ m_IgnoredInfinitePixelCount= std::vector<unsigned long>(numberOfThreads, 0);
+ }
+
+ if (m_IgnoreUserDefinedValue)
+ {
+ m_IgnoredUserPixelCount= std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
+ }
+}
+
+
+
+/**
+ * Method OutputRegionToInputRegion for GenerateInputRequestedRegion
+ */
+template<class TInputImage >
+typename PersistentMaximumMinimumImageFilter< TInputImage >::RegionType
+PersistentMaximumMinimumImageFilter< TInputImage >
+::OutputRegionToInputRegion(const RegionType& outputRegion) const
+{
+ // Compute the input requested region (size and start index)
+ // Use the image transformations to insure an input requested region
+ // that will provide the proper range
+ const SizeType & outputRequestedRegionSize = outputRegion.GetSize();
+ const IndexType & outputRequestedRegionIndex = outputRegion.GetIndex();
+
+ // Compute input index
+ IndexType inputRequestedRegionIndex = outputRequestedRegionIndex;
+
+ // Compute input size
+ SizeType inputRequestedRegionSize = outputRequestedRegionSize;
+
+ // Input region
+ RegionType inputRequestedRegion = outputRegion;
+ inputRequestedRegion.SetIndex(inputRequestedRegionIndex);
+ inputRequestedRegion.SetSize(inputRequestedRegionSize);
+
+ return inputRequestedRegion;
+}
+/**
+ * Method OutputRegionToInputRegion
+ */
+template<class TInputImage>
+void
+PersistentMaximumMinimumImageFilter< TInputImage >
+::GenerateInputRequestedRegion()
+{
+ RegionType outputRequestedRegion = this->GetOutput()->GetRequestedRegion();
+ RegionType inputRequestedRegion = OutputRegionToInputRegion(outputRequestedRegion);
+
+ InputImagePointer inputPtr = const_cast< ImageType * >( this->GetInput() );
+
+ inputPtr->SetRequestedRegion(inputRequestedRegion);
+}
+
+
+
+template<class TInputImage>
+void
+PersistentMaximumMinimumImageFilter<TInputImage>
+::ThreadedGenerateData(const RegionType& outputRegionForThread,
+ itk::ThreadIdType threadId)
+{
+ // Grab the input
+ InputImagePointer inputPtr = const_cast<TInputImage *>(this->GetInput(0));
+ // Compute corresponding input region
+ RegionType inputRegionForThread = OutputRegionToInputRegion(outputRegionForThread);
+
+ // support progress methods/callbacks
+ itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
+
+ // Define/declare an iterator that will walk the input region for this
+ // thread.
+ typedef itk::ImageScanlineConstIterator< ImageType > InputIterator;
+ InputIterator inIt(this->GetInput(), inputRegionForThread);
+
+ inIt.GoToBegin();
+
+ ValueType valueCurrent;
+
+ // For each line
+ while ( !inIt.IsAtEnd())
+ {
+ inIt.GoToBeginOfLine();
+
+ // For each column
+ while (!inIt.IsAtEndOfLine())
+ {
+ //The input image value
+ valueCurrent= inIt.Get();
+
+ bool valueCurrentIsFinite = (vnl_math_isfinite(valueCurrent));
+
+ // Assigne if finite
+ if (valueCurrentIsFinite)
+ {
+ // Conditions for min/max
+ if (m_ThreadMin[threadId] > valueCurrent)
+ {
+ m_ThreadMin[threadId] = valueCurrent;
+ m_ThreadIndexMinC[threadId] = inIt.GetIndex()[0];
+ m_ThreadIndexMinL[threadId] = inIt.GetIndex()[1];
+ }
+ if (m_ThreadMax[threadId] < valueCurrent)
+ {
+ m_ThreadMax[threadId] = valueCurrent;
+ m_ThreadIndexMaxC[threadId] = inIt.GetIndex()[0];
+ m_ThreadIndexMaxL[threadId] = inIt.GetIndex()[1];
+ }
+ }
+
+
+ // Next colunm
+ ++inIt;
+ }
+ // Next Line
+ inIt.NextLine();
+
+ }
+
+}
+
+template <class TInputImage>
+void
+PersistentMaximumMinimumImageFilter<TInputImage>
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+ Superclass::PrintSelf(os, indent);
+
+ os << indent << "Max: "
+ << static_cast<typename itk::NumericTraits<ValueType>::PrintType>(this->GetMax()) << std::endl;
+ os << indent << "Min: "
+ << static_cast<typename itk::NumericTraits<ValueType>::PrintType>(this->GetMin()) << std::endl;
+ os << indent << "Index of Max: " << this->GetIndexOfMax() << std::endl;
+ os << indent << "Index of Min: " << this->GetIndexOfMin() << std::endl;
+}
+} // end namespace otb
+#endif
--
GitLab