Commit 032cac4e authored by Cyrille Valladeau's avatar Cyrille Valladeau

ENH : add FineCorrelation class + tests

parent b5b894e7
/*=========================================================================
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 __otbFineCorrelationImageFilter_h
#define __otbFineCorrelationImageFilter_h
#include "itkImageToImageFilter.h"
#include "itkConstNeighborhoodIterator.h"
#include "itkImageRegionIterator.h"
#include "itkInterpolateImageFunction.h"
#include "itkContinuousIndex.h"
namespace otb
{
/** \class FineCorrelationImageFilter
* \brief Computes a deformation field between two images
*
* This class compute a correlation field and the associated
* deformation field between two images.
*
* The m_Radius parameter defines the size of the window to compute
* the normalized correlation.
*
* The m_SearchRadius parameter defines the size of the area where to
* search for a correlation peak in the moving image.
*
* Once the correlation peak has been found, it can be further refined
* by trying to fit a quadric surface in the premises of the
* maxima, thus obtaining sub-pixel precision. This option is
* activated by setting m_RefineLocation to true.
*
* TOutputCorrelation is supposed to be a scalar image with floating
* point precision.
*
* TOutputDeformationField is supposed to be an image whose pixel is a
* fixed size vector of size 2, with floating point precision.
*
* Images spacing is ignored during computation and in the output
* deformation field.
*
* All computation are done in double precision.
*
* \ingroup IntensityImageFilters Multithreaded, Streamed
*
*/
template <class TInputImage, class T0utputCorrelation, class TOutputDeformationField>
class ITK_EXPORT FineCorrelationImageFilter : public itk::ImageToImageFilter<TInputImage,T0utputCorrelation>
{
public:
/** Standard class typedefs. */
typedef FineCorrelationImageFilter Self;
typedef itk::ImageToImageFilter<TInputImage,T0utputCorrelation> 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(FineCorrelationImageFilter, ImageToImageFilter);
/** Some convenient typedefs. */
typedef typename T0utputCorrelation::RegionType OutputImageRegionType;
typedef typename TOutputDeformationField::PixelType DeformationValueType;
typedef typename TInputImage::Pointer InputImagePointerType;
typedef typename TInputImage::RegionType InputImageRegionType;
typedef typename TInputImage::SizeType SizeType;
typedef typename TInputImage::IndexType IndexType;
typedef typename itk::ImageRegionIterator<TInputImage> InputImageRegionIteratorType;
typedef typename itk::InterpolateImageFunction
<TInputImage, double> InterpolatorType;
typedef typename InterpolatorType::Pointer InterpolatorPointerType;
typedef typename itk::ContinuousIndex<double, 2> ContinuousIndexType;
/** Typedef for refinements mode */
typedef enum {COARSE = 0, LSQR_QUADFIT = 1, SUBPIXEL = 2} RefinementModeType;
/** Connect one of the operands for pixel-wise addition */
void SetFixedInput( const TInputImage * image);
/** Connect one of the operands for pixel-wise addition */
void SetMovingInput( const TInputImage * image);
/** Get the inputs */
const TInputImage * GetFixedInput();
const TInputImage * GetMovingInput();
/** Get the output deformation field */
TOutputDeformationField * GetOutputDeformationField();
/** Typedef of the input iterator */
typedef itk::ConstNeighborhoodIterator<TInputImage> NeighborhoodIteratorType;
typedef typename NeighborhoodIteratorType::RadiusType RadiusType;
typedef typename NeighborhoodIteratorType::NeighborhoodType NeighborhoodType;
typedef typename NeighborhoodIteratorType::OffsetType OffsetType;
itkSetMacro(Radius, RadiusType);
itkGetMacro(Radius, RadiusType);
itkSetMacro(SearchRadius,RadiusType);
itkGetMacro(SearchRadius,RadiusType);
/** Set the subpixel precision */
itkSetMacro(SubPixelPrecision,unsigned int);
itkGetMacro(SubPixelPrecision,unsigned int);
/** Set/Get Refinement mode */
itkSetEnumMacro(RefinementMode,RefinementModeType);
itkGetEnumMacro(RefinementMode,RefinementModeType);
void SetRefinementModeToCoarse()
{
m_RefinementMode = COARSE;
}
void SetRefinementModeToLSQRQuadFit()
{
m_RefinementMode = LSQR_QUADFIT;
}
void SetRefinementModeToSubPixel()
{
m_RefinementMode = SUBPIXEL;
}
/** Set unsigned int radius */
void SetRadius(unsigned int radius)
{
m_Radius.Fill(radius);
}
/** Set unsigned int radius */
void SetSearchRadius(unsigned int radius)
{
m_SearchRadius.Fill(radius);
}
/** Set the refinement mode */
protected:
/** Constructor */
FineCorrelationImageFilter();
/** Destructor */
virtual ~FineCorrelationImageFilter() {};
/** Threaded generate data */
virtual void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,int threadId );
/** Generate the input requested regions */
virtual void GenerateInputRequestedRegion(void);
/** Things to do before threaded computation */
virtual void BeforeThreadedGenerateData();
private:
FineCorrelationImageFilter(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
/** Compute correlation from offset */
inline double Correlation(const NeighborhoodType & fixed, const NeighborhoodType & moving, const OffsetType & offset, unsigned int scale) const;
/** Refine the location of the correlation maximum */
inline double RefineLocation(const NeighborhoodType & correlMap,const OffsetType & maxPos, DeformationValueType & value) const;
/** Compute subpixel neighborhood */
const NeighborhoodType ComputeSubPixelNeighborhood(const IndexType & location, unsigned int scale) const;
typedef std::pair<double,OffsetType> PairType;
typedef std::vector<PairType> PairVectorType;
/** Compare pairs by correlation values */
static bool CompPairs(const PairType &p1, const PairType &p2)
{
return p1.first>p2.first;
}
/** The radius for correlation */
RadiusType m_Radius;
/** The search radius */
RadiusType m_SearchRadius;
/** The interpolator */
InterpolatorPointerType m_Interpolator;
/** Choose the refinement type */
RefinementModeType m_RefinementMode;
/** Precision used in subpixel mode */
unsigned int m_SubPixelPrecision;
};
} // end namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
#include "otbFineCorrelationImageFilter.txx"
#endif
#endif
This diff is collapsed.
......@@ -34,6 +34,7 @@ SET(FEATUREEXTRACTION_TESTS12 ${CXX_TEST_PATH}/otbFeatureExtractionTests12)
SET(FEATUREEXTRACTION_TESTS13 ${CXX_TEST_PATH}/otbFeatureExtractionTests13)
SET(FEATUREEXTRACTION_TESTS14 ${CXX_TEST_PATH}/otbFeatureExtractionTests14)
SET(FEATUREEXTRACTION_TESTS15 ${CXX_TEST_PATH}/otbFeatureExtractionTests15)
SET(FEATUREEXTRACTION_TESTS16 ${CXX_TEST_PATH}/otbFeatureExtractionTests16)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ otbFeatureExtractionTests1 ~~~~~~~~~~~~~~~~~~~~~
......@@ -1889,6 +1890,86 @@ ADD_TEST(feTvScalarImageToPanTexTextureFilter ${FEATUREEXTRACTION_TESTS15}
${TEMP}/feTvScalarImageToPanTexTextureFilterOutput
8 5)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ otbFeatureExtractionTests16 ~~~~~~~~~~~~~~~~~~~~~
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
ADD_TEST(feTuFineCorrelationImageFilterNew ${FEATUREEXTRACTION_TESTS16}
otbFineCorrelationImageFilterNew
)
ADD_TEST(feTuFineCorrelationImageFilterTest_Coarse ${FEATUREEXTRACTION_TESTS16}
--compare-n-images ${EPSILON_8} 2
${BASELINE}/feTvFineCorrelationImageFilterTest_Coarse_Correl.tif
${TEMP}/feTvFineCorrelationImageFilterTest_Coarse_Correl.tif
${BASELINE}/feTvFineCorrelationImageFilterTest_Coarse_Field.tif
${TEMP}/feTvFineCorrelationImageFilterTest_Coarse_Field.tif
otbFineCorrelationImageFilterTest
${INPUTDATA}/ROI_IKO_PAN_LesHalles_sub.tif # fixedFileName
${INPUTDATA}/ROI_IKO_PAN_LesHalles_sub_warped_translation.tif # movingFileName
${TEMP}/feTvFineCorrelationImageFilterTest_Coarse_Correl.tif # output correlFileName
${TEMP}/feTvFineCorrelationImageFilterTest_Coarse_Field.tif # output fieldFileName
7 # radius
2 # sradius
0 # refine : Coarse
10 # precision
)
ADD_TEST(feTuFineCorrelationImageFilterTest_LSQR ${FEATUREEXTRACTION_TESTS16}
--compare-n-images ${EPSILON_8} 2
${BASELINE}/feTvFineCorrelationImageFilterTest_LSQR_Correl.tif
${TEMP}/feTvFineCorrelationImageFilterTest_LSQR_Correl.tif
${BASELINE}/feTvFineCorrelationImageFilterTest_LSQR_Field.tif
${TEMP}/feTvFineCorrelationImageFilterTest_LSQR_Field.tif
otbFineCorrelationImageFilterTest
${INPUTDATA}/ROI_IKO_PAN_LesHalles_sub.tif # fixedFileName
${INPUTDATA}/ROI_IKO_PAN_LesHalles_sub_warped_translation.tif # movingFileName
${TEMP}/feTvFineCorrelationImageFilterTest_LSQR_Correl.tif # output correlFileName
${TEMP}/feTvFineCorrelationImageFilterTest_LSQR_Field.tif # output fieldFileName
7 # radius
2 # sradius
1 # refine : LSQR
10 # precision
)
ADD_TEST(feTuFineCorrelationImageFilterTest_SUBPIX ${FEATUREEXTRACTION_TESTS16}
--compare-n-images ${EPSILON_8} 2
${BASELINE}/feTvFineCorrelationImageFilterTest_SUBPIX_Correl.tif
${TEMP}/feTvFineCorrelationImageFilterTest_SUBPIX_Correl.tif
${BASELINE}/feTvFineCorrelationImageFilterTest_SUBPIX_Field.tif
${TEMP}/feTvFineCorrelationImageFilterTest_SUBPIX_Field.tif
otbFineCorrelationImageFilterTest
${INPUTDATA}/ROI_IKO_PAN_LesHalles_sub.tif # fixedFileName
${INPUTDATA}/ROI_IKO_PAN_LesHalles_sub_warped_translation.tif # movingFileName
${TEMP}/feTvFineCorrelationImageFilterTest_SUBPIX_Correl.tif # output correlFileName
${TEMP}/feTvFineCorrelationImageFilterTest_SUBPIX_Field.tif # output fieldFileName
7 # radius
2 # sradius
2 # refine : SUBPIX
10 # precision
)
ADD_TEST(feTuFineCorrelationImageFilterTest_LSQR ${FEATUREEXTRACTION_TESTS16}
--compare-n-images ${EPSILON_8} 2
${BASELINE}/feTvFineCorrelationImageFilterTest_LSQR_Correl.tif
${TEMP}/feTvFineCorrelationImageFilterTest_LSQR_Correl.tif
${BASELINE}/feTvFineCorrelationImageFilterTest_LSQR_Field.tif
${TEMP}/feTvFineCorrelationImageFilterTest_LSQR_Field.tif
otbFineCorrelationImageFilterTest
${INPUTDATA}/ROI_IKO_PAN_LesHalles_sub.tif # fixedFileName
${INPUTDATA}/ROI_IKO_PAN_LesHalles_sub_warped_translation.tif # movingFileName
${TEMP}/feTvFineCorrelationImageFilterTest_LSQR_Correl.tif # output correlFileName
${TEMP}/feTvFineCorrelationImageFilterTest_LSQR_Field.tif # output fieldFileName
7 # radius
2 # sradius
1 # refine : LSQR
10 # precision
)
# A enrichir
SET(BasicFeatureExtraction_SRCS1
otbFeatureExtractionTests1.cxx
......@@ -2079,6 +2160,13 @@ otbScalarImageToHigherOrderTexturesFilterNew.cxx
otbScalarImageToHigherOrderTexturesFilter.cxx
)
SET(BasicFeatureExtraction_SRCS16
otbFeatureExtractionTests16.cxx
otbFineCorrelationImageFilterNew.cxx
otbFineCorrelationImageFilterTest.cxx
)
OTB_ADD_EXECUTABLE(otbFeatureExtractionTests1 "${BasicFeatureExtraction_SRCS1}" "OTBFeatureExtraction;OTBIO;OTBTesting")
OTB_ADD_EXECUTABLE(otbFeatureExtractionTests2 "${BasicFeatureExtraction_SRCS2}" "OTBFeatureExtraction;OTBIO;OTBTesting")
OTB_ADD_EXECUTABLE(otbFeatureExtractionTests3 "${BasicFeatureExtraction_SRCS3}" "OTBFeatureExtraction;OTBIO;OTBTesting")
......@@ -2094,6 +2182,7 @@ OTB_ADD_EXECUTABLE(otbFeatureExtractionTests12 "${BasicFeatureExtraction_SRCS12}
OTB_ADD_EXECUTABLE(otbFeatureExtractionTests13 "${BasicFeatureExtraction_SRCS13}" "OTBFeatureExtraction;OTBIO;OTBTesting")
OTB_ADD_EXECUTABLE(otbFeatureExtractionTests14 "${BasicFeatureExtraction_SRCS14}" "OTBFeatureExtraction;OTBIO;OTBTesting")
OTB_ADD_EXECUTABLE(otbFeatureExtractionTests15 "${BasicFeatureExtraction_SRCS15}" "OTBFeatureExtraction;OTBIO;OTBTesting")
OTB_ADD_EXECUTABLE(otbFeatureExtractionTests16 "${BasicFeatureExtraction_SRCS16}" "OTBFeatureExtraction;OTBIO;OTBTesting")
ENDIF( NOT OTB_DISABLE_CXX_TESTING AND BUILD_TESTING )
/*=========================================================================
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.
=========================================================================*/
// this file defines the otbCommonTest for the test driver
// and all it expects is that you have a function called RegisterTests
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#include "otbTestMain.h"
void RegisterTests()
{
REGISTER_TEST(otbFineCorrelationImageFilterNew);
REGISTER_TEST(otbFineCorrelationImageFilterTest);
}
/*=========================================================================
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.
=========================================================================*/
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#include "otbImage.h"
#include "itkFixedArray.h"
#include "otbFineCorrelationImageFilter.h"
int otbFineCorrelationImageFilterNew( int argc, char * argv[] )
{
typedef double PixelType;
typedef itk::FixedArray<PixelType> DeformationValueType;
typedef otb::Image< PixelType> ImageType;
typedef otb::Image<DeformationValueType,2> FieldImageType;
typedef otb::FineCorrelationImageFilter<ImageType,ImageType,FieldImageType> CorrelationFilterType;
CorrelationFilterType::Pointer correlation = CorrelationFilterType::New();
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.
=========================================================================*/
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#include "itkFixedArray.h"
#include "otbImage.h"
#include "otbImageFileReader.h"
#include "otbImageFileWriter.h"
#include "otbFineCorrelationImageFilter.h"
#include "otbStandardFilterWatcher.h"
#include "itkTimeProbe.h"
int otbFineCorrelationImageFilterTest( int argc, char * argv[] )
{
if(argc!=9)
{
std::cerr<<"Usage: "<<argv[0]<<" fixed_fname moving_fname output_correl output_field radius search_radius refineMethod subpixPrecision"<<std::endl;
return EXIT_FAILURE;
}
const char * fixedFileName = argv[1];
const char * movingFileName = argv[2];
const char * correlFileName = argv[3];
const char * fieldFileName = argv[4];
const unsigned int radius = atoi(argv[5]);
const unsigned int sradius = atoi(argv[6]);
const int refine = atoi(argv[7]);
const unsigned int precision = atoi(argv[8]);
typedef double PixelType;
const unsigned int Dimension = 2;
typedef itk::FixedArray<PixelType,Dimension> DeformationValueType;
typedef otb::Image< PixelType, Dimension > ImageType;
typedef otb::Image<DeformationValueType,Dimension> FieldImageType;
typedef otb::ImageFileReader< ImageType > ReaderType;
typedef otb::ImageFileWriter< ImageType > CorrelWriterType;
typedef otb::ImageFileWriter< FieldImageType> FieldWriterType;
typedef otb::FineCorrelationImageFilter<ImageType,ImageType,FieldImageType> CorrelationFilterType;
ReaderType::Pointer freader = ReaderType::New();
freader->SetFileName(fixedFileName);
//freader->Update();
ReaderType::Pointer mreader = ReaderType::New();
mreader->SetFileName(movingFileName);
//mreader->Update();
CorrelationFilterType::Pointer correlation = CorrelationFilterType::New();
correlation->SetFixedInput(freader->GetOutput());
correlation->SetMovingInput(mreader->GetOutput());
correlation->SetRadius(radius);
correlation->SetSearchRadius(sradius);
if(refine == 0)
{
correlation->SetRefinementModeToCoarse();
}
else if(refine == 1)
{
correlation->SetRefinementModeToLSQRQuadFit();
}
else
{
correlation->SetRefinementModeToSubPixel();
}
correlation->SetSubPixelPrecision(precision);
correlation->SetNumberOfThreads(1);
// otb::StandardFilterWatcher watcher(correlation,"Correlation");
itk::TimeProbe chrono;
chrono.Start();
correlation->Update();
chrono.Stop();
std::cout<<chrono.GetMeanTime()<<"\t";
CorrelWriterType::Pointer correlWriter = CorrelWriterType::New();
correlWriter->SetFileName(correlFileName);
correlWriter->SetInput(correlation->GetOutput());
correlWriter->Update();
FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
fieldWriter->SetFileName(fieldFileName);
fieldWriter->SetInput(correlation->GetOutputDeformationField());
fieldWriter->Update();
return EXIT_SUCCESS;
}
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