Commit 05db7f05 authored by Cresson Remi's avatar Cresson Remi

WIP: SPS

parent 03bc4a2c
......@@ -12,4 +12,8 @@ OTB_CREATE_APPLICATION(NAME Harmonizer
SOURCES otbHarmonizer.cxx
LINK_LIBRARIES ${${otb-module}_LIBRARIES})
OTB_CREATE_APPLICATION(NAME SquarePatchesSelection
SOURCES otbSquarePatchesSelection.cxx
LINK_LIBRARIES ${${otb-module}_LIBRARIES})
/*=========================================================================
Copyright (c) Remi Cresson (IRSTEA). All rights reserved.
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 "otbWrapperApplicationFactory.h"
#include "otbSquarePatchesStatisticsFilter.h"
namespace otb
{
namespace Wrapper
{
class SquarePatchesSelection : public Application
{
public:
/** Standard class typedefs. */
typedef SquarePatchesSelection Self;
typedef Application Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
/** Standard macro */
itkNewMacro(Self);
itkTypeMacro(SquarePatchesSelection, Application);
typedef otb::SquarePatchesStatisticsFilter<FloatVectorImageType, UInt16VectorImageType> FilterType;
private:
void DoInit() override
{
SetName("SquarePatchesSelection");
SetDescription("Compute stats on patches");
// Documentation
SetDocLongDescription("This application computes some stats over patches.");
SetDocLimitations("None");
SetDocAuthors("Remi Cresson");
SetDocSeeAlso(" ");
AddParameter(ParameterType_InputImage, "in", "input");
AddParameter(ParameterType_OutputImage, "out", "Output image");
AddParameter(ParameterType_Int, "patchsize", "Patch size");
SetDefaultParameterInt ("patchsize", 64);
SetMinimumParameterIntValue ("patchsize", 2);
AddParameter(ParameterType_Float, "nodata", "Nodata value");
SetDefaultParameterFloat("nodata", 0.);
AddRAMParameter();
}
void DoUpdateParameters() override
{
}
void DoExecute() override
{
// Get the input
FloatVectorImageType* img = GetParameterImage("in");
// Filter
m_Filter = FilterType::New();
m_Filter->SetInput(img);
m_Filter->SetNoDataValue(GetParameterFloat("nodata"));
m_Filter->SetPatchSize(GetParameterInt("patchsize"));
// Set the output image
SetParameterOutputImage("out", m_Filter->GetOutput());
}
FilterType::Pointer m_Filter;
};
} // end namespace Wrapper
} // end namespace otb
OTB_APPLICATION_EXPORT(otb::Wrapper::SquarePatchesSelection)
/*=========================================================================
Copyright (c) Remi Cresson (IRSTEA). All rights reserved.
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 __SquarePatchesStatisticsFilter_H
#define __SquarePatchesStatisticsFilter_H
#include "itkImageToImageFilter.h"
#include <vector>
#include "vnl/vnl_matrix.h"
namespace otb
{
/** \class SquarePatchesStatisticsFilter
* \brief Computes some stats in the patches.
* Patches are sampled as a matrix of adjacent patches.
*
* \ingroup MLUtils
*/
template <class TInputImage, class TOutputImage>
class ITK_EXPORT SquarePatchesStatisticsFilter : public itk::ImageToImageFilter<TInputImage, TOutputImage>
{
public:
/** Standard Self typedef */
typedef SquarePatchesStatisticsFilter Self;
typedef itk::ImageToImageFilter<TInputImage, TOutputImage> 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(SquarePatchesStatisticsFilter, itk::ImageToImageFilter);
/** Input image typedefs. */
typedef TInputImage InputImageType;
typedef typename InputImageType::PixelType InputImagePixelType;
typedef typename InputImageType::InternalPixelType InputImageValueType;
typedef typename itk::ImageRegionConstIterator<TInputImage> IteratorType;
typedef typename InputImageType::RegionType InputImageRegionType;
typedef typename InputImageType::IndexType InputImageIndexType;
typedef typename InputImageType::SizeType InputImageSizeType;
typedef typename InputImageType::PointType InputImagePointType;
typedef typename InputImageType::Pointer InputImagePointer;
typedef typename InputImageType::SpacingType InputImageSpacingType;
/** Output image typedefs. */
typedef TOutputImage OutputImageType;
typedef typename OutputImageType::PointType OutputImagePointType;
typedef typename OutputImageType::PixelType OutputImagePixelType;
typedef typename OutputImageType::InternalPixelType OutputImageInternalPixelType;
typedef typename OutputImageType::RegionType OutputImageRegionType;
typedef typename itk::ImageRegionIterator<TOutputImage> OutIteratorType;
typedef vnl_matrix<OutputImageInternalPixelType> CountMatrixType;
itkSetMacro(PatchSize, unsigned int);
itkSetMacro(NoDataValue, InputImageValueType);
protected:
SquarePatchesStatisticsFilter() {
m_NoDataValue = 0;
m_PatchSize = 0;
}
virtual ~SquarePatchesStatisticsFilter() {
}
/** Overrided methods */
virtual void GenerateOutputInformation();
virtual void GenerateInputRequestedRegion();
virtual void BeforeThreadedGenerateData();
virtual void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId );
virtual void AfterThreadedGenerateData();
private:
SquarePatchesStatisticsFilter(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
unsigned int m_PatchSize;
InputImageValueType m_NoDataValue;
std::vector<CountMatrixType> m_ThreadsMasks;
}; // end of class
} // end namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
#include "otbSquarePatchesStatisticsFilter.hxx"
#endif
#endif
/*=========================================================================
Copyright (c) Remi Cresson (IRSTEA). All rights reserved.
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 __SquarePatchesStatisticsFilter_txx
#define __SquarePatchesStatisticsFilter_txx
#include "otbSquarePatchesStatisticsFilter.h"
namespace otb {
template <class TInputImage, class TOutputImage>
void
SquarePatchesStatisticsFilter<TInputImage, TOutputImage>
::GenerateOutputInformation()
{
Superclass::GenerateOutputInformation();
//////////////////////////////////////////////////////////////////////////////////////////
// Compute the output image extent
//////////////////////////////////////////////////////////////////////////////////////////
// Spacing
const InputImageSpacingType m_InputSpacing = this->GetInput()->GetSignedSpacing();
InputImageSpacingType m_OutputSpacing = m_InputSpacing;
for (unsigned int dim = 0 ; dim < 2 ; dim++)
m_OutputSpacing[dim] *= m_PatchSize;
// Origin
InputImagePointType m_OutputOrigin = this->GetInput()->GetOrigin();
for (unsigned int dim = 0 ; dim < 2 ; dim++)
{
m_OutputOrigin[dim] -= 0.5 * m_InputSpacing[dim];
m_OutputOrigin[dim] += 0.5 * m_OutputSpacing[dim];
}
// Size
InputImageSizeType m_OutputSize = this->GetInput()->GetLargestPossibleRegion().GetSize();
for (unsigned int dim = 0 ; dim < 2 ; dim++)
m_OutputSize[dim] /= m_PatchSize;
InputImageRegionType largestPossibleRegion;
largestPossibleRegion.SetSize(m_OutputSize);
// Number of components per pixel
unsigned int outputPixelSize = 1;
// Copy input image projection
InputImageType * inputImage = static_cast<InputImageType * >( Superclass::ProcessObject::GetInput(0) );
const std::string projectionRef = inputImage->GetProjectionRef();
// Set output image origin/spacing/size/projection
OutputImageType * outputPtr = this->GetOutput();
outputPtr->SetNumberOfComponentsPerPixel(outputPixelSize);
outputPtr->SetProjectionRef ( projectionRef );
outputPtr->SetOrigin ( m_OutputOrigin );
outputPtr->SetSignedSpacing ( m_OutputSpacing );
outputPtr->SetLargestPossibleRegion( largestPossibleRegion );
}
template <class TInputImage, class TOutputImage>
void
SquarePatchesStatisticsFilter<TInputImage, TOutputImage>
::GenerateInputRequestedRegion()
{
Superclass::GenerateInputRequestedRegion();
// Output requested region
OutputImageRegionType requestedRegion = this->GetOutput()->GetRequestedRegion();
// Input requested region
for (unsigned int dim = 0 ; dim < 2 ; dim++)
{
requestedRegion.GetModifiableIndex()[dim] *= m_PatchSize;
requestedRegion.GetModifiableSize()[dim] *= m_PatchSize;
}
InputImageType * inputImage = static_cast<InputImageType * >(Superclass::ProcessObject::GetInput(0) );
inputImage->SetRequestedRegion(requestedRegion);
}
/*
* Before threaded generation
*/
template <class TInputImage, class TOutputImage>
void
SquarePatchesStatisticsFilter<TInputImage, TOutputImage>
::BeforeThreadedGenerateData()
{
m_ThreadsMasks.clear();
const InputImageRegionType reqRegion = this->GetInput()->GetRequestedRegion();
const InputImageSizeType reqSize = reqRegion.GetSize();
for (unsigned int i = 0; i < this->GetNumberOfThreads(); i++)
m_ThreadsMasks.push_back(CountMatrixType(reqSize[0], reqSize[1], 0));
}
/**
* Processing
*/
template <class TInputImage, class TOutputImage>
void
SquarePatchesStatisticsFilter<TInputImage, TOutputImage>
::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
{
// Output region --> Input region
InputImageIndexType newStart = outputRegionForThread.GetIndex();
InputImageSizeType newSize = outputRegionForThread.GetSize();
for (unsigned int dim = 0; dim < 2; dim++)
{
newStart[dim] *= m_PatchSize;
newSize[dim] *= m_PatchSize;
}
InputImageRegionType inputRegion(newStart, newSize);
// Iterate through the input image
const InputImagePointer inputPtr = const_cast<TInputImage*>(this->GetInput());
IteratorType it(inputPtr, inputRegion);
for (it.GoToBegin(); !it.IsAtEnd(); ++it)
{
const InputImageValueType val = it.Get()[0];
const InputImageIndexType idx = it.GetIndex();
m_ThreadsMasks[threadId][idx[0]/m_PatchSize][idx[1]/m_PatchSize] += static_cast<OutputImageInternalPixelType>(val != m_NoDataValue);
}
}
template <class TInputImage, class TOutputImage>
void
SquarePatchesStatisticsFilter<TInputImage, TOutputImage>
::AfterThreadedGenerateData()
{
const InputImagePointer inputPtr = const_cast<TInputImage*>(this->GetInput());
const InputImageRegionType reqRegion = inputPtr->GetRequestedRegion();
const InputImageSizeType reqSize = reqRegion.GetSize();
CountMatrixType result(reqSize[0], reqSize[1], 0);
// Reduce threads
unsigned int i,j;
for (auto& threadMask: m_ThreadsMasks)
for (i=0; i<reqSize[0]; i++)
for (j=0; j<reqSize[1]; j++)
result[i][j] += threadMask[i][j];
// Set output
typename TOutputImage::Pointer outputPtr = this->GetOutput();
const OutputImageRegionType outputReqRegion = outputPtr->GetRequestedRegion();
OutIteratorType it(outputPtr, outputReqRegion);
OutputImagePixelType outPix;
outPix.SetSize(1);
for (it.GoToBegin(); !it.IsAtEnd(); ++it)
{
const InputImageIndexType idx = it.GetIndex();
outPix[0] = result[(unsigned int) idx[0]][(unsigned int) idx[1]];
it.Set(outPix);
}
}
} // end namespace gtb
#endif
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment