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

Intégration du code d'interpolation par splines de points éparts, contribution de Insight Journal.

parent cd14374f
IF( NOT OTB_DISABLE_CXX_TESTING )
# ------- ossim integration test ----------
SET(BASELINE ${OTB_DATA_ROOT}/Baseline/OTB/Images)
SET(BASELINE_FILES ${OTB_DATA_ROOT}/Baseline/OTB/Files)
SET(INPUTDATA ${OTB_DATA_ROOT}/Input)
#Images de teledetection (grosses images )
SET(IMAGEDATA ${OTB_DATA_ROOT}/LargeInput )
SET(TEMP ${OTBTesting_BINARY_DIR}/Temporary)
#Tolerance sur diff pixel image
SET(TOL 0.0)
SET(UTILITIES_TESTS ${CXX_TEST_PATH}/otbUtilitiesTests)
ADD_TEST(utOssimIntegrationTest ${UTILITIES_TESTS}
ossimIntegrationTest
)
ADD_TEST(utIjBSplineScatteredDataPointSetToImageFilterTest ${UTILITIES_TESTS}
--compare-image ${TOL}
${BASELINE}/ijBSplineScatteredDataPointSetToImageFilterTestOutput.hdr
${TEMP}/ijBSplineScatteredDataPointSetToImageFilterTestOutput.hdr
ijBSplineScatteredDataPointSetToImageFilterTest
${INPUTDATA}/brain.png
${TEMP}/ijBSplineScatteredDataPointSetToImageFilterTestOutput.hdr
)
ADD_TEST(utTuOSSIMIntegrationTest ${CXX_TEST_PATH}/utOSSIMIntegrationTest)
# ------- Fichiers sources CXX -----------------------------------
SET(UtilitiesTests_SRCS
ossimIntegrationTest.cxx
ijBSplineScatteredDataPointSetToImageFilterTest.cxx
)
INCLUDE_DIRECTORIES("${OTBTesting_BINARY_DIR}")
ADD_EXECUTABLE(utOSSIMIntegrationTest ossimIntegrationTest.cxx)
TARGET_LINK_LIBRARIES(utOSSIMIntegrationTest gdal ossim)
ADD_EXECUTABLE(otbUtilitiesTests otbUtilitiesTests.cxx ${UtilitiesTests_SRCS})
TARGET_LINK_LIBRARIES(otbUtilitiesTests OTBIO OTBCommon gdal ITKIO ITKAlgorithms ITKStatistics ITKCommon ossim)
ENDIF( NOT OTB_DISABLE_CXX_TESTING )
#include "otbImage.h"
#include "otbImageFileReader.h"
#include "otbImageFileWriter.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkPointSet.h"
#include "ijBSplineScatteredDataPointSetToImageFilter.h"
/**
* In this test, we approximate a 2-D scalar field.
* The scattered data is derived from a segmented
* image. We write the output to an image for
* comparison.
*/
int ijBSplineScatteredDataPointSetToImageFilterTest( int argc, char *argv[] )
{
const unsigned int ParametricDimension = 2;
const unsigned int DataDimension = 1;
typedef int PixelType;
typedef otb::Image<PixelType, ParametricDimension> InputImageType;
typedef float RealType;
typedef itk::Vector<RealType, DataDimension> VectorType;
typedef otb::Image<VectorType, ParametricDimension> VectorImageType;
typedef itk::PointSet
<VectorImageType::PixelType, ParametricDimension> PointSetType;
PointSetType::Pointer pointSet = PointSetType::New();
typedef otb::ImageFileReader<InputImageType> ReaderType;
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName( argv[1] );
reader->Update();
itk::ImageRegionIteratorWithIndex<InputImageType>
It( reader->GetOutput(), reader->GetOutput()->GetLargestPossibleRegion() );
// Iterate through the input image which consists of multivalued
// foreground pixels (=nonzero) and background values (=zero).
// The foreground pixels comprise the input point set.
for ( It.GoToBegin(); !It.IsAtEnd(); ++It )
{
if ( It.Get() != itk::NumericTraits<PixelType>::Zero )
{
// We extract both the 2-D location of the point
// and the pixel value of that point.
PointSetType::PointType point;
reader->GetOutput()->TransformIndexToPhysicalPoint( It.GetIndex(), point );
unsigned long i = pointSet->GetNumberOfPoints();
pointSet->SetPoint( i, point );
PointSetType::PixelType V( DataDimension );
V[0] = static_cast<RealType>( It.Get() );
pointSet->SetPointData( i, V );
}
}
// Instantiate the B-spline filter and set the desired parameters.
typedef ij::BSplineScatteredDataPointSetToImageFilter
<PointSetType, VectorImageType> FilterType;
FilterType::Pointer filter = FilterType::New();
filter->SetSplineOrder( 3 );
FilterType::ArrayType ncps;
ncps.Fill( 4 );
filter->SetNumberOfControlPoints( ncps );
filter->SetNumberOfLevels( 10 );
// Define the parametric domain.
filter->SetOrigin( reader->GetOutput()->GetOrigin() );
filter->SetSpacing( reader->GetOutput()->GetSpacing() );
filter->SetSize( reader->GetOutput()->GetLargestPossibleRegion().GetSize() );
filter->SetInput( pointSet );
try
{
filter->Update();
}
catch (...)
{
std::cerr << "Test 1: itkBSplineScatteredDataImageFilter exception thrown"
<< std::endl;
return EXIT_FAILURE;
}
// Write the output to an image.
typedef otb::Image<RealType, ParametricDimension> RealImageType;
RealImageType::Pointer image = RealImageType::New();
image->SetRegions( reader->GetOutput()->GetLargestPossibleRegion() );
image->Allocate();
itk::ImageRegionIteratorWithIndex<RealImageType>
Itt( image, image->GetLargestPossibleRegion() );
for ( Itt.GoToBegin(); !Itt.IsAtEnd(); ++Itt )
{
Itt.Set( filter->GetOutput()->GetPixel( Itt.GetIndex() )[0] );
}
typedef otb::ImageFileWriter<RealImageType> WriterType;
WriterType::Pointer writer = WriterType::New();
writer->SetInput( image );
writer->SetFileName( argv[2] );
writer->Update();
return EXIT_SUCCESS;
};
......@@ -50,7 +50,7 @@ void printOutputTypes();
ossimProjection* newUtmView(const ossimGpt& centerGround,
const ossimDpt& metersPerPixel);
int main(int argc, char* argv[])
int ossimIntegrationTest(int argc, char* argv[])
{
ossimInit::instance()->initialize(argc, argv);
......@@ -163,7 +163,7 @@ ossimProjection* newUtmView(const ossimGpt& centerGround,
ossimUtmProjection* utm = new ossimUtmProjection;
// we will make it a square pixel in meters
double averageGsd = (metersPerPixel.x + metersPerPixel.y)*.5;
//double averageGsd = (metersPerPixel.x + metersPerPixel.y)*.5;
utm->setZone(centerGround);
utm->setMetersPerPixel(ossimDpt(metersPerPixel));
......
/*=========================================================================
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 <iostream>
#include "otbTestMain.h"
void RegisterTests()
{
REGISTER_TEST(ossimIntegrationTest);
REGISTER_TEST(ijBSplineScatteredDataPointSetToImageFilterTest);
}
......@@ -4,7 +4,7 @@ IF(NOT OTB_USE_EXTERNAL_ITK)
SUBDIRS( ITK )
ENDIF(NOT OTB_USE_EXTERNAL_ITK)
SUBDIRS(BGL otbsvm dxflib)
SUBDIRS(BGL otbsvm dxflib InsightJournal)
OPTION(OTB_OSSIM_COMPILATION "EXPERIMENTAL: Compile OSSIM library embedded within OTB" OFF)
MARK_AS_ADVANCED(OTB_OSSIM_COMPILATION)
......
# Sources of non-templated classes.
FILE(GLOB InsightJournal_SRCS "*.cxx" )
ADD_LIBRARY(InsightJournal ${InsightJournal_SRCS})
TARGET_LINK_LIBRARIES (InsightJournal ITKCommon)
INSTALL(TARGETS InsightJournal
RUNTIME DESTINATION ${OTB_INSTALL_BIN_DIR} COMPONENT RuntimeLibraries
LIBRARY DESTINATION ${OTB_INSTALL_LIB_DIR} COMPONENT RuntimeLibraries
ARCHIVE DESTINATION ${OTB_INSTALL_LIB_DIR} COMPONENT Development)
FILE(GLOB __files1 "${CMAKE_CURRENT_SOURCE_DIR}/*.h")
FILE(GLOB __files2 "${CMAKE_CURRENT_SOURCE_DIR}/*.txx")
INSTALL(FILES ${__files1} ${__files2}
DESTINATION ${OTB_INSTALL_INCLUDE_DIR}/Utilities/InsightJournal
COMPONENT Development)
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: ijBSplineKernelFunction.h,v $
Language: C++
Date: $Date: $
Version: $Revision: $
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm 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 _ijBSplineKernelFunction_h
#define _ijBSplineKernelFunction_h
#include "itkKernelFunction.h"
#include "vnl/vnl_math.h"
#include "vnl/vnl_real_polynomial.h"
#include "vnl/vnl_matrix.h"
namespace ij
{
/** \class BSplineKernelFunction
* \brief BSpline kernel used for density estimation and nonparameteric
* regression.
*
* This class enscapsulates BSpline kernel for
* density estimation or nonparameteric regression.
* See documentation for KernelFunction for more details.
*
* This class is templated over the spline order to cohere with
* the previous incarnation of this class. One can change the
* order during an instantiation's existence. Note that
* other authors have defined the B-spline order as being the
* degree of spline + 1. In the ITK context (e.g. in this
* class), the spline order is equivalent to the degree of
* the spline.
*
* \sa KernelFunction
*
* \ingroup Functions
*/
template <unsigned int VSplineOrder = 3>
class ITK_EXPORT BSplineKernelFunction
: public itk::KernelFunction
{
public:
/** Standard class typedefs. */
typedef BSplineKernelFunction Self;
typedef itk::KernelFunction Superclass;
typedef itk::SmartPointer<Self> Pointer;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro( BSplineKernelFunction, KernelFunction );
typedef double RealType;
typedef vnl_vector<RealType> VectorType;
typedef vnl_real_polynomial PolynomialType;
typedef vnl_matrix<RealType> MatrixType;
/** Get/Sets the Spline Order */
void SetSplineOrder( unsigned int );
itkGetMacro( SplineOrder, unsigned int );
/** Evaluate the function. */
inline RealType Evaluate( const RealType & u ) const
{
RealType absValue = vnl_math_abs( u );
int which;
if ( this->m_SplineOrder % 2 == 0 )
{
which = static_cast< int>( absValue+0.5 );
}
else
{
which = static_cast< int>( absValue );
}
if ( which < this->m_BSplineShapeFunctions.rows() )
{
return PolynomialType(
this->m_BSplineShapeFunctions.get_row( which ) ).evaluate( absValue );
}
else
{
return itk::NumericTraits<RealType>::Zero;
}
}
/** Evaluate the derivative. */
inline RealType EvaluateDerivative( const double & u ) const
{
RealType absValue = vnl_math_abs( u );
int which;
if ( this->m_SplineOrder % 2 == 0 )
{
which = static_cast<unsigned int>( absValue+0.5 );
}
else
{
which = static_cast<unsigned int>( absValue );
}
if ( which < this->m_BSplineShapeFunctions.rows() )
{
RealType der = PolynomialType(
this->m_BSplineShapeFunctions.get_row( which ) ).devaluate( absValue );
if ( u < itk::NumericTraits<RealType>::Zero )
{
return -der;
}
else
{
return der;
}
}
else
{
return itk::NumericTraits<RealType>::Zero;
}
}
/**
* For a specific order, return the (this->m_SplineOrder+1) pieces of
* the single basis function centered at zero.
*/
MatrixType GetShapeFunctions();
/**
* For a specific order, generate and return the (this->m_SplineOrder+1)
* pieces of the different basis functions in the [0, 1] interval.
*/
MatrixType GetShapeFunctionsInZeroToOneInterval();
protected:
BSplineKernelFunction();
~BSplineKernelFunction();
void PrintSelf( std::ostream& os, itk::Indent indent ) const;
private:
BSplineKernelFunction( const Self& ); //purposely not implemented
void operator=( const Self& ); //purposely not implemented
/**
* For a specific order, generate the (this->m_SplineOrder+1) pieces of
* the single basis function centered at zero.
*/
void GenerateBSplineShapeFunctions( unsigned int );
/**
* Use the CoxDeBoor recursion relation to generate the piecewise
* polynomials which compose the basis function.
* See, for example, L. Piegl, L. Tiller, "The NURBS Book,"
* Springer 1997, p. 50.
*/
PolynomialType CoxDeBoor( unsigned short, VectorType, unsigned int, unsigned int );
MatrixType m_BSplineShapeFunctions;
unsigned int m_SplineOrder;
};
} // end namespace itk
#ifndef OTB_MANUAL_INSTANTIATION
#include "ijBSplineKernelFunction.txx"
#endif
#endif
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: itkBSplineKernelFunction.txx,v $
Language: C++
Date: $Date: $
Version: $Revision: $
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm 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 _ijBSplineKernelFunction_txx
#define _ijBSplineKernelFunction_txx
#include "ijBSplineKernelFunction.h"
namespace ij
{
template <unsigned int VSplineOrder>
BSplineKernelFunction<VSplineOrder>
::BSplineKernelFunction()
{
this->m_SplineOrder = VSplineOrder;
this->GenerateBSplineShapeFunctions( this->m_SplineOrder+1 );
}
template <unsigned int VSplineOrder>
BSplineKernelFunction<VSplineOrder>
::~BSplineKernelFunction()
{
}
template <unsigned int VSplineOrder>
void
BSplineKernelFunction<VSplineOrder>
::SetSplineOrder( unsigned int order )
{
if ( order != this->m_SplineOrder )
{
this->m_SplineOrder = order;
this->GenerateBSplineShapeFunctions( this->m_SplineOrder+1 );
this->Modified();
}
}
template <unsigned int VSplineOrder>
void
BSplineKernelFunction<VSplineOrder>
::GenerateBSplineShapeFunctions( unsigned int order )
{
unsigned int NumberOfPieces = static_cast<unsigned int>( 0.5*( order+1 ) );
this->m_BSplineShapeFunctions.set_size( NumberOfPieces, order );
VectorType knots( order+1 );
for ( unsigned int i = 0; i < knots.size(); i++)
{
knots[i] = -0.5*static_cast<RealType>( order ) + static_cast<RealType>( i );
}
for ( unsigned int i = 0; i < NumberOfPieces; i++ )
{
PolynomialType poly = this->CoxDeBoor(order, knots,
0, static_cast<unsigned int>( 0.5*( order ) ) + i );
this->m_BSplineShapeFunctions.set_row( i, poly.coefficients() );
}
}
template <unsigned int VSplineOrder>
typename BSplineKernelFunction<VSplineOrder>::PolynomialType
BSplineKernelFunction<VSplineOrder>
::CoxDeBoor( unsigned short order, VectorType knots,
unsigned int whichBasisFunction, unsigned int whichPiece )
{
VectorType tmp(2);
PolynomialType poly1(0.0), poly2(0.0);
RealType den;
unsigned short p = order-1;
unsigned short i = whichBasisFunction;
if ( p == 0 && whichBasisFunction == whichPiece )
{
PolynomialType poly(1.0);
return poly;
}
// Term 1
den = knots(i+p)-knots(i);
if ( den == itk::NumericTraits<RealType>::Zero )
{
PolynomialType poly(0.0);
poly1 = poly;
}
else
{
tmp(0) = 1.0;
tmp(1) = -knots(i);
tmp /= den;
poly1 = PolynomialType(tmp) * this->CoxDeBoor( order-1, knots, i, whichPiece );
}
// Term 2
den = knots(i+p+1)-knots(i+1);
if ( den == itk::NumericTraits<RealType>::Zero )
{
PolynomialType poly(0.0);
poly2 = poly;
}
else
{
tmp(0) = -1.0;
tmp(1) = knots(i+p+1);
tmp /= den;
poly2 = PolynomialType(tmp) * this->CoxDeBoor( order-1, knots, i+1, whichPiece );
}
return ( poly1 + poly2 );
}
template <unsigned int VSplineOrder>
typename BSplineKernelFunction<VSplineOrder>::MatrixType
BSplineKernelFunction<VSplineOrder>
::GetShapeFunctionsInZeroToOneInterval()
{
int order = this->m_SplineOrder+1;
unsigned int NumberOfPieces = static_cast<unsigned int>( order );
MatrixType ShapeFunctions( NumberOfPieces, order );
VectorType knots( 2*order );
for ( unsigned int i = 0; i < knots.size(); i++ )
{
knots[i] = -static_cast<RealType>( this->m_SplineOrder )
+ static_cast<RealType>( i );
}
for (unsigned int i = 0; i < NumberOfPieces; i++ )
{
PolynomialType poly = this->CoxDeBoor( order, knots, i, order-1 );
ShapeFunctions.set_row( i, poly.coefficients() );
}
return ShapeFunctions;
}
template <unsigned int VSplineOrder>
void
BSplineKernelFunction<VSplineOrder>
::PrintSelf( std::ostream& os, itk::Indent indent ) const
{
Superclass::PrintSelf( os, indent );
os << indent << "Spline Order: " << this->m_SplineOrder << std::endl;
os << indent << "Piecewise Polynomial Pieces: " << std::endl;
RealType a, b;
for ( unsigned int i = 0; i < this->m_BSplineShapeFunctions.rows(); i++ )
{
os << indent << indent;
PolynomialType( this->m_BSplineShapeFunctions.get_row( i ) ).print( os );
if ( i == 0 )
{
a = 0.0;
if ( this->m_SplineOrder % 2 == 0 )
{
b = 0.5;
}
else
{
b = 1.0;
}
}
else
{
a = b;
b += 1.0;
}
os << ", X \\in [" << a << ", " << b << "]" << std::endl;
}
}