[Insight-developers] itkImageMapper
Lydia Ng
lng@insightful.com
Mon, 3 Dec 2001 12:45:46 -0800
This is a multi-part message in MIME format.
------=_NextPart_000_0013_01C17BF8.7AD834A0
Content-Type: text/plain;
charset="us-ascii"
Content-Transfer-Encoding: 7bit
Hi Josh,
Attached is a 2D test for MI/Affine/GradientDescent
( It basically went it and strip the existing 3D test
and made it 2D ).
It was able to recover a 20 percent dilation -
so I presume all the components are working in 2D ...
Lydia
> -----Original Message-----
> From: insight-developers-admin@public.kitware.com
> [mailto:insight-developers-admin@public.kitware.com]On Behalf
> Of Joshua
> Cates
> Sent: Monday, December 03, 2001 10:21 AM
> To: Insight-Developers
> Subject: [Insight-developers] itkImageMapper
>
>
> Hi,
>
> Can someone verify that itkImageMapper is working properly with 2d
> transforms on 2d images?
>
> What inputs does itkImageMapper require besides a transform
> and an image
> pointer?
>
> Thanks,
>
>
> Josh.
>
>
> ______________________________
> Josh Cates
> School of Computer Science
> University of Utah
> Email: cates@cs.utah.edu
> Phone: (801) 587-7697
> URL: www.cs.utk.edu/~cates
>
>
> _______________________________________________
> Insight-developers mailing list
> Insight-developers@public.kitware.com
> http://public.kitware.com/mailman/listinfo/insight-developers
------=_NextPart_000_0013_01C17BF8.7AD834A0
Content-Type: text/plain;
name="CMakeLists.txt"
Content-Transfer-Encoding: quoted-printable
Content-Disposition: attachment;
filename="CMakeLists.txt"
PROJECT(TestOfAffineMI2D)
#=20
# Find ITK
#
FIND_PATH( ITK_BINARY_DIR itkConfigure.h )
# Load in the values from ITK if found
IF ( ITK_BINARY_DIR )
LOAD_CACHE(${ITK_BINARY_DIR})
INCLUDE (${ITK_SOURCE_DIR}/itkCMakeOptions.cmake)
ENDIF (ITK_BINARY_DIR )
INCLUDE_DIRECTORIES(
${ITK_SOURCE_DIR}/Code/Common
${ITK_SOURCE_DIR}/Code/BasicFilters
${ITK_SOURCE_DIR}/Code/Algorithms
)
LINK_DIRECTORIES(
${ITK_BINARY_DIR}/Code/Common
${ITK_BINARY_DIR}/Code/BasicFilters
${ITK_BINARY_DIR}/Code/Algorithms
)
LINK_LIBRARIES (
VXLNumerics
ITKCommon
ITKBasicFilters=20
)
ADD_EXECUTABLE( TestOfAffineMI2D =
itkImageToImageAffineMutualInformationGradientDescentRegistrationTest )
------=_NextPart_000_0013_01C17BF8.7AD834A0
Content-Type: application/octet-stream;
name="itkImageToImageAffineMutualInformationGradientDescentRegistrationTest.cxx"
Content-Transfer-Encoding: quoted-printable
Content-Disposition: attachment;
filename="itkImageToImageAffineMutualInformationGradientDescentRegistrationTest.cxx"
/*=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D
Program: Insight Segmentation & Registration Toolkit (ITK)
Module: $RCSfile: =
itkImageToImageAffineMutualInformationGradientDescentRegistrationTest.cxx=
,v $
Language: C++
Date: $Date: 2001/11/30 04:19:44 $
Version: $Revision: 1.24 $
Copyright (c) 2001 Insight Consortium
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are =
met:
* Redistributions of source code must retain the above copyright =
notice,
this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright =
notice,
this list of conditions and the following disclaimer in the =
documentation
and/or other materials provided with the distribution.
* The name of the Insight Consortium, nor the names of any consortium =
members,
nor of any contributors, may be used to endorse or promote products =
derived
from this software without specific prior written permission.
* Modified source versions must be plainly marked as such, and must =
not be
misrepresented as being the original software.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS ``AS =
IS''
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, =
THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR =
PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE =
FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS =
OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) =
HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT =
LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF =
THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D*/
#include "itkImage.h"
#include "itkImageRegionIterator.h"
#include =
"itkImageToImageAffineMutualInformationGradientDescentRegistration.h"
#include "vnl/vnl_math.h"
#include "itkExceptionObject.h"
#include "itkCommandIterationUpdate.h"
#include <iostream>
/**
* This function defines the test image pattern.
* The pattern is a 2D gaussian in the middle
* and some directional pattern on the outside.
*/
double F( double x, double y )
{
const double s =3D 50;
double value =3D 200.0 * exp( - ( x*x + y*y )/(s*s) );
x -=3D 8; y +=3D 3;=20
double r =3D vnl_math_sqrt( x*x + y*y );
if( r > 35 )
{
value =3D 2 * ( vnl_math_abs( x ) +
0.8 * vnl_math_abs( y ) );
}
if( r < 4 )
{
value =3D 400;
}
return value;
}
int main()
{
//------------------------------------------------------------
// Create two simple images
// Translate and dilate one of the image
//------------------------------------------------------------
// Allocate Images
typedef float PixelType;
typedef itk::Image<PixelType,2> ReferenceType;
typedef itk::Image<PixelType,2> TargetType;
enum { ImageDimension =3D ReferenceType::ImageDimension };
ReferenceType::SizeType size =3D {{100,100}};
ReferenceType::IndexType index =3D {{0,0}};
ReferenceType::RegionType region;
region.SetSize( size );
region.SetIndex( index );
ReferenceType::Pointer imgReference =3D ReferenceType::New();
imgReference->SetLargestPossibleRegion( region );
imgReference->SetBufferedRegion( region );
imgReference->SetRequestedRegion( region );
imgReference->Allocate();
TargetType::Pointer imgTarget =3D TargetType::New();
imgTarget->SetLargestPossibleRegion( region );
imgTarget->SetBufferedRegion( region );
imgTarget->SetRequestedRegion( region );
imgTarget->Allocate();
// Fill images with a 2D gaussian with some directional pattern
// in the background
typedef itk::ImageRegionIterator<ReferenceType>
ReferenceIteratorType;
typedef itk::ImageRegionIterator<TargetType>
TargetIteratorType;
itk::Point<double,2> center;
center[0] =3D (double)region.GetSize()[0]/2.0;
center[1] =3D (double)region.GetSize()[1]/2.0;
itk::Point<double,2> p;
itk::Vector<double,2> d;
// Set the displacement and scale
itk::Vector<double,2> displacement;
displacement[0] =3D 7;
displacement[1] =3D 3;
// Twenty percent dilation
double scale[3] =3D { 0.80, 1.0 };
ReferenceIteratorType ri(imgReference,region);
TargetIteratorType ti(imgTarget,region);
while(!ri.IsAtEnd())
{
p[0] =3D ri.GetIndex()[0];
p[1] =3D ri.GetIndex()[1];
d =3D p-center;
const double x =3D d[0] * scale[0] + displacement[0];
const double y =3D d[1] * scale[1] + displacement[1];
ri.Set( (PixelType) F(x,y) );
++ri;
}
while(!ti.IsAtEnd())
{
p[0] =3D ti.GetIndex()[0];
p[1] =3D ti.GetIndex()[1];
d =3D p-center;
const double x =3D d[0];
const double y =3D d[1];
ti.Set( (PixelType) F(x,y) );
++ti;
}
// set image origin to be center of the image
double transCenter[2];
for( unsigned int j =3D 0; j < 2; j++ )
{
transCenter[j] =3D -0.5 * double(size[j]);
}
imgReference->SetOrigin( transCenter );
imgTarget->SetOrigin( transCenter );
//-----------------------------------------------------------
// Set up a the registrator
//-----------------------------------------------------------
typedef =
itk::ImageToImageAffineMutualInformationGradientDescentRegistration<
ReferenceType,TargetType> RegistrationType;
RegistrationType::Pointer registrationMethod =3D =
RegistrationType::New();
// connect the images
registrationMethod->SetReference(imgReference);
registrationMethod->SetTarget(imgTarget);
// set translation scale
typedef RegistrationType::OptimizerType OptimizerType;
typedef OptimizerType::TransformType::ParametersType ScaleType;
ScaleType scales;
scales.Fill( 1.0 );
for( unsigned j =3D 4; j < 6; j++ )
{
scales[j] =3D 0.0001;
}
registrationMethod->GetOptimizer()->GetTransform()->SetScale( scales =
);
// set metric related parameters
registrationMethod->GetMetric()->SetTargetStandardDeviation( 5.0 );
registrationMethod->GetMetric()->SetReferenceStandardDeviation( 5.0 );
registrationMethod->GetMetric()->SetNumberOfSpatialSamples( 50 );
// Define the type for the observer command to monitor progress
typedef itk::CommandIterationUpdate< RegistrationType::OptimizerType =
>
=
CommandIterationType;
CommandIterationType::Pointer iterationCommand =3D =
CommandIterationType::New();
iterationCommand->SetOptimizer( registrationMethod->GetOptimizer() );
registrationMethod->GetOptimizer()->AddObserver( =
itk::IterationEvent(),
iterationCommand );
// do the registration
// reduce learning rate as we go
unsigned int iter[3] =3D {300,300,350};
double rates[3] =3D {1e-3, 5e-4, 1e-4};
for( unsigned int i =3D 0; i < 3; i++ )
{
registrationMethod->SetNumberOfIterations( iter[i] );
registrationMethod->SetLearningRate( rates[i] );
try
{
registrationMethod->StartRegistration();
}
catch(itk::ExceptionObject& err)
{
// caught an exception object
std::cout << "Caught an ExceptionObject" << std::endl;
std::cout << err << std::endl;
std::cout << "Test failed." << std::endl;
return EXIT_FAILURE;
}
catch(...)
{
// caught some other error
std::cout << "Caught unknown exception" << std::endl;
std::cout << "Test failed. " << std::endl;
return EXIT_FAILURE;
}
}
// get the results
RegistrationType::ParametersType solution =3D
registrationMethod->GetParameters();
std::cout << "Solution is: " << solution << std::endl;
//
// check results to see if it is within range
//
bool pass =3D true;
double trueParameters[12] =3D { 1, 0, 0, 1, 0, 0 };
trueParameters[0] =3D 1/scale[0];
trueParameters[3] =3D 1/scale[1];
trueParameters[4] =3D - displacement[0]/scale[0];
trueParameters[5] =3D - displacement[1]/scale[1];
std::cout << "True solution is: ";
for ( unsigned int j =3D 0; j < 6; j++)
std::cout << trueParameters[j] << " ";
std::cout << std::endl;
for( unsigned int j =3D 0; j < 4; j++ )
{
if( vnl_math_abs( solution[j] - trueParameters[j] ) > 0.025 )
{
pass =3D false;
}
}
for( unsigned int j =3D 4; j < 6; j++ )
{
if( vnl_math_abs( solution[j] - trueParameters[j] ) > 1.0 )
{
pass =3D false;
}
}
if( !pass )
{
std::cout << "Test failed." << std::endl;
return EXIT_FAILURE;
}
// check for parzen window exception
double oldValue =3D =
registrationMethod->GetMetric()->GetReferenceStandardDeviation();
registrationMethod->GetMetric()->SetReferenceStandardDeviation( 0.005 =
);
try
{
pass =3D false;
registrationMethod->StartRegistration();
}
catch(itk::ExceptionObject& err)
{
std::cout << "Caught expected ExceptionObject" << std::endl;
std::cout << err << std::endl;
pass =3D true;
}
if( !pass )
{
std::cout << "Should have caught an exception" << std::endl;
std::cout << "Test failed." << std::endl;
return EXIT_FAILURE;
}
// check for mapped out of image error
registrationMethod->GetMetric()->SetReferenceStandardDeviation( =
oldValue );
solution[5] =3D 1000;
registrationMethod->SetParameters(solution);
try
{
pass =3D false;
registrationMethod->StartRegistration();
}
catch(itk::ExceptionObject& err)
{
std::cout << "Caught expected ExceptionObject" << std::endl;
std::cout << err << std::endl;
pass =3D true;
}
if( !pass )
{
std::cout << "Should have caught an exception" << std::endl;
std::cout << "Test failed." << std::endl;
return EXIT_FAILURE;
}
std::cout << "Test passed." << std::endl;
return EXIT_SUCCESS;
}
------=_NextPart_000_0013_01C17BF8.7AD834A0
Content-Type: application/octet-stream;
name="itkCommandIterationUpdate.h"
Content-Transfer-Encoding: quoted-printable
Content-Disposition: attachment;
filename="itkCommandIterationUpdate.h"
/*=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D
Program: Insight Segmentation & Registration Toolkit (ITK)
Module: $RCSfile: itkCommandIterationUpdate.h,v $
Language: C++
Date: $Date: 2001/11/12 13:39:40 $
Version: $Revision: 1.5 $
Copyright (c) 2001 Insight Consortium
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are =
met:
* Redistributions of source code must retain the above copyright =
notice,
this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright =
notice,
this list of conditions and the following disclaimer in the =
documentation
and/or other materials provided with the distribution.
* The name of the Insight Consortium, nor the names of any consortium =
members,
nor of any contributors, may be used to endorse or promote products =
derived
from this software without specific prior written permission.
* Modified source versions must be plainly marked as such, and must =
not be
misrepresented as being the original software.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS ``AS =
IS''
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, =
THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR =
PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE =
FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS =
OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) =
HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT =
LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF =
THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D*/
#ifndef __itkCommandIterationUpdate_h
#define __itkCommandIterationUpdate_h
#include <itkCommand.h>
#include "itkWeakPointer.h"
namespace itk {
/**
* Implementation of the Command Pattern to be invoked every iteration
*/
template < class TOptimizer >
class ITK_EXPORT CommandIterationUpdate : public Command=20
{
public:
/**
* Standard "Self" typedef.
*/
typedef CommandIterationUpdate Self;
/**
* Standard "Superclass" typedef.
*/
typedef itk::Command Superclass;
/**
* Smart pointer typedef support.
*/
typedef itk::SmartPointer<Self> Pointer;
/**=20
* ConstSmart pointer typedef support.
*/
typedef itk::SmartPointer<const Self> ConstPointer;
/**
* Execute method will print data at each iteration
*/
void Execute(itk::Object *caller, const itk::EventObject & event)
{
Execute( (const itk::Object *)caller, event);
}
void Execute(const itk::Object *caller, const itk::EventObject & =
event)
{
if( typeid( event ) =3D=3D typeid( itk::StartEvent ) )
{
std::cout << std::endl << "Position Value";
std::cout << std::endl << std::endl;
} =20
else if( typeid( event ) =3D=3D typeid( itk::IterationEvent ) )
{
std::cout << m_Optimizer->GetCurrentIteration() << " =3D ";
std::cout << m_Optimizer->GetValue() << " : ";
std::cout << m_Optimizer->GetCurrentPosition() << std::endl;
}
else if( typeid( event ) =3D=3D typeid( itk::EndEvent ) )
{
std::cout << std::endl << std::endl;
std::cout << "After " << m_Optimizer->GetCurrentIteration();
std::cout << " iterations " << std::endl;
std::cout << "Solution is =3D " << =
m_Optimizer->GetCurrentPosition();
std::cout << std::endl;
std::cout << "With value =3D " << m_Optimizer->GetValue();
std::cout << std::endl;
std::cout << "Stop condition =3D " << =
m_Optimizer->GetStopCondition();
std::cout << std::endl;
}
}
/**=20
* Run-time type information (and related methods).
*/
itkTypeMacro( CommandIterationUpdate, ::itk::Command );
/**
* Method for creation through the object factory.
*/
itkNewMacro( Self );
=20
/**
* Type defining the optimizer
*/
typedef TOptimizer OptimizerType;
/**
* Set Optimizer
*/
void SetOptimizer( OptimizerType * optimizer )
{ m_Optimizer =3D optimizer; }
protected:
/**
* Constructor
*/
CommandIterationUpdate() {};
private:
/**
* WeakPointer to the Optimizer
*/
WeakPointer<OptimizerType> m_Optimizer;
=20
=20
};
} // end namespace itk
#endif
------=_NextPart_000_0013_01C17BF8.7AD834A0--