Skip to content
GitLab
Select Git revision
  • master default
Find file BlameHistoryPermalink
otbHarmonizer.cxx 18.00 KiB
1 
/*=========================================================================
2
3 
  Copyright (c) Remi Cresson (IRSTEA). All rights reserved.
4
5
6 
     This software is distributed WITHOUT ANY WARRANTY; without even
7 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
8 
     PURPOSE.  See the above copyright notices for more information.
9
10 
=========================================================================*/
11 
#include "otbWrapperApplicationFactory.h"
12
13 
#include "otbRegionComparator.h"
14 
#include "itkUnaryFunctorImageFilter.h"
15 
#include "otbExtractROI.h"
16 
#include "otbMultiChannelExtractROI.h"
17 
#include "otbShiftScaleVectorImageFilter.h"
18
19 
// Masking
20 
#include "itkBinaryErodeImageFilter.h"
21 
#include "itkFlatStructuringElement.h"
22 
#include "itkMaskImageFilter.h"
23
24 
// Smoothing
25 
#include "otbGridResampleImageFilter.h"
26 
#include "itkNearestNeighborInterpolateImageFunction.h"
27 
#include "otbBCOInterpolateImageFunction.h"
28
29 
// Stats
30 
#include "otbStreamingStatisticsMosaicFilter.h"
31
32 
namespace otb
33 
{
34
35
36 
namespace Wrapper
37 
{
38
39 
class Harmonizer : public Application
40 
{
41 
public:
42
43 
  // Functor for linear stretching
44 
  template<class TPixel, class OutputPixel>
45 
  class LinearStretch
46 
  {
47 
  public:
48 
    LinearStretch(){
49 
      m_Gain = 1.0;
50 
      m_Offset = 0.0;
51 
    }
52 
    ~LinearStretch(){}
53
54 
    inline OutputPixel operator()( const TPixel & in ) const
55 
    {
56 
      OutputPixel out(in);
57 
      for (unsigned int band = 0 ; band < in.Size() ; band++)
58 
      {
59 
        if (in[band] == m_NoDataValue)
60 
          out[band] = m_NoDataValue;
61 
        else
62 
        {
63 
          out[band] *= m_Gain   [band];
64 
          out[band] += m_Offset [band];
65 
        }
66 
      }
67 
      return out;
68 
    }
69
70 
    void SetNoDataValue(typename TPixel::ValueType value)
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
{ 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)
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
{ 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();
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
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 SetDocLongDescription("This application performs the harmonization of an image radiometry from a reference."); SetDocLimitations("None"); SetDocAuthors("Remi Cresson"); SetDocSeeAlso(" "); AddDocTag(Tags::Calibration); AddDocTag("Harmonization"); AddParameter(ParameterType_InputImage, "inr", "Reference input"); SetParameterDescription ("inr","The input reference image.");
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350
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); AddParameter(ParameterType_Choice, "zeroy", "Zero-y intercept mode in linear transform"); AddChoice("zeroy.on", "Yes (y=ax)"); AddChoice("zeroy.off", "No (y=ax+b)"); 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());
351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420
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()); std::stringstream exp; exp << "{"; bool zeroy = (GetParameterInt("zeroy") == 0); if (zeroy) { otbAppLogINFO("Using zero-y intercept (y=ax)"); } else {
421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489
otbAppLogINFO("Using linear correction model (y=ax+b)"); } 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 b0, b1; if (zeroy) { b1 = meanXY / meanX2; b0 = 0; } else { b1 = (meanX * meanY - meanXY) / (mean2X - meanX2); b0 = meanY - b1 * meanX; } otbAppLogINFO("Band " << band << " gain: " << b1 << " bias: " << b0); if (band>0) exp << ";"; exp << "im1b" << (band+1) << "*" << b1; if (b0>0) exp << "+" << b0; else exp << "" << b0; gain[band] = b1; offset[band] = b0; } exp << "}"; otbAppLogINFO("BandMathX expression: " << exp.str()); 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)