Commit 2549defd authored by Julien Michel's avatar Julien Michel
Browse files

Merge branch 'grid-resampling' into develop

parents 1f8e5995 8668075c
......@@ -19,6 +19,7 @@
#include "otbMultiToMonoChannelExtractROI.h"
#include "otbGenericRSResampleImageFilter.h"
#include "otbGridResampleImageFilter.h"
#include "otbImportGeoInformationImageFilter.h"
#include "otbBCOInterpolateImageFunction.h"
#include "otbSimpleRcsPanSharpeningFusionImageFilter.h"
......@@ -133,7 +134,7 @@ private:
typedef otb::BCOInterpolateImageFunction<FloatVectorImageType> InterpolatorType;
typedef otb::GenericRSResampleImageFilter<FloatVectorImageType, FloatVectorImageType> ResamplerType;
typedef otb::StreamingResampleImageFilter<FloatVectorImageType, FloatVectorImageType> BasicResamplerType;
typedef otb::GridResampleImageFilter<FloatVectorImageType, FloatVectorImageType> BasicResamplerType;
typedef otb::ImportGeoInformationImageFilter<FloatVectorImageType,InternalImageType> ImportGeoInformationFilterType;
typedef otb::SimpleRcsPanSharpeningFusionImageFilter<InternalImageType, FloatVectorImageType, FloatVectorImageType> FusionFilterType;
......@@ -204,13 +205,22 @@ private:
{
otbAppLogINFO("Using the PHR mode");
otb::PleiadesPToXSAffineTransformCalculator::TransformType::Pointer transform
= otb::PleiadesPToXSAffineTransformCalculator::Compute(panchro, xs);
otb::PleiadesPToXSAffineTransformCalculator::TransformType::OffsetType offset
= otb::PleiadesPToXSAffineTransformCalculator::ComputeOffset(GetParameterImage("inp"),
GetParameterImage("inxs"));
origin+=offset;
origin[0]=origin[0]/4;
origin[1]=origin[1]/4;
basicResampler->SetOutputOrigin(origin);
basicResampler->SetInput(xs);
basicResampler->SetTransform(transform);
basicResampler->SetOutputOrigin(origin);
basicResampler->SetOutputSpacing(spacing);
FloatVectorImageType::SpacingType xsSpacing = GetParameterImage("inxs")->GetSpacing();
xsSpacing*=0.25;
basicResampler->SetOutputSpacing(xsSpacing);
basicResampler->SetOutputSize(size);
basicResampler->SetOutputStartIndex(start);
basicResampler->SetEdgePaddingValue(defaultValue);
......
......@@ -25,11 +25,11 @@
#include "otbCompositeTransform.h"
#include "itkScalableAffineTransform.h"
#include "itkTranslationTransform.h"
#include "itkIdentityTransform.h"
#include "itkScaleTransform.h"
#include "itkCenteredRigid2DTransform.h"
#include "otbStreamingResampleImageFilter.h"
#include "otbGridResampleImageFilter.h"
namespace otb
{
......@@ -64,8 +64,7 @@ public:
typedef itk::TranslationTransform<double, FloatVectorImageType::ImageDimension> TransformType;
typedef otb::StreamingResampleImageFilter<FloatVectorImageType, FloatVectorImageType, double> ResampleFilterType;
typedef itk::IdentityTransform<double, FloatVectorImageType::ImageDimension> IdentityTransformType;
typedef otb::GridResampleImageFilter<FloatVectorImageType,FloatVectorImageType> GridResampleFilterType;
typedef itk::ScalableAffineTransform<double, FloatVectorImageType::ImageDimension> ScalableTransformType;
typedef ScalableTransformType::OutputVectorType OutputVectorType;
......@@ -182,7 +181,9 @@ private:
FloatVectorImageType* inputImage = GetParameterImage("in");
m_Resampler = ResampleFilterType::New();
m_GridResampler = GridResampleFilterType::New();
m_Resampler->SetInput(inputImage);
m_GridResampler->SetInput(inputImage);
// Get Interpolator
switch ( GetParameterInt("interpolator") )
......@@ -218,9 +219,7 @@ private:
{
case Transform_Identity:
{
IdentityTransformType::Pointer transform = IdentityTransformType::New();
m_Resampler->SetOutputParametersFromImage( inputImage );
m_GridResampler->SetOutputParametersFromImage( inputImage );
// Scale Transform
OutputVectorType scale;
scale[0] = 1.0 / GetParameterFloat("transform.type.id.scalex");
......@@ -232,23 +231,21 @@ private:
OutputSpacing[0] = spacing[0] * scale[0];
OutputSpacing[1] = spacing[1] * scale[1];
m_Resampler->SetOutputSpacing(OutputSpacing);
m_GridResampler->SetOutputSpacing(OutputSpacing);
FloatVectorImageType::PointType origin = inputImage->GetOrigin();
FloatVectorImageType::PointType outputOrigin;
outputOrigin[0] = origin[0] + 0.5 * spacing[0] * (scale[0] - 1.0);
outputOrigin[1] = origin[1] + 0.5 * spacing[1] * (scale[1] - 1.0);
m_Resampler->SetOutputOrigin(outputOrigin);
m_Resampler->SetTransform(transform);
m_GridResampler->SetOutputOrigin(outputOrigin);
// Evaluate size
ResampleFilterType::SizeType recomputedSize;
recomputedSize[0] = inputImage->GetLargestPossibleRegion().GetSize()[0] / scale[0];
recomputedSize[1] = inputImage->GetLargestPossibleRegion().GetSize()[1] / scale[1];
m_Resampler->SetOutputSize(recomputedSize);
m_GridResampler->SetOutputSize(recomputedSize);
otbAppLogINFO( << "Output image size : " << recomputedSize );
}
break;
......@@ -437,6 +434,8 @@ private:
}
ResampleFilterType::Pointer m_Resampler;
GridResampleFilterType::Pointer m_GridResampler;
}; //class
......
......@@ -18,6 +18,7 @@
#include "otbWrapperApplicationFactory.h"
#include "otbGenericRSResampleImageFilter.h"
#include "otbGridResampleImageFilter.h"
#include "otbImportGeoInformationImageFilter.h"
#include "otbBCOInterpolateImageFunction.h"
......@@ -72,7 +73,7 @@ public:
typedef itk::ScalableAffineTransform<double, 2> TransformType;
typedef otb::StreamingResampleImageFilter
typedef otb::GridResampleImageFilter
<FloatVectorImageType,
FloatVectorImageType> BasicResamplerType;
......@@ -251,19 +252,23 @@ private:
{
otbAppLogINFO("Using the PHR mode");
otb::PleiadesPToXSAffineTransformCalculator::TransformType::Pointer transform
= otb::PleiadesPToXSAffineTransformCalculator::Compute(GetParameterImage("inr"),
GetParameterImage("inm"));
m_BasicResampler->SetTransform(transform);
otb::PleiadesPToXSAffineTransformCalculator::TransformType::OffsetType offset
= otb::PleiadesPToXSAffineTransformCalculator::ComputeOffset(GetParameterImage("inr"),
GetParameterImage("inm"));
m_BasicResampler->SetInput(movingImage);
origin+=offset;
origin[0]=origin[0]/4;
origin[1]=origin[1]/4;
m_BasicResampler->SetOutputOrigin(origin);
m_BasicResampler->SetOutputSpacing(spacing);
m_BasicResampler->SetOutputSize(size);
m_BasicResampler->SetOutputStartIndex(start);
FloatVectorImageType::SpacingType xsSpacing = GetParameterImage("inm")->GetSpacing();
xsSpacing*=0.25;
m_BasicResampler->SetOutputSpacing(xsSpacing);
m_BasicResampler->SetOutputSize(size);
m_Resampler->SetOutputStartIndex(start);
m_BasicResampler->SetEdgePaddingValue(defaultValue);
m_GeoImport->SetInput(m_BasicResampler->GetOutput());
......
/*=========================================================================
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 __otbGridResampleImageFilter_h
#define __otbGridResampleImageFilter_h
#include "itkImageToImageFilter.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkDefaultConvertPixelTraits.h"
#include "otbMacro.h"
namespace otb
{
/** \class GridResampleImageFilter
* \brief Resample input image on a new origin/spacing/size grid
*
* This filter resample the input image on a new grid defined by
* OutputOrigin, OutputSpacing and OutputSize, using the provided
* interpolator.
*
* This is equivalent to a itk::ResampleImageFilter using an
* itk::IdentityTransform, except that it is slightly more efficient
* and that in this simplified case it is possible to explicitely
* compute the input requested region. The GridResampleImageFilter
* therefore supports streaming, contrary to the
* itk::ResampleImageFilter.
*
* When grid position is outside of the input image domain, the
* default EdgePaddingValue is used.
*
* If CheckOutputBounds flag is set to true (default value), the
* interpolated value will be checked for output pixel type range
* prior to casting.
*
* \ingroup OTBImageManipulation
* \ingroup Streamed
* \ingroup Threaded
**/
template <typename TInputImage, typename TOutputImage,
typename TInterpolatorPrecision = double>
class ITK_EXPORT GridResampleImageFilter :
public itk::ImageToImageFilter<TInputImage, TOutputImage>
{
public:
/** Standard class typedefs. */
typedef GridResampleImageFilter 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);
/** Run-time type information (and related methods). */
itkTypeMacro(GridResampleImageFilter, itk::ImageToImageFilter);
/** Number of dimensions. */
itkStaticConstMacro(ImageDimension, unsigned int,
TOutputImage::ImageDimension);
itkStaticConstMacro(InputImageDimension, unsigned int,
TInputImage::ImageDimension);
/** Typedef parameters*/
typedef TInputImage InputImageType;
typedef TOutputImage OutputImageType;
typedef typename OutputImageType::RegionType OutputImageRegionType;
typedef typename TOutputImage::PixelType OutputPixelType;
typedef itk::DefaultConvertPixelTraits<OutputPixelType> OutputPixelConvertType;
typedef typename OutputPixelConvertType::ComponentType OutputPixelComponentType;
/** Interpolator type */
typedef itk::InterpolateImageFunction<InputImageType,
TInterpolatorPrecision> InterpolatorType;
typedef typename InterpolatorType::Pointer InterpolatorPointerType;
typedef itk::LinearInterpolateImageFunction<InputImageType,
TInterpolatorPrecision> DefaultInterpolatorType;
typedef typename InterpolatorType::OutputType InterpolatorOutputType;
typedef itk::DefaultConvertPixelTraits< InterpolatorOutputType > InterpolatorConvertType;
typedef typename InterpolatorConvertType::ComponentType InterpolatorComponentType;
/** Input pixel continuous index typdef */
typedef typename itk::ContinuousIndex<double,InputImageDimension > ContinuousInputIndexType;
/** ImageBase typedef */
typedef itk::ImageBase<OutputImageType::ImageDimension> ImageBaseType;
typedef typename ImageBaseType::SpacingType SpacingType;
typedef typename ImageBaseType::SizeType SizeType;
typedef typename ImageBaseType::PointType PointType;
typedef typename ImageBaseType::IndexType IndexType;
itkSetMacro(OutputStartIndex,IndexType);
itkGetConstReferenceMacro(OutputStartIndex,IndexType);
itkSetMacro(OutputSize,SizeType);
itkGetConstReferenceMacro(OutputSize,SizeType);
itkSetMacro(OutputOrigin,PointType);
itkGetConstReferenceMacro(OutputOrigin,PointType);
itkSetMacro(OutputSpacing,SpacingType);
itkGetConstReferenceMacro(OutputSpacing,SpacingType);
itkSetMacro(EdgePaddingValue,OutputPixelType);
itkGetConstReferenceMacro(EdgePaddingValue,OutputPixelType);
itkSetMacro(CheckOutputBounds,bool);
itkGetMacro(CheckOutputBounds,bool);
itkBooleanMacro(CheckOutputBounds);
itkSetObjectMacro(Interpolator,InterpolatorType);
itkGetObjectMacro(Interpolator,InterpolatorType);
/** Import output parameters from a given image */
void SetOutputParametersFromImage(const ImageBaseType * image);
/** Method Compute the Modified Time based on changed to the components. */
itk::ModifiedTimeType GetMTime(void) const;
protected:
GridResampleImageFilter();
/** Destructor */
virtual ~GridResampleImageFilter() {};
virtual void GenerateOutputInformation();
virtual void GenerateInputRequestedRegion();
virtual void BeforeThreadedGenerateData();
virtual void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId);
virtual void AfterThreadedGenerateData();
inline void CastPixelWithBoundsChecking( const InterpolatorOutputType& value,
const InterpolatorComponentType& minComponent,
const InterpolatorComponentType& maxComponent,
OutputPixelType& outputValue ) const
{
// Method imported from itk::ResampleImageFilter
const unsigned int nComponents = InterpolatorConvertType::GetNumberOfComponents(value);
itk::NumericTraits<OutputPixelType>::SetLength( outputValue, nComponents );
for (unsigned int n=0; n<nComponents; n++)
{
InterpolatorComponentType component = InterpolatorConvertType::GetNthComponent( n, value );
if ( m_CheckOutputBounds && component < minComponent )
{
OutputPixelConvertType::SetNthComponent( n, outputValue, static_cast<OutputPixelComponentType>( minComponent ) );
}
else if ( m_CheckOutputBounds && component > maxComponent )
{
OutputPixelConvertType::SetNthComponent( n, outputValue, static_cast<OutputPixelComponentType>( maxComponent ) );
}
else
{
OutputPixelConvertType::SetNthComponent(n, outputValue,
static_cast<OutputPixelComponentType>( component ) );
}
}
}
void PrintSelf(std::ostream& os, itk::Indent indent) const;
private:
GridResampleImageFilter(const Self &); //purposely not implemented
void operator =(const Self&); //purposely not implemented
IndexType m_OutputStartIndex; // output image start index
SizeType m_OutputSize; // output image size
PointType m_OutputOrigin; // output image origin
SpacingType m_OutputSpacing; // output image spacing
OutputPixelType m_EdgePaddingValue; // Default pixel value
bool m_CheckOutputBounds; // Shall we check
// output bounds when
// casting?
InterpolatorPointerType m_Interpolator; // Interpolator used
// for resampling
OutputImageRegionType m_ReachableOutputRegion; // Internal
// variable for
// speed-up. Computed
// in BeforeThreadedGenerateData
};
} // namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
#include "otbGridResampleImageFilter.txx"
#endif
#endif
/*=========================================================================
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 __otbGridResampleImageFilter_txx
#define __otbGridResampleImageFilter_txx
#include "otbGridResampleImageFilter.h"
#include "otbStreamingTraits.h"
#include "itkNumericTraits.h"
#include "itkProgressReporter.h"
#include "itkImageScanlineIterator.h"
#include "itkContinuousIndex.h"
namespace otb
{
template <typename TInputImage, typename TOutputImage,
typename TInterpolatorPrecision>
GridResampleImageFilter<TInputImage, TOutputImage, TInterpolatorPrecision>
::GridResampleImageFilter()
: m_OutputStartIndex(),
m_OutputSize(),
m_OutputOrigin(),
m_OutputSpacing(),
m_EdgePaddingValue(),
m_CheckOutputBounds(true),
m_Interpolator(),
m_ReachableOutputRegion()
{
// Set linear interpolator as default
m_Interpolator = dynamic_cast<InterpolatorType *>(DefaultInterpolatorType::New().GetPointer());
// Initialize EdgePaddingValue
m_EdgePaddingValue = itk::NumericTraits<OutputPixelType>::ZeroValue(m_EdgePaddingValue);
// Initialize origin and spacing
m_OutputOrigin.Fill(0.);
m_OutputSpacing.Fill(1.);
m_OutputStartIndex.Fill(0);
m_OutputSize.Fill(0);
}
/** Import output parameters from a given image */
template <typename TInputImage, typename TOutputImage,
typename TInterpolatorPrecision>
void
GridResampleImageFilter<TInputImage, TOutputImage, TInterpolatorPrecision>
::SetOutputParametersFromImage(const ImageBaseType * image)
{
this->SetOutputOrigin ( image->GetOrigin() );
this->SetOutputSpacing ( image->GetSpacing() );
this->SetOutputStartIndex ( image->GetLargestPossibleRegion().GetIndex() );
this->SetOutputSize ( image->GetLargestPossibleRegion().GetSize() );
}
template <typename TInputImage, typename TOutputImage,
typename TInterpolatorPrecision>
void
GridResampleImageFilter<TInputImage, TOutputImage, TInterpolatorPrecision>
::GenerateOutputInformation()
{
// call the superclass' implementation of this method
Superclass::GenerateOutputInformation();
// get pointers to the input and output
OutputImageType *outputPtr = this->GetOutput();
if ( !outputPtr )
{
return;
}
// Fill output image data
typename TOutputImage::RegionType outputLargestPossibleRegion;
outputLargestPossibleRegion.SetSize(m_OutputSize);
outputLargestPossibleRegion.SetIndex(m_OutputStartIndex);
outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion);
outputPtr->SetSpacing(m_OutputSpacing);
outputPtr->SetOrigin(m_OutputOrigin);
// TODO: Report no data value here
}
template <typename TInputImage, typename TOutputImage,
typename TInterpolatorPrecision>
void
GridResampleImageFilter<TInputImage, TOutputImage, TInterpolatorPrecision>
::GenerateInputRequestedRegion()
{
// call the superclass's implementation of this method
Superclass::GenerateInputRequestedRegion();
// Check for input
if ( !this->GetInput() )
{
return;
}
// get pointers to the input and output
typename TInputImage::Pointer inputPtr =
const_cast< TInputImage * >( this->GetInput() );
// check for output
OutputImageType *outputPtr = this->GetOutput();
if ( !outputPtr )
{
return;
}
typename TOutputImage::RegionType outputRequestedRegion = outputPtr->GetRequestedRegion();
IndexType outULIndex,outLRIndex;
typedef itk::ContinuousIndex<double,TInputImage::ImageDimension> ContinuousIndexType;
ContinuousIndexType inULCIndex,inLRCIndex;
// Get output image requested region corners as Index
outULIndex = outputRequestedRegion.GetIndex();
outLRIndex = outULIndex+outputRequestedRegion.GetSize();
outLRIndex[0]-=1;
outLRIndex[1]-=1;
// Transform to physiscal points
PointType outULPoint,outLRPoint;
outputPtr->TransformIndexToPhysicalPoint(outULIndex,outULPoint);
outputPtr->TransformIndexToPhysicalPoint(outLRIndex,outLRPoint);
// Transform to input image Index
inputPtr->TransformPhysicalPointToContinuousIndex(outULPoint,inULCIndex);
inputPtr->TransformPhysicalPointToContinuousIndex(outLRPoint,inLRCIndex);
SizeType inSize;
IndexType inULIndex,inLRIndex;
// Reorder index properly and compute size
for(unsigned int dim = 0; dim < ImageDimension;++dim)
{
if(inULCIndex[dim] > inLRCIndex[dim])
{
// swap
typename ContinuousIndexType::ValueType tmp(inULCIndex[dim]);
inULCIndex[dim]=inLRCIndex[dim];
inLRCIndex[dim]=tmp;
}
// Ensure correct rounding of coordinates
inULIndex[dim] = vcl_floor(inULCIndex[dim]);
inLRIndex[dim] = vcl_ceil(inLRCIndex[dim]);
inSize[dim] = static_cast<typename SizeType::SizeValueType>(inLRIndex[dim]-inULIndex[dim])+1;
}
// Build the input requested region
typename TInputImage::RegionType inputRequestedRegion;
inputRequestedRegion.SetIndex(inULIndex);
inputRequestedRegion.SetSize(inSize);
// Compute the padding due to the interpolator
unsigned int interpolatorRadius =
StreamingTraits<typename Superclass::InputImageType>::CalculateNeededRadiusForInterpolator(this->GetInterpolator());
inputRequestedRegion.PadByRadius(interpolatorRadius);
// crop the input requested region at the input's largest possible region
if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
{
inputPtr->SetRequestedRegion(inputRequestedRegion);
}
else
{
// store what we tried to request (prior to trying to crop)
inputPtr->SetRequestedRegion(inputRequestedRegion);
// build an exception
itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
e.SetLocation(ITK_LOCATION);
e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
e.SetDataObject(inputPtr);
throw e;
}
}
template <typename TInputImage, typename TOutputImage,