Commit 5b16b4a2 authored by Cresson Remi's avatar Cresson Remi
Browse files

ADD: new application

parent 680b62bc
......@@ -8,3 +8,8 @@ OTB_CREATE_APPLICATION(NAME SuperSuperimpose
SOURCES otbSuperSuperimpose.cxx
LINK_LIBRARIES OTBCommon)
OTB_CREATE_APPLICATION(NAME Harmonizer
SOURCES otbHarmonizer.cxx
LINK_LIBRARIES OTBCommon)
/*=========================================================================
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 "otbRegionComparator.h"
#include "itkUnaryFunctorImageFilter.h"
#include "otbExtractROI.h"
#include "otbMultiChannelExtractROI.h"
#include "otbShiftScaleVectorImageFilter.h"
// Masking
#include "itkBinaryErodeImageFilter.h"
#include "itkFlatStructuringElement.h"
#include "itkMaskImageFilter.h"
// Smoothing
#include "otbGridResampleImageFilter.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
#include "otbBCOInterpolateImageFunction.h"
// Stats
#include "otbStreamingStatisticsMosaicFilter.h"
namespace otb
{
namespace Wrapper
{
class Harmonizer : public Application
{
public:
// Functor for linear stretching
template<class TPixel, class OutputPixel>
class LinearStretch
{
public:
LinearStretch(){
m_Gain = 1.0;
m_Offset = 0.0;
}
~LinearStretch(){}
inline OutputPixel operator()( const TPixel & in ) const
{
OutputPixel out(in);
for (unsigned int band = 0 ; band < in.Size() ; band++)
{
if (in[band] == m_NoDataValue)
out[band] = m_NoDataValue;
else
{
out[band] *= m_Gain [band];
out[band] += m_Offset [band];
}
}
return out;
}
void SetNoDataValue(typename TPixel::ValueType value)
{
m_NoDataValue = value;
}
void SetGain(TPixel pix)
{
m_Gain = pix;
}
void SetOffset(TPixel pix)
{
m_Offset = pix;
}
private:
typename TPixel::ValueType m_NoDataValue;
TPixel m_Gain;
TPixel m_Offset;
};
/** Standard class typedefs. */
typedef Harmonizer Self;
typedef Application Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
/** Standard macro */
itkNewMacro(Self);
itkTypeMacro(Harmonizer, Application);
// ROI
typedef otb::MultiChannelExtractROI<FloatVectorImageType::InternalPixelType,
FloatVectorImageType::InternalPixelType> ExtractROIFilterType;
typedef itk::NearestNeighborInterpolateImageFunction<FloatVectorImageType> NNInterpolatorType;
typedef otb::GridResampleImageFilter<FloatVectorImageType, FloatVectorImageType> ResampleFilterType;
typedef itk::NearestNeighborInterpolateImageFunction<UInt8ImageType> MaskNNInterpolatorType;
typedef otb::GridResampleImageFilter<UInt8ImageType, UInt8ImageType> MaskResampleFilterType;
typedef itk::MaskImageFilter<FloatVectorImageType, UInt8ImageType, FloatVectorImageType> MaskImageFilterType;
typedef otb::ExtractROI<UInt8ImageType::PixelType, UInt8ImageType::PixelType> MaskExtractROIFilterType;
// Statistics
typedef otb::StreamingStatisticsMosaicFilter<FloatVectorImageType, FloatVectorImageType, double> StatsFilterType;
// Shiftscale
typedef LinearStretch<FloatVectorImageType::PixelType, FloatVectorImageType::PixelType> LinStretchFunctorType;
typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatVectorImageType,
LinStretchFunctorType> LinStretchFilterType;
private:
// Functor to retrieve nodata
template<class TPixel, class OutputPixel>
class IsNoData
{
public:
IsNoData(){}
~IsNoData(){}
inline OutputPixel operator()( const TPixel & A ) const
{
for (unsigned int band = 0 ; band < A.Size() ; band++)
{
if (A[band] != m_NoDataValue)
return 1;
}
return 0;
}
void SetNoDataValue(typename TPixel::ValueType value)
{
m_NoDataValue = value;
}
private:
typename TPixel::ValueType m_NoDataValue;
};
class SmoothValidPipeline
{
public:
// Smoothing
typedef otb::BCOInterpolateImageFunction<FloatVectorImageType> BCOInterpolatorType;
typedef otb::GridResampleImageFilter<FloatVectorImageType, FloatVectorImageType> BCOResampleFilterType;
// Masking
typedef IsNoData<FloatVectorImageType::PixelType, UInt8ImageType::PixelType> IsNoDataFunctorType;
typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, UInt8ImageType,
IsNoDataFunctorType> IsNoDataFilterType;
typedef itk::FlatStructuringElement<2> StructuringType;
typedef StructuringType::RadiusType RadiusType;
typedef itk::BinaryErodeImageFilter<UInt8ImageType, UInt8ImageType, StructuringType> ErodeFilterType;
// Mini-pipeline that decimates efficiently an image and which takes cares of no-data
SmoothValidPipeline(FloatVectorImageType::Pointer & input, float nodata, int radius){
// Image --> mask
m_NoDataFilter = IsNoDataFilterType::New();
m_NoDataFilter->GetFunctor().SetNoDataValue(nodata);
m_NoDataFilter->SetInput(input);
m_NoDataFilter->UpdateOutputInformation();
// mask --> padded mask
m_MaskResampleFilter1 = MaskResampleFilterType::New();
m_MaskResampleFilter1->SetInput(m_NoDataFilter->GetOutput() );
UInt8ImageType::SizeType outputSize = m_NoDataFilter->GetOutput()->GetLargestPossibleRegion().GetSize();
UInt8ImageType::SpacingType outputSpacing = m_NoDataFilter->GetOutput()->GetSignedSpacing();
UInt8ImageType::PointType outputOrigin = m_NoDataFilter->GetOutput()->GetOrigin();
for (int dim = 0 ; dim < 2 ; dim++)
{
outputSize[dim] += radius * 2 ;
outputOrigin[dim] -= radius * outputSpacing[dim];
}
m_MaskResampleFilter1->SetOutputSize(outputSize);
m_MaskResampleFilter1->SetOutputStartIndex(m_NoDataFilter->GetOutput()->GetLargestPossibleRegion().GetIndex());
m_MaskResampleFilter1->SetOutputSpacing(outputSpacing);
m_MaskResampleFilter1->SetOutputOrigin(outputOrigin);
MaskNNInterpolatorType::Pointer nnInterpolator = MaskNNInterpolatorType::New();
m_MaskResampleFilter1->SetInterpolator(nnInterpolator);
m_MaskResampleFilter1->UpdateOutputInformation();
// padded mask --> eroded mask
RadiusType rad;
rad.Fill(radius);
StructuringType se = StructuringType::Ball(rad);
m_ErodeFilter = ErodeFilterType::New();
m_ErodeFilter->SetInput(m_MaskResampleFilter1->GetOutput());
m_ErodeFilter->SetKernel(se);
m_ErodeFilter->SetForegroundValue(1);
m_ErodeFilter->SetBackgroundValue(0);
// Smoothing
m_SmoothingFilter = BCOResampleFilterType::New();
m_SmoothingFilter->SetInput(input);
BCOInterpolatorType::Pointer interpolator = BCOInterpolatorType::New();
interpolator->SetRadius(radius);
m_SmoothingFilter->SetInterpolator(interpolator);
FloatVectorImageType::SpacingType spacing = input->GetSignedSpacing();
FloatVectorImageType::IndexType start = input->GetLargestPossibleRegion().GetIndex();
FloatVectorImageType::SizeType size = input->GetLargestPossibleRegion().GetSize();
FloatVectorImageType::PointType origin = input->GetOrigin();
for (int dim = 0 ; dim < 2 ; dim++)
{
spacing[dim] *= radius;
size[dim] /= radius;
origin[dim] += 0.5*(spacing[dim] - input->GetSignedSpacing()[dim]);
}
m_SmoothingFilter->SetOutputOrigin(origin);
m_SmoothingFilter->SetOutputSpacing(spacing);
m_SmoothingFilter->SetOutputSize(size);
m_SmoothingFilter->SetOutputStartIndex(start);
m_SmoothingFilter->UpdateOutputInformation();
// Masking the smoothed image with the eroded mask
m_MaskResampleFilter2 = MaskResampleFilterType::New();
m_MaskResampleFilter2->SetInput(m_ErodeFilter->GetOutput() );
m_MaskResampleFilter2->SetInterpolator(nnInterpolator);
m_MaskResampleFilter2->SetOutputOrigin(m_SmoothingFilter->GetOutput()->GetOrigin());
m_MaskResampleFilter2->SetOutputSpacing(m_SmoothingFilter->GetOutput()->GetSignedSpacing());
m_MaskResampleFilter2->SetOutputSize(m_SmoothingFilter->GetOutput()->GetLargestPossibleRegion().GetSize());
m_MaskResampleFilter2->SetOutputStartIndex(m_SmoothingFilter->GetOutput()->GetLargestPossibleRegion().GetIndex());
m_MaskResampleFilter2->UpdateOutputInformation();
m_MaskFilter = MaskImageFilterType::New();
m_MaskFilter->SetInput(m_SmoothingFilter->GetOutput());
m_MaskFilter->SetMaskImage(m_MaskResampleFilter2->GetOutput());
m_MaskFilter->UpdateOutputInformation();
};
virtual ~SmoothValidPipeline (){};
FloatVectorImageType::Pointer GetOutput()
{
return m_MaskFilter->GetOutput();
}
UInt8ImageType::Pointer GetOutputMask()
{
return m_MaskResampleFilter2->GetOutput();
}
private:
IsNoDataFilterType::Pointer m_NoDataFilter;
MaskResampleFilterType::Pointer m_MaskResampleFilter1;
ErodeFilterType::Pointer m_ErodeFilter;
BCOResampleFilterType::Pointer m_SmoothingFilter;
MaskResampleFilterType::Pointer m_MaskResampleFilter2;
MaskImageFilterType::Pointer m_MaskFilter;
};
void DoInit() override
{
SetName("Harmonizer");
SetDescription("Harmonize a target image with a reference image");
// Documentation
SetDocName("Harmonizer");
SetDocLongDescription("This application performs the harmonization of an image radiometry from a reference.");
SetDocLimitations("None");
SetDocAuthors("Remi Cresson");
SetDocSeeAlso(" ");
AddDocTag(Tags::Geometry);
AddDocTag("Harmonization");
AddParameter(ParameterType_InputImage, "inr", "Reference input");
SetParameterDescription ("inr","The input reference image.");
AddParameter(ParameterType_InputImage, "inm", "The image to correct");
SetParameterDescription ("inm","The image to correct into the radiometry of the reference input.");
AddParameter(ParameterType_OutputImage, "out", "Output image");
SetParameterDescription ("out", "Output corrected image.");
MandatoryOff ("out");
AddParameter(ParameterType_Float, "inrnodata", "No-data value (reference)");
SetDefaultParameterFloat ("inrnodata", 0);
MandatoryOff ("inrnodata");
AddParameter(ParameterType_Float, "inmnodata", "No-data value (image to correct)");
SetDefaultParameterFloat ("inmnodata", 0);
MandatoryOff ("inmnodata");
AddParameter(ParameterType_Int, "smoothingradius", "Smoothing radius (in pixels of the reference image)");
SetDefaultParameterInt ("smoothingradius", 10);
AddRAMParameter();
// Doc example parameter settings
SetDocExampleParameterValue("inr", "QB_Toulouse_Ortho_XS1.tif");
SetDocExampleParameterValue("inm", "QB_Toulouse_Ortho_XS2.tif");
SetDocExampleParameterValue("out", "Harmonizerd_XS2_to_XS1.tif");
SetOfficialDocLink();
}
void DoUpdateParameters() override
{
}
void DoExecute() override
{
// Get the inputs
FloatVectorImageType::Pointer refImage = GetParameterImage("inr");
FloatVectorImageType::Pointer movImage = GetParameterImage("inm");
// Compute the radius for the moving image
m_Pipelines.clear();
int radiusForRef = GetParameterInt("smoothingradius");
float refSpacing = 0;
float movSpacing = 0;
for (unsigned int dim = 0 ; dim < 2 ; dim++)
{
float refV = refImage->GetSignedSpacing()[dim];
float movV = movImage->GetSignedSpacing()[dim];
refSpacing += refV*refV;
movSpacing += movV*movV;
}
refSpacing = std::sqrt(refSpacing);
movSpacing = std::sqrt(movSpacing);
unsigned int radiusForMov = ((float) radiusForRef) / movSpacing * refSpacing;
otbAppLogINFO("Estimated smoothing radius for the image to correct: " << radiusForMov);
SmoothValidPipeline pipelineRef (refImage, GetParameterFloat("inrnodata"), radiusForRef);
SmoothValidPipeline pipelineMov (movImage, GetParameterFloat("inmnodata"), radiusForMov);
m_Pipelines.push_back(pipelineRef);
m_Pipelines.push_back(pipelineMov);
// Compute rasters intersection region, check overlap
otb::RegionComparator<FloatVectorImageType, FloatVectorImageType> comparator;
comparator.SetImage1(pipelineRef.GetOutput());
comparator.SetImage2(pipelineMov.GetOutput());
if (!comparator.DoesOverlap())
{
otbAppLogFATAL("Inputs do not overlap!");
}
// Initialize ROI extract filters
m_ExtractROIFilterForRef = ExtractROIFilterType::New();
m_ExtractROIFilterForRef->SetInput(pipelineRef.GetOutput());
m_ExtractROIFilterForRef->SetExtractionRegion(comparator.GetOverlapInImage1Indices());
m_ExtractROIFilterForRef->UpdateOutputInformation();
m_MaskExtractROIForRef = MaskExtractROIFilterType::New();
m_MaskExtractROIForRef->SetInput(pipelineRef.GetOutputMask());
m_MaskExtractROIForRef->SetExtractionRegion(comparator.GetOverlapInImage1Indices());
m_MaskExtractROIForRef->UpdateOutputInformation();
m_ResampleFilterForMov = ResampleFilterType::New();
m_ResampleFilterForMov->SetInput(pipelineMov.GetOutput());
m_ResampleFilterForMov->SetOutputOrigin (m_ExtractROIFilterForRef->GetOutput()->GetOrigin());
m_ResampleFilterForMov->SetOutputSpacing(m_ExtractROIFilterForRef->GetOutput()->GetSignedSpacing());
m_ResampleFilterForMov->SetOutputStartIndex(m_ExtractROIFilterForRef->GetOutput()->GetLargestPossibleRegion().GetIndex());
m_ResampleFilterForMov->SetOutputSize(m_ExtractROIFilterForRef->GetOutput()->GetLargestPossibleRegion().GetSize());
NNInterpolatorType::Pointer nnInterpolator = NNInterpolatorType::New();
m_ResampleFilterForMov->SetInterpolator(nnInterpolator);
m_ResampleFilterForMov->UpdateOutputInformation();
m_ResampleFilterForMovMsk = MaskResampleFilterType::New();
m_ResampleFilterForMovMsk->SetInput(pipelineMov.GetOutputMask());
m_ResampleFilterForMovMsk->SetOutputOrigin(m_ExtractROIFilterForRef->GetOutput()->GetOrigin());
m_ResampleFilterForMovMsk->SetOutputSpacing(m_ExtractROIFilterForRef->GetOutput()->GetSignedSpacing());
m_ResampleFilterForMovMsk->SetOutputStartIndex(m_ExtractROIFilterForRef->GetOutput()->GetLargestPossibleRegion().GetIndex());
m_ResampleFilterForMovMsk->SetOutputSize(m_ExtractROIFilterForRef->GetOutput()->GetLargestPossibleRegion().GetSize());
MaskNNInterpolatorType::Pointer masknnInterpolator = MaskNNInterpolatorType::New();
m_ResampleFilterForMovMsk->SetInterpolator(masknnInterpolator);
m_ResampleFilterForMovMsk->UpdateOutputInformation();
m_MaskFilterForMov = MaskImageFilterType::New();
m_MaskFilterForMov->SetInput(m_ResampleFilterForMov->GetOutput());
m_MaskFilterForMov->SetMaskImage(m_MaskExtractROIForRef->GetOutput());
m_MaskFilterForMov->UpdateOutputInformation();
m_MaskFilterForRef = MaskImageFilterType::New();
m_MaskFilterForRef->SetInput(m_ExtractROIFilterForRef->GetOutput());
m_MaskFilterForRef->SetMaskImage(m_ResampleFilterForMovMsk->GetOutput());
m_MaskFilterForRef->UpdateOutputInformation();
m_StatsFilter = StatsFilterType::New();
m_StatsFilter->PushBackInput(m_MaskFilterForRef->GetOutput());
m_StatsFilter->PushBackInput(m_MaskFilterForMov->GetOutput());
// SetParameterOutputImage("out", m_MaskFilterForMov->GetOutput());
m_StatsFilter->GetStreamer()->SetAutomaticAdaptativeStreaming(GetParameterInt("ram"));
AddProcess(m_StatsFilter->GetStreamer(), "Computing statistics");
m_StatsFilter->Update();
otbAppLogINFO("Statistics");
FloatVectorImageType::PixelType gain, offset;
gain.SetSize(refImage->GetNumberOfComponentsPerPixel());
offset.SetSize(refImage->GetNumberOfComponentsPerPixel());
for (unsigned int band = 0 ; band < refImage->GetNumberOfComponentsPerPixel() ; band++)
{
float meanY = m_StatsFilter->GetMeans().at(band)[0][0];
float meanX = m_StatsFilter->GetMeans().at(band)[1][1];
float meanXY = m_StatsFilter->GetMeansOfProducts().at(band)[0][1];
float meanX2 = m_StatsFilter->GetMeansOfProducts().at(band)[1][1];
float mean2X = meanX * meanX;
float b1 = (meanX * meanY - meanXY) / (mean2X - meanX2);
float b0 = meanY - b1 * meanX;
otbAppLogINFO("Band " << band << " gain: " << b1 << " bias: " << b0);
gain[band] = b1;
offset[band] = b0;
}
if (this->HasValue("out"))
{
m_LinStretchFilter = LinStretchFilterType::New();
m_LinStretchFilter->SetInput(movImage);
m_LinStretchFilter->GetFunctor().SetNoDataValue(GetParameterFloat("inmnodata"));
m_LinStretchFilter->GetFunctor().SetGain(gain);
m_LinStretchFilter->GetFunctor().SetOffset(offset);
SetParameterOutputImage("out", m_LinStretchFilter->GetOutput());
}
}
std::vector<SmoothValidPipeline> m_Pipelines;
ExtractROIFilterType::Pointer m_ExtractROIFilterForRef;
MaskResampleFilterType::Pointer m_ResampleFilterForMovMsk;
ResampleFilterType::Pointer m_ResampleFilterForMov;
StatsFilterType::Pointer m_StatsFilter;
LinStretchFilterType::Pointer m_LinStretchFilter;
MaskImageFilterType::Pointer m_MaskFilterForMov;
MaskImageFilterType::Pointer m_MaskFilterForRef;
MaskExtractROIFilterType::Pointer m_MaskExtractROIForRef;
};
} // end namespace Wrapper
} // end namespace otb
OTB_APPLICATION_EXPORT(otb::Wrapper::Harmonizer)
......@@ -6,6 +6,7 @@ otb_module(MLUtils
OTBCommon
OTBApplicationEngine
Mosaic
SimpleExtractionTools
TEST_DEPENDS
OTBTestKernel
......
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