Commit bdc3e82a authored by Julien Michel's avatar Julien Michel
Browse files

Ajout de l'interpolation linéaire et de l'interpolation par spline de la déformation.

parent d964f4bb
/*=========================================================================
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 _otbBSplinesInterpolateDeformationFieldGenerator_h
#define _otbBSplinesInterpolateDeformationFieldGenerator_h
#include "otbPointSetToDeformationFieldGenerator.h"
namespace otb
{
/** \class BSplinesInterpolateDeformationFieldGenerator
* \brief This filters encapsulate the itk::DeformationFieldSource to produce a BSpline interpolation of the point in point set whose metric values are
* sufficient.
* \sa itk::DeformationFieldSource
* \ingroup
* \ingroup
*/
template <class TPointSet, class TDeformationField>
class ITK_EXPORT BSplinesInterpolateDeformationFieldGenerator
: public PointSetToDeformationFieldGenerator<TPointSet, TDeformationField>
{
public:
/** Standard typedefs */
typedef BSplinesInterpolateDeformationFieldGenerator Self;
typedef PointSetToDeformationFieldGenerator<TPointSet,TDeformationField> Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
/** Type macro */
itkNewMacro(Self);
/** Creation through object factory macro */
itkTypeMacro(BSplinesInterpolateDeformationFieldGenerator,PointSetToDeformationFieldGenerator);
/** Template parameters typedefs */
typedef typename Superclass::PointSetType PointSetType;
typedef typename Superclass::PointSetPointerType PointSetPointerType;
typedef typename Superclass::DeformationFieldType DeformationFieldType;
typedef typename Superclass::DeformationFieldPointerType DeformationFieldPointerType;
typedef typename Superclass::IndexType IndexType;
typedef typename DeformationFieldType::PixelType PixelType;
typedef typename Superclass::ValueType ValueType;
typedef typename Superclass::PointType PointType;
typedef typename Superclass::IndexVectorType IndexVectorType;
typedef typename Superclass::DistanceVectorType DistanceVectorType;
protected:
/** Constructor */
BSplinesInterpolateDeformationFieldGenerator() {};
/** Destructor */
virtual ~BSplinesInterpolateDeformationFieldGenerator() {};
/**PrintSelf method */
virtual void PrintSelf(std::ostream& os, itk::Indent indent) const;
/** Main computation method */
virtual void GenerateData();
private:
BSplinesInterpolateDeformationFieldGenerator(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
};
}// End namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
#include "otbBSplinesInterpolateDeformationFieldGenerator.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 _otbBSplinesInterpolateDeformationFieldGenerator_txx
#define _otbBSplinesInterpolateDeformationFieldGenerator_txx
#include "otbBSplinesInterpolateDeformationFieldGenerator.h"
#include "itkDeformationFieldSource.h"
#include "otbMacro.h"
#include "otbImage.h"
#include "itkImageRegionIterator.h"
namespace otb
{
/** Main computation method */
template <class TPointSet,class TDeformationField>
void
BSplinesInterpolateDeformationFieldGenerator<TPointSet, TDeformationField>
::GenerateData(void)
{
DeformationFieldPointerType outputPtr = this->GetOutput();
typedef itk::Vector< ValueType,2 > VectorType;
typedef otb::Image< VectorType,DeformationFieldType::ImageDimension> ImageType;
typedef itk::DeformationFieldSource<ImageType> DeformationFieldSourceType;
typedef typename DeformationFieldSourceType::LandmarkContainerPointer LandmarkContainerPointer;
typedef typename DeformationFieldSourceType::LandmarkContainer LandmarkContainerType;
typedef typename DeformationFieldSourceType::LandmarkPointType LandmarkPointType;
typedef typename PointSetType::PointsContainer PointsContainer;
typedef typename PointsContainer::ConstIterator PointsIterator;
typedef typename PointSetType::PointDataContainer PointDataContainer;
typedef typename PointDataContainer::ConstIterator PointDataIterator;
typename LandmarkContainerType::Pointer sourceLandmarks = LandmarkContainerType::New();
typename LandmarkContainerType::Pointer targetLandmarks = LandmarkContainerType::New();
LandmarkPointType sourcePoint;
LandmarkPointType targetPoint;
PointsIterator pointIterator = this->GetPointSet()->GetPoints()->Begin();
PointsIterator end = this->GetPointSet()->GetPoints()->End();
unsigned int pointId = 0;
PointDataIterator pointDataIterator = this->GetPointSet()->GetPointData()->Begin();
while( pointIterator != end )
{
typename PointDataContainer::Element valueAndDeformations = pointDataIterator.Value();
if(vcl_abs(valueAndDeformations[0])>=this->GetMetricThreshold())
{
typename PointSetType::PointType p = pointIterator.Value(); // access the point
sourcePoint[0] = p[0];
sourcePoint[1] = p[1];
targetPoint[0] = p[0] + valueAndDeformations[1];
targetPoint[1] = p[1] + valueAndDeformations[2];
otbMsgDebugMacro(<<"Adding landmark "<<pointId<<", source point: "<<sourcePoint<<", targetpoint: "<<targetPoint);
sourceLandmarks->InsertElement( pointId, sourcePoint );
targetLandmarks->InsertElement( pointId, targetPoint );
++pointId;
}
++pointIterator;
++pointDataIterator;
}
typename DeformationFieldSourceType::Pointer deformer = DeformationFieldSourceType::New();
deformer->SetOutputSpacing(this->GetOutputSpacing());
deformer->SetOutputOrigin(this->GetOutputOrigin());
deformer->SetOutputRegion(outputPtr->GetRequestedRegion());
deformer->SetSourceLandmarks(sourceLandmarks.GetPointer());
deformer->SetTargetLandmarks(targetLandmarks.GetPointer());
deformer->Update();
outputPtr->Allocate();
PixelType defaultPixel;
defaultPixel.SetSize(2);
defaultPixel.Fill(this->GetDefaultValue());
outputPtr->FillBuffer(defaultPixel);
typedef itk::ImageRegionIterator<ImageType> ImageIteratorType;
typedef itk::ImageRegionIterator<DeformationFieldType> OutputIteratorType;
ImageIteratorType inIt(deformer->GetOutput(),outputPtr->GetRequestedRegion());
OutputIteratorType outIt(outputPtr,outputPtr->GetRequestedRegion());
int i=0;
// Casting otb::Image<itt::Vector<ValueType,2>,2> to otb::VectorImage<ValueType,2>
for(inIt.GoToBegin(),outIt.GoToBegin();(!inIt.IsAtEnd())&&(!outIt.IsAtEnd());++inIt,++outIt,++i)
{
typename ImageType::PixelType inPixel;
inPixel = inIt.Get();
PixelType outPixel;
outPixel.SetSize(2);
outPixel[0]=inPixel[0];
outPixel[1]=inPixel[1];
outIt.Set(outPixel);
}
}
/**
* PrintSelf Method
*/
template <class TPointSet,class TDeformationField>
void
BSplinesInterpolateDeformationFieldGenerator<TPointSet, TDeformationField>
::PrintSelf(std::ostream& os, itk::Indent indent) const
{
Superclass::PrintSelf(os, indent);
}
} // End namespace otb
#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 _otbNNearestPointsLinearInterpolateDeformationFieldGenerator_h
#define _otbNNearestPointsLinearInterpolateDeformationFieldGenerator_h
#include "otbPointSetToDeformationFieldGenerator.h"
namespace otb
{
/** \class NNearestPointsLinearInterpolateDeformationFieldGenerator
* \brief This class generate the deformation field by performing a linear interpolation of the deformation induced by the n nearest point.
* \ingroup
* \ingroup
*/
template <class TPointSet, class TDeformationField>
class ITK_EXPORT NNearestPointsLinearInterpolateDeformationFieldGenerator
: public PointSetToDeformationFieldGenerator<TPointSet, TDeformationField>
{
public:
/** Standard typedefs */
typedef NNearestPointsLinearInterpolateDeformationFieldGenerator Self;
typedef PointSetToDeformationFieldGenerator<TPointSet,TDeformationField> Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
/** Type macro */
itkNewMacro(Self);
/** Creation through object factory macro */
itkTypeMacro(NNearestPointsLinearInterpolateDeformationFieldGenerator,PointSetToDeformationFieldGenerator);
/** Template parameters typedefs */
typedef typename Superclass::PointSetType PointSetType;
typedef typename Superclass::PointSetPointerType PointSetPointerType;
typedef typename Superclass::DeformationFieldType DeformationFieldType;
typedef typename Superclass::DeformationFieldPointerType DeformationFieldPointerType;
typedef typename Superclass::IndexType IndexType;
typedef typename DeformationFieldType::PixelType PixelType;
typedef typename Superclass::ValueType ValueType;
typedef typename Superclass::PointType PointType;
typedef typename Superclass::IndexVectorType IndexVectorType;
typedef typename Superclass::DistanceVectorType DistanceVectorType;
itkSetMacro(NumberOfPoints,unsigned int);
itkGetMacro(NumberOfPoints,unsigned int);
protected:
/** Constructor */
NNearestPointsLinearInterpolateDeformationFieldGenerator() {};
/** Destructor */
virtual ~NNearestPointsLinearInterpolateDeformationFieldGenerator() {};
/**PrintSelf method */
virtual void PrintSelf(std::ostream& os, itk::Indent indent) const;
/** Main computation method */
virtual void GenerateData();
private:
NNearestPointsLinearInterpolateDeformationFieldGenerator(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
unsigned int m_NumberOfPoints;
};
}// End namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
#include "otbNNearestPointsLinearInterpolateDeformationFieldGenerator.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 _otbNNearestPointsLinearInterpolateDeformationFieldGenerator_txx
#define _otbNNearestPointsLinearInterpolateDeformationFieldGenerator_txx
#define EPSILON 1e-15
#include "otbNNearestPointsLinearInterpolateDeformationFieldGenerator.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "otbMacro.h"
namespace otb
{
/** Main computation method */
template <class TPointSet,class TDeformationField>
void
NNearestPointsLinearInterpolateDeformationFieldGenerator<TPointSet, TDeformationField>
::GenerateData(void)
{
DeformationFieldPointerType outputPtr = this->GetOutput();
PixelType defaultValue(2);
defaultValue.Fill(this->GetDefaultValue());
outputPtr->Allocate();
outputPtr->FillBuffer(defaultValue);
typedef itk::ImageRegionIteratorWithIndex<DeformationFieldType> IteratorType;
IteratorType it(outputPtr,outputPtr->GetRequestedRegion());
for(it.GoToBegin();!it.IsAtEnd();++it)
{
IndexType index = it.GetIndex();
IndexVectorType indexVector = this->GenerateNearestValidPointsPointSet(it.GetIndex(),m_NumberOfPoints);
PixelType pixel(2);
double xdeformation, ydeformation,normalization;
xdeformation = 0;
ydeformation = 0;
normalization = 0;
for(typename IndexVectorType::iterator indexIt=indexVector.begin();indexIt!=indexVector.end();++indexIt)
{
PointType point;
point[0] = static_cast<double>(this->GetPointSet()->GetPoints()->GetElement(*indexIt)[0]);
point[1] = static_cast<double>(this->GetPointSet()->GetPoints()->GetElement(*indexIt)[1]);
double distance = this->EuclideanDistance(index,point);
if(distance<EPSILON)
{
distance = EPSILON;
}
xdeformation += this->GetPointSet()->GetPointData()->GetElement((*indexIt))[1]/distance;
ydeformation += this->GetPointSet()->GetPointData()->GetElement((*indexIt))[2]/distance;
normalization+=1/distance;
}
if(normalization>0)
{
pixel[0] = static_cast<ValueType>(xdeformation/normalization);
pixel[1] = static_cast<ValueType>(ydeformation/normalization);
}
else
{
pixel=defaultValue;
}
it.Set(pixel);
}
}
/**
* PrintSelf Method
*/
template <class TPointSet,class TDeformationField>
void
NNearestPointsLinearInterpolateDeformationFieldGenerator<TPointSet, TDeformationField>
::PrintSelf(std::ostream& os, itk::Indent indent) const
{
Superclass::PrintSelf(os, indent);
}
} // End namespace otb
#endif
......@@ -132,12 +132,13 @@ protected:
*/
IndexVectorType GenerateNearestValidPointsPointSet(IndexType index, unsigned int n = 1);
/** Euclidean distance of point to index */
double EuclideanDistance(IndexType index, PointType p);
private:
PointSetToDeformationFieldGenerator(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
/** Euclidean distance of point to index */
double EuclideanDistance(IndexType index, PointType p);
/**
* The threshold of metric value.
*/
......
......@@ -50,6 +50,33 @@ ADD_TEST(dmTvNearestPointDeformationFieldGenerator ${DISPARITYMAP_TESTS}
${TEMP}/dmTvNearestPointDeformationField.hdr
)
# ------- otb::NNearestPointsLinearInterpolateDeformationFieldGenerator ----------
ADD_TEST(dmTuNNearestPointsLinearInterpolateDeformationFieldGeneratorNew ${DISPARITYMAP_TESTS}
otbNNearestPointsLinearInterpolateDeformationFieldGeneratorNew)
ADD_TEST(dmTvNNearestPointsLinearInterpolateDeformationFieldGenerator ${DISPARITYMAP_TESTS}
--compare-image ${TOL}
${BASELINE}/dmTvNNearestPointsLinearInterpolateDeformationField.hdr
${TEMP}/dmTvNNearestPointsLinearInterpolateDeformationField.hdr
otbNNearestPointsLinearInterpolateDeformationFieldGenerator
${TEMP}/dmTvNNearestPointsLinearInterpolateDeformationField.hdr
)
# ------- otb::BSplinesInterpolateDeformationFieldGenerator ----------
ADD_TEST(dmTuBSplinesInterpolateDeformationFieldGeneratorNew ${DISPARITYMAP_TESTS}
otbNNearestPointsLinearInterpolateDeformationFieldGeneratorNew)
ADD_TEST(dmTvBSplinesInterpolateDeformationFieldGenerator ${DISPARITYMAP_TESTS}
--compare-image ${TOL}
${BASELINE}/dmTvBSplinesInterpolateDeformationField.hdr
${TEMP}/dmTvBSplinesInterpolateDeformationField.hdr
otbBSplinesInterpolateDeformationFieldGenerator
${TEMP}/dmTvBSplinesInterpolateDeformationField.hdr
)
# ------- Fichiers sources CXX -----------------------------------
SET(BasicDisparityMap_SRCS
otbDisparityMapEstimationMethodNew.cxx
......@@ -57,6 +84,11 @@ otbDisparityMapEstimationMethod.cxx
otbPointSetToDeformationFieldGeneratorNew.cxx
otbNearestPointDeformationFieldGeneratorNew.cxx
otbNearestPointDeformationFieldGenerator.cxx
otbNNearestPointsLinearInterpolateDeformationFieldGeneratorNew.cxx
otbNNearestPointsLinearInterpolateDeformationFieldGenerator.cxx
otbBSplinesInterpolateDeformationFieldGeneratorNew.cxx
otbBSplinesInterpolateDeformationFieldGenerator.cxx
)
INCLUDE_DIRECTORIES("${OTBTesting_BINARY_DIR}")
......
/*=========================================================================
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.
=========================================================================*/
#include "itkExceptionObject.h"
#include "itkPointSet.h"
#include "otbVectorImage.h"
#include "otbBSplinesInterpolateDeformationFieldGenerator.h"
#include "otbImageFileWriter.h"
int otbBSplinesInterpolateDeformationFieldGenerator(int argc, char * argv[])
{
try
{
const unsigned int Dimension = 2;
const char * outfname = argv[1];
typedef double PixelType;
typedef otb::VectorImage<PixelType,Dimension> ImageType;
typedef itk::Array<double> ParamType;
typedef itk::PointSet<ParamType,Dimension> PointSetType;
typedef PointSetType::PointType PointType;
typedef otb::BSplinesInterpolateDeformationFieldGenerator<PointSetType,ImageType> FilterType;
typedef otb::ImageFileWriter<ImageType> WriterType;
ImageType::SizeType size;
size.Fill(100);
double thresh = 0.9;
// Preparing point set
PointSetType::Pointer ps = PointSetType::New();
PointType p1,p2,p3,p4,p5;
ParamType pd1(3),pd2(3),pd3(3),pd4(3),pd5(3);
p1[0] = 10;
p1[1] = 10;
p2[0] = 75;
p2[1] = 10;
p3[0] = 50;
p3[1] = 50;
p4[0] = 10;
p4[1] = 60;
p5[0] = 85;
p5[1] = 70;
pd1[0] = 0.95;
pd1[1] = 10;
pd1[2] = -5;
pd2[0] = 0.98;
pd2[1] = 2;
pd2[2] = 5;
pd3[0] = 0.5;
pd3[1] = 20;
pd3[2] = -20;
pd4[0] = 0.91;
pd4[1] = 15;
pd4[2] = -5;
pd5[0] = 0.91;
pd5[1] = 5;
pd5[2] = 5;
ps->SetPoint(0,p1);
ps->SetPointData(0,pd1);
ps->SetPoint(1,p2);
ps->SetPointData(1,pd2);
ps->SetPoint(2,p3);
ps->SetPointData(2,pd3);
ps->SetPoint(3,p4);
ps->SetPointData(3,pd4);
ps->SetPoint(4,p5);
ps->SetPointData(4,pd5);
// Instantiating object
FilterType::Pointer filter = FilterType::New();
filter->SetOutputSize(size);
filter->SetMetricThreshold(thresh);
filter->SetPointSet(ps);
WriterType::Pointer writer = WriterType::New();
writer->SetInput(filter->GetOutput());
writer->SetFileName(outfname);
writer->Update();
}
catch( itk::ExceptionObject & err )
{
std::cout << "Exception itk::ExceptionObject thrown !" << std::endl;
std::cout << err << std::endl;
return EXIT_FAILURE;
}
catch( ... )
{
std::cout << "Unknown exception thrown !" << std::endl;
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
/*=========================================================================
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.
=========================================================================*/
#include "itkExceptionObject.h"
#include "itkPointSet.h"
#include "otbVectorImage.h"
#include "otbBSplinesInterpolateDeformationFieldGenerator.h"
int otbBSplinesInterpolateDeformationFieldGeneratorNew(int argc, char * argv[])
{
try
{
const unsigned int Dimension = 2;
typedef double PixelType;
typedef otb::VectorImage<PixelType,Dimension> ImageType;
typedef ImageType::PointType PointType;
typedef itk::PointSet<PointType,Dimension> PointSetType;
typedef otb::BSplinesInterpolateDeformationFieldGenerator<PointSetType,ImageType> FilterType;
// Instantiating object
FilterType::Pointer filter = FilterType::New();
}