[Insight-users] Query on Resampling
Lydia Ng
lng@insightful.com
Tue, 12 Nov 2002 12:44:47 -0800
This is a multi-part message in MIME format.
------_=_NextPart_001_01C28A8C.56D95676
Content-Type: text/plain;
charset="iso-8859-1"
Content-Transfer-Encoding: quoted-printable
I had a look through the snippet and I can't see anything that is =
obviously wrong.
=20
The MultiResMIRegistration App also uses the ResampleImageFilter
(with affine transform and linear interpolation) and I just verified =
that
it is working ok - so I am confident that the ResampleImageFilter is
working as expected.
=20
Other suggestions:
- Are you setting spacing correctly?
- Is it an IO problem? Endian-ness?
=20
I notice that you convert your volume once for registration and once for =
the
resampling - are you doing anything different between the two?
What happens if you do the resampling with exactly the same into as
for registration (ie. without conversion)?
=20
You code contains other non-itk stuff so I can't compile and debug.
If you can reproduce the problem with ITK only code, package with
small image then we can have a better chance in isolating the problem.
=20
- Lydia=20
=20
-----Original Message-----
From: cspl [mailto:affable@hd2.dot.net.in]
Sent: Monday, November 11, 2002 10:26 PM
To: Lydia Ng
Cc: insight-users@public.kitware.com
Subject: Re: [Insight-users] Query on Resampling
Hi!
=20
Thankyou for your response.
Please find herewith the code snippet of Registration and Resample =
volume.
We set all the reigstration components. For registration and resampling =
we set the transformation componentto Affine Transform and interpolation =
component to Linear interpolation.
=20
=20
Registration:=20
=20
//Converting volume data to ITK image data
typedef BufferToImageConversion<short> ConverterType;=20
typedef ConverterType::ImageType FixedImageType;
typedef itk::StatisticsImageFilter<FixedImageType> =
StatisticsImageFilterType;
=20
//Resgitartion components
typedef itk::AffineTransform<double> TransformType;
typedef itk::RegularStepGradientDescentOptimizer OptimizerType;=20
typedef itk::LinearInterpolateImageFunction <FixedImageType, double> =
InterpolatorType;
typedef itk::ImageRegistrationMethod <FixedImageType , FixedImageType =
> RegistrationType1; =20
=20
//Get Fixed and Moving images
ConverterType converter;
FixedImageType::Pointer fixedImage =3D converter.GetImage(pFixed);
FixedImageType::Pointer movingImage =3D =
converter.GetImage(pMoving);
=20
=20
////////////////////////////////////////////////////////////////////
//Set up the statistic filter=20
////////////////////////////////////////////////////////////////////
=20
StatisticsImageFilterType::Pointer fixedImageStatisticsFilter =3D=20
StatisticsImageFilterType::New();
fixedImageStatisticsFilter->SetInput( fixedImage );
fixedImageStatisticsFilter->Update();=20
=20
StatisticsImageFilterType::Pointer movingImageStatisticsFilter =3D=20
StatisticsImageFilterType::New();
movingImageStatisticsFilter->SetInput( movingImage );
movingImageStatisticsFilter->Update();=20
=20
=20
=20
////////////////////////////////////////////////////////////////////
//Set up the metric.=20
////////////////////////////////////////////////////////////////////
=20
/*metric->SetMovingImageStandardDeviation( =
movingImageStatisticsFilter->GetSigma() * 0.4 );
metric->SetFixedImageStandardDeviation( =
fixedImageStatisticsFilter->GetSigma() * 0.4 );
metric->SetNumberOfSpatialSamples(m_spatialsamples);
metric->SetFixedImageRegion( fixedImage->GetBufferedRegion() );*/
=20
RegistrationType1::Pointer registration =3D =
RegistrationType1::New(); =20
=20
if (m_metrictype=3D=3D1) =20
{
typedef itk::MutualInformationImageToImageMetric<FixedImageType, =
FixedImageType> MetricType;
MetricType::Pointer metric =3D MetricType::New();
metric->SetMovingImageStandardDeviation( =
movingImageStatisticsFilter->GetSigma() * 0.4 );
metric->SetFixedImageStandardDeviation( =
fixedImageStatisticsFilter->GetSigma() * 0.4 );
metric->SetNumberOfSpatialSamples(m_spatialsamples);
metric->SetFixedImageRegion( fixedImage->GetBufferedRegion() );
registration->SetMetric( metric );
}
else if (m_metrictype=3D=3D2)=20
{
typedef itk::NormalizedCorrelationImageToImageMetric<FixedImageType, =
FixedImageType> MetricType;
MetricType::Pointer metric =3D MetricType::New();
metric->SetFixedImageRegion( fixedImage->GetBufferedRegion() );
registration->SetMetric( metric );
}
else =20
{
typedef itk::MeanSquaresImageToImageMetric<FixedImageType, =
FixedImageType> MetricType;
MetricType::Pointer metric =3D MetricType::New();
metric->SetFixedImageRegion( fixedImage->GetBufferedRegion() );
registration->SetMetric( metric );
}
=20
=20
Matrix *mat =3D NULL;
double result[12];
=20
TransformType::Pointer transform =3D TransformType::New();
OptimizerType::Pointer optimizer =3D OptimizerType::New();
InterpolatorType::Pointer interpolator =3D =
InterpolatorType::New();
=20
=20
////////////////////////////////////////////////////////////////////
// Set up initial transform parameters
//////////////////////////////////////////////////////////////////// =
=20
RegistrationType1::ParametersType initialParameters(=20
transform->GetNumberOfParameters() );
=20
//Set up transform parameters
// transform->SetIdentity();
=20
TransformType::OutputVectorType shift;
shift[0] =3D1;
shift[1] =3D0;
shift[2] =3D0; // translation of (10,15,20)
// transform->Translate( shift );
=20
TransformType::OutputVectorType axis;
axis[0] =3D 1.0;
axis[1] =3D 0.0;
axis[2] =3D 0.0; // axis=3D vector parallel to the Z axis
const double angle =3D 100;
// transform->Rotate3D( axis, angle,0 );
=20
// TransformType::OutputVectorType shift;
shift[0] =3D1;
shift[1] =3D 0;
shift[2] =3D 0; // translation of (10,15,20)
// transform->Translate( shift );=20
/* TransformType::OutputVectorType shr;
shr[0]=3D0.5;
shr[1]=3D0.0;
shr[2]=3D0.0;
transform->Scale(shr,0);*/
=20
=20
=20
registration->SetInitialTransformParameters( =
transform->GetParameters());
=20
=20
/* TransformType::ParametersType parameter;=20
parameter=3Dtransform->GetParameters();
CString str;
for(int i=3D0;i<12;i++)
{ =20
double parameters[12];
parameters[i]=3D parameter[i];
str.Format("%lf", parameters[i]);
AfxMessageBox(str);
=20
}
*/
=20
////////////////////////////////////////////////////////////////////
// Set up the optimizer.
//////////////////////////////////////////////////////////////////// =
=20
=20
typedef OptimizerType::ScalesType ScalesType;
ScalesType parametersScales( transform->GetNumberOfParameters() );
=20
parametersScales.Fill( 1.0 );
=20
double scale =3D 1.0 / vnl_math_sqr(m_translatescale);
=20
for (int j =3D 9; j < 12; j++ )
{
=20
parametersScales[j] =3D scale;
}
=20
optimizer->SetScales( parametersScales );
=20
optimizer->SetNumberOfIterations(m_iterations);
=20
=20
////////////////////////////////////////////////////////////////////
// Set up the registrator.
////////////////////////////////////////////////////////////////////
=20
// connect up the components
=20
registration->SetOptimizer( optimizer );
registration->SetTransform( transform );
registration->SetFixedImage(fixedImage );
registration->SetMovingImage( movingImage );
registration->SetInterpolator( interpolator );
=20
ProgressUpdate(75, " Start registration", MRIVolumeName);
try
{
registration->StartRegistration(); =20
=20
}catch(itk::ExceptionObject& Eo)
{ =20
AfxMessageBox(Eo.GetDescription());
return NULL;
}
=20
=20
TransformType::ParametersType lastParams =3D =
registration->GetLastTransformParameters(); // Get the results
=20
for (int pi=3D0;pi<12;pi++)
result[pi] =3D lastParams[pi];
=20
// result[12] =3D 0.0; result[13] =3D 0.0; result[14] =3D 0.0; =
result[15] =3D 1.0; =20
mat =3D new Matrix;
mat->Mat =3D (double**)malloc(sizeof (double) * 16);
mat->Rows =3D 4;
mat->Columns =3D 4;
memmove(mat->Mat, result, sizeof(double) * 16); =20
=20
return mat;
=20
//tranformation code
=20
if (pResult =3D=3D NULL || pResult->Mem =3D=3D NULL){=20
AfxMessageBox("Out volume can not be NULL");
return NULL;
}
=20
typedef BufferToImageConversion <short> ConverterType; // Volume - =
Image converetr
typedef ConverterType::ImageType UCharImage;
=20
typedef itk::AffineTransform<double, 3> TransformType;
itk::SmartPointer< TransformType> transform;
=20
typedef itk::ResampleImageFilter<UCharImage, UCharImage> =
ResampleFilter;
typedef UCharImage::SizeType ImageSizeType; =20
=20
typedef itk::LinearInterpolateImageFunction<UCharImage, double> =
InterpolatorType;
itk::SmartPointer< InterpolatorType> interpolator;
=
/////////////////////////////////////////////////////////////////////////=
//////////////////
//
// initializations
//
=
/////////////////////////////////////////////////////////////////////////=
//////////////////
=20
double result[12];
memcpy(result, mat->Mat, sizeof(double)*12);
ConverterType converter ;
ConverterType ::ImagePointer fixedImage =3D converter.GetImage(pFixed); =
=20
ConverterType ::ImagePointer movingImage =3D =
converter.GetImage(pMoving); =20
=20
=20
ResampleFilter::Pointer resampleFilter =3D ResampleFilter::New(); =20
=20
transform =3D TransformType::New();
interpolator=3DInterpolatorType::New();
TransformType::ParametersType =
params(transform->GetNumberOfParameters()) ;
=20
for (int pi=3D0;pi<12;pi++)
{
params[pi] =3D result[pi];
=20
}
=20
=20
UCharImage::Pointer outImage =3D NULL; =20
=20
transform->SetParameters( params);
=20
resampleFilter->SetInput(movingImage ); =20
resampleFilter ->SetTransform(transform.GetPointer());
resampleFilter ->SetInterpolator(interpolator.GetPointer());=20
resampleFilter->SetSize( =
fixedImage->GetLargestPossibleRegion().GetSize());
resampleFilter->SetOutputOrigin( fixedImage->GetOrigin() );
resampleFilter->SetOutputSpacing( fixedImage->GetSpacing() );
resampleFilter->SetDefaultPixelValue( 0 );
=20
try
{
resampleFilter->Update();
=20
}catch(itk::ExceptionObject &Eo){
=20
AfxMessageBox(Eo.GetDescription());
}
=20
outImage =3D resampleFilter->GetOutput();
=20
ConverterType::PixelType* pMem =3D outImage->GetBufferPointer();=20
=20
if (pMem =3D=3D NULL){
AfxMessageBox("NUll");
}=20
=20
=20
memcpy(pResult->Mem, pMem, pFixed->bytesPerVolume);
=20
return pResult;=20
}
=20
=20
Regards,
CSPL
----- Original Message -----=20
From: Lydia Ng <mailto:lng@insightful.com> =20
To: cspl <mailto:affable@hd2.dot.net.in> ; =
insight-users@public.kitware.com=20
Sent: Monday, November 11, 2002 11:17 PM
Subject: RE: [Insight-users] Query on Resampling
CSPL,
=20
Your code snippet looks ok - you might need to post more
of your code to see if everything is setup correctly.
I wonder if it is a templating issue?
Are you using the ITK registration components?=20
Have you set the transform and interpolation template parameter
of the ResampleImageFilter to the appropriate classes?
=20
-Lydia
-----Original Message-----
From: cspl [mailto:affable@hd2.dot.net.in]
Sent: Sunday, November 10, 2002 8:39 PM
To: insight-users@public.kitware.com
Subject: [Insight-users] Query on Resampling
Dear Friends,
I have some doubts regarding Resample image filter.After=20
=20
registring fixed and moving volumes, I applied registration output =
matrix on=20
=20
Moving volume.Output volume is not as expected.In some other tool,I=20
=20
have applied the same matrix on same moving volume.Output volume is=20
=20
totally different from resample output.Can u please tell what component=20
=20
is missing in resampling.I am enclosing the resample code. =20
=20
Please click on this url to see the pictures=20
=20
http://www.cspl.org/registration.jpg
=20
//code
ResampleFilter::Pointer resampleFilter =3D ResampleFilter::New(); =20
=20
transform =3D TransformType::New();
interpolator=3DInterpolatorType::New();
TransformType::ParametersType=20
=20
params(transform->GetNumberOfParameters()) ;
=20
for (int pi=3D0;pi<12;pi++)
{
params[pi] =3D result[pi];
=20
}
=20
=20
UCharImage::Pointer outImage =3D NULL; =20
=20
transform->SetParameters( params);
=20
resampleFilter->SetInput(movingImage ); =20
resampleFilter ->SetTransform(transform.GetPointer());
resampleFilter ->SetInterpolator(interpolator.GetPointer());=20
resampleFilter->SetSize(=20
=20
fixedImage->GetLargestPossibleRegion().GetSize());
resampleFilter->SetOutputOrigin( fixedImage->GetOrigin() );
resampleFilter->SetOutputSpacing( fixedImage->GetSpacing() );
=20
try
{
resampleFilter->Update();
=20
}catch(itk::ExceptionObject &Eo){
=20
AfxMessageBox(Eo.GetDescription());
}
=20
outImage =3D resampleFilter->GetOutput();
=20
Regards,
CSPL=20
------_=_NextPart_001_01C28A8C.56D95676
Content-Type: text/html;
charset="iso-8859-1"
Content-Transfer-Encoding: quoted-printable
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<HTML><HEAD>
<META HTTP-EQUIV=3D"Content-Type" CONTENT=3D"text/html; =
charset=3Diso-8859-1">
<META content=3D"MSHTML 5.00.3315.2870" name=3DGENERATOR>
<STYLE></STYLE>
</HEAD>
<BODY bgColor=3D#ffffff>
<DIV><FONT face=3D"Century Gothic" size=3D2><SPAN =
class=3D837483420-12112002>I had a=20
look through the snippet and I can't see anything that is obviously=20
wrong.</SPAN></FONT></DIV>
<DIV> </DIV>
<DIV><FONT face=3D"Century Gothic" size=3D2><SPAN =
class=3D837483420-12112002>The=20
MultiResMIRegistration App also uses the =
ResampleImageFilter</SPAN></FONT></DIV>
<DIV><FONT face=3D"Century Gothic" size=3D2><SPAN =
class=3D837483420-12112002>(with=20
affine transform and linear interpolation) and I just verified=20
that</SPAN></FONT></DIV>
<DIV><FONT face=3D"Century Gothic" size=3D2><SPAN =
class=3D837483420-12112002>it is=20
working ok - so I am confident that the ResampleImageFilter=20
is</SPAN></FONT></DIV>
<DIV><FONT face=3D"Century Gothic" size=3D2><SPAN =
class=3D837483420-12112002>working=20
as expected.</SPAN></FONT></DIV>
<DIV><FONT face=3D"Century Gothic" size=3D2><SPAN=20
class=3D837483420-12112002></SPAN></FONT> </DIV>
<DIV><FONT face=3D"Century Gothic" size=3D2><SPAN =
class=3D837483420-12112002>Other=20
suggestions:</SPAN></FONT></DIV>
<DIV><FONT face=3D"Century Gothic" size=3D2><SPAN =
class=3D837483420-12112002>- Are you=20
setting spacing correctly?</SPAN></FONT></DIV>
<DIV><FONT face=3D"Century Gothic" size=3D2><SPAN =
class=3D837483420-12112002>- Is it=20
an IO problem? Endian-ness?</SPAN></FONT></DIV>
<DIV><FONT face=3D"Century Gothic" size=3D2><SPAN=20
class=3D837483420-12112002></SPAN></FONT> </DIV>
<DIV><FONT face=3D"Century Gothic" size=3D2><SPAN =
class=3D837483420-12112002>I notice=20
that you convert your volume once for registration and once for=20
the</SPAN></FONT></DIV>
<DIV><FONT face=3D"Century Gothic" size=3D2><SPAN=20
class=3D837483420-12112002>resampling - are you doing anything different =
between=20
the two?</SPAN></FONT></DIV>
<DIV><FONT face=3D"Century Gothic" size=3D2><SPAN =
class=3D837483420-12112002>What=20
happens if you do the resampling with exactly the same into=20
as</SPAN></FONT></DIV>
<DIV><FONT face=3D"Century Gothic" size=3D2><SPAN =
class=3D837483420-12112002>for=20
registration (ie. without conversion)?</SPAN></FONT></DIV>
<DIV><FONT face=3D"Century Gothic" size=3D2><SPAN=20
class=3D837483420-12112002></SPAN></FONT> </DIV>
<DIV><FONT face=3D"Century Gothic" size=3D2><SPAN =
class=3D837483420-12112002>You code=20
contains other non-itk stuff so I can't compile and=20
debug.</SPAN></FONT></DIV>
<DIV><FONT face=3D"Century Gothic" size=3D2><SPAN =
class=3D837483420-12112002>If you=20
can reproduce the problem with ITK only code, package =
with</SPAN></FONT></DIV>
<DIV><FONT face=3D"Century Gothic" size=3D2><SPAN =
class=3D837483420-12112002>small=20
image then we can have a better chance in isolating the=20
problem.</SPAN></FONT></DIV>
<DIV><FONT face=3D"Century Gothic" size=3D2><SPAN=20
class=3D837483420-12112002></SPAN></FONT> </DIV>
<DIV><FONT face=3D"Century Gothic" size=3D2><SPAN =
class=3D837483420-12112002>-=20
Lydia </SPAN></FONT></DIV>
<DIV><FONT color=3D#808000 face=3D"Century Gothic" size=3D2><SPAN=20
class=3D837483420-12112002></SPAN></FONT> </DIV>
<BLOCKQUOTE dir=3Dltr=20
style=3D"BORDER-LEFT: #808000 2px solid; MARGIN-LEFT: 5px; MARGIN-RIGHT: =
0px; PADDING-LEFT: 5px">
<DIV align=3Dleft class=3DOutlookMessageHeader dir=3Dltr><FONT =
face=3DTahoma=20
size=3D2>-----Original Message-----<BR><B>From:</B> cspl=20
[mailto:affable@hd2.dot.net.in]<BR><B>Sent:</B> Monday, November 11, =
2002=20
10:26 PM<BR><B>To:</B> Lydia Ng<BR><B>Cc:</B>=20
insight-users@public.kitware.com<BR><B>Subject:</B> Re: =
[Insight-users] Query=20
on Resampling<BR><BR></DIV></FONT>
<DIV>
<DIV><FONT face=3DArial size=3D2>Hi!</FONT></DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV><FONT face=3DArial size=3D2>Thankyou for your =
response.</FONT></DIV>
<DIV><FONT face=3DArial size=3D2>Please find herewith the =
code snippet of=20
Registration and Resample volume.</FONT></DIV>
<DIV><FONT face=3DArial size=3D2>We set all the reigstration =
components. For=20
registration and resampling we set the transformation componentto =
Affine=20
Transform and interpolation component to Linear =
interpolation.</FONT></DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV><FONT face=3DArial size=3D2>Registration:</FONT> </DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV><FONT face=3DArial size=3D2>//Converting volume data to ITK image =
data<BR> typedef BufferToImageConversion<short>=20
ConverterType; <BR> typedef ConverterType::ImageType=20
FixedImageType;<BR> typedef=20
itk::StatisticsImageFilter<FixedImageType>=20
StatisticsImageFilterType;</FONT></DIV>
<DIV> </DIV><FONT face=3DArial size=3D2>
<DIV><BR> //Resgitartion components<BR> typedef=20
itk::AffineTransform<double> =
TransformType;<BR> typedef=20
itk::RegularStepGradientDescentOptimizer OptimizerType;=20
<BR> typedef itk::LinearInterpolateImageFunction=20
<FixedImageType, double> =
InterpolatorType;<BR> typedef=20
itk::ImageRegistrationMethod <FixedImageType , FixedImageType >=20
RegistrationType1; </DIV>
<DIV> </DIV>
<DIV><BR> //Get Fixed and Moving =
images<BR> =20
ConverterType converter;<BR> =
FixedImageType::Pointer=20
fixedImage =3D converter.GetImage(pFixed);<BR> =
FixedImageType::Pointer movingImage =3D=20
converter.GetImage(pMoving);<BR> </DIV>
<DIV> </DIV>
<DIV><BR> =20
=
////////////////////////////////////////////////////////////////////<BR>&=
nbsp; =20
//Set up the statistic filter <BR> =20
=
////////////////////////////////////////////////////////////////////<BR>&=
nbsp; =20
<BR> StatisticsImageFilterType::Pointer=20
fixedImageStatisticsFilter =3D=20
=
<BR> StatisticsImageFilter=
Type::New();<BR> =20
fixedImageStatisticsFilter->SetInput( fixedImage =
);<BR> =20
fixedImageStatisticsFilter->Update(); </DIV>
<DIV> </DIV>
<DIV> StatisticsImageFilterType::Pointer=20
movingImageStatisticsFilter =3D=20
=
<BR> StatisticsImageFilter=
Type::New();<BR> =20
movingImageStatisticsFilter->SetInput( movingImage =
);<BR> =20
movingImageStatisticsFilter->Update(); </DIV>
<DIV> </DIV>
<DIV> </DIV>
<DIV> </DIV>
=
<DIV> ////////////////////////////////////////////////////////=
////////////<BR> =20
//Set up the metric. <BR> =20
=
////////////////////////////////////////////////////////////////////<BR>&=
nbsp; =20
<BR> =
/*metric->SetMovingImageStandardDeviation( =20
movingImageStatisticsFilter->GetSigma() * 0.4 =
);<BR> =20
metric->SetFixedImageStandardDeviation(=20
fixedImageStatisticsFilter->GetSigma() * 0.4 =
);<BR> =20
=
metric->SetNumberOfSpatialSamples(m_spatialsamples);<BR> &n=
bsp;=20
metric->SetFixedImageRegion( fixedImage->GetBufferedRegion()=20
);*/<BR> =
<BR> RegistrationType1::Pointer =20
registration =3D RegistrationType1::New(); </DIV>
<DIV> </DIV>
<DIV> if (m_metrictype=3D=3D1) =20
<BR> {<BR> typedef=20
itk::MutualInformationImageToImageMetric<FixedImageType, =
FixedImageType>=20
=
MetricType;<BR> MetricType::Pointer &nb=
sp; =20
metric =3D=20
=
MetricType::New();<BR> metric->SetMovingImageStandard=
Deviation( =20
movingImageStatisticsFilter->GetSigma() * 0.4=20
);<BR> metric->SetFixedImageStandardDeviation(=20
fixedImageStatisticsFilter->GetSigma() * 0.4=20
=
);<BR> metric->SetNumberOfSpatialSamples(m_spatialsam=
ples);<BR> metric->SetFixedImageRegion(=20
fixedImage->GetBufferedRegion()=20
);<BR> registration->SetMetric( metric=20
);<BR> }<BR> else if (m_metrictype=3D=3D2)=20
<BR> {<BR> typedef=20
itk::NormalizedCorrelationImageToImageMetric<FixedImageType,=20
FixedImageType>=20
=
MetricType;<BR> MetricType::Pointer &nb=
sp; =20
metric =3D=20
=
MetricType::New();<BR> metric->SetFixedImageRegion(=20
fixedImage->GetBufferedRegion()=20
);<BR> registration->SetMetric( metric=20
);<BR> }<BR> else =20
<BR> {<BR> typedef=20
itk::MeanSquaresImageToImageMetric<FixedImageType, =
FixedImageType>=20
=
MetricType;<BR> MetricType::Pointer &nb=
sp; =20
metric =3D=20
=
MetricType::New();<BR> metric->SetFixedImageRegion(=20
fixedImage->GetBufferedRegion()=20
);<BR> registration->SetMetric( metric=20
);<BR> }</DIV>
<DIV> </DIV>
<DIV> <BR> Matrix *mat =3D=20
NULL;<BR> double result[12];</DIV>
<DIV> </DIV>
<DIV> =
TransformType::Pointer =20
transform =3D=20
TransformType::New();<BR> =20
OptimizerType::Pointer =20
optimizer =3D=20
OptimizerType::New();<BR> =20
InterpolatorType::Pointer interpolator =3D=20
InterpolatorType::New();<BR> </DIV>
<DIV> </DIV>
=
<DIV> //////////////////////////////////////////////////////////////=
//////<BR> //=20
Set up initial transform=20
=
parameters<BR> /////////////////////////////////////////////////////=
/////////////// =20
</DIV>
<DIV> </DIV>
<DIV> RegistrationType1::ParametersType initialParameters( =
=
<BR> transform->GetNumberOfParameters()=20
);</DIV>
<DIV> </DIV>
<DIV> //Set up transform parameters<BR> =
// =20
transform->SetIdentity();<BR> <BR> =20
TransformType::OutputVectorType shift;<BR> shift[0]=20
=3D1;<BR> shift[1] =3D0;<BR> shift[2] =
=3D0; //=20
translation of (10,15,20)<BR>// transform->Translate( =
shift=20
);</DIV>
<DIV> </DIV>
<DIV> TransformType::OutputVectorType =
axis;<BR> =20
axis[0] =3D 1.0;<BR> axis[1] =3D=20
0.0;<BR> axis[2] =3D 0.0; // =
axis=3D vector=20
parallel to the Z axis<BR> const double =
angle =3D=20
100;<BR> // transform->Rotate3D( axis, =
angle,0=20
);</DIV>
<DIV> </DIV>
<DIV> // TransformType::OutputVectorType =
shift;<BR> =20
shift[0] =3D1;<BR> shift[1] =3D 0;<BR> =
shift[2] =3D 0; =20
// translation of (10,15,20)<BR>// =
transform->Translate( shift=20
); <BR> /* TransformType::OutputVectorType =
shr;<BR> =20
shr[0]=3D0.5;<BR> shr[1]=3D0.0;<BR> =20
shr[2]=3D0.0;<BR> transform->Scale(shr,0);*/</DIV>
<DIV> </DIV>
<DIV> </DIV>
<DIV> </DIV>
<DIV> registration->SetInitialTransformParameters(=20
transform->GetParameters());<BR> <BR> =
<BR>/* =20
TransformType::ParametersType parameter; <BR> =20
parameter=3Dtransform->GetParameters();<BR> CString=20
str;<BR> for(int i=3D0;i<12;i++)<BR> =20
{ <BR> double=20
parameters[12];<BR> parameters[i]=3D=20
parameter[i];<BR> str.Format("%lf",=20
parameters[i]);<BR> =20
AfxMessageBox(str);<BR> <BR> =
}<BR> =20
*/</DIV>
<DIV> </DIV>
<DIV> =20
=
////////////////////////////////////////////////////////////////////<BR>&=
nbsp; =20
// Set up the optimizer.<BR> =20
=
//////////////////////////////////////////////////////////////////// =
; =20
<BR> <BR> typedef=20
OptimizerType::ScalesType ScalesType;<BR> ScalesType =
parametersScales( transform->GetNumberOfParameters()=20
);<BR> <BR> parametersScales.Fill( =
1.0=20
);<BR> <BR> double scale =3D 1.0 / =
vnl_math_sqr(m_translatescale);<BR> <BR> =
for=20
(int j =3D 9; j < 12; j++ )<BR> {</DIV>
<DIV> </DIV>
<DIV> parametersScales[j] =3D=20
scale;<BR> }<BR> =
<BR> =20
optimizer->SetScales( parametersScales =
);<BR> =20
<BR> =20
optimizer->SetNumberOfIterations(m_iterations);<BR> =20
<BR> <BR> =20
=
////////////////////////////////////////////////////////////////////<BR>&=
nbsp; =20
// Set up the registrator.<BR> =20
=
////////////////////////////////////////////////////////////////////</DIV=
>
<DIV> </DIV>
<DIV> // connect up the components<BR> <BR> =20
registration->SetOptimizer( optimizer );<BR> =20
registration->SetTransform( transform );<BR> =20
registration->SetFixedImage(fixedImage );<BR> =20
registration->SetMovingImage( movingImage );<BR> =20
registration->SetInterpolator( interpolator );</DIV>
<DIV> </DIV>
<DIV><BR> ProgressUpdate(75, " Start registration",=20
MRIVolumeName);<BR> try<BR> =
{<BR> =20
registration->StartRegistration(); </DIV>
<DIV> </DIV>
<DIV> }catch(itk::ExceptionObject& Eo)<BR> =
{ <BR> =20
AfxMessageBox(Eo.GetDescription());<BR> return=20
NULL;<BR> }</DIV>
<DIV> </DIV>
<DIV> <BR> TransformType::ParametersType=20
lastParams =3D =
registration->GetLastTransformParameters(); =20
// Get the results</DIV>
<DIV> </DIV>
<DIV><BR> for (int =
pi=3D0;pi<12;pi++)<BR> =20
result[pi] =3D lastParams[pi];</DIV>
<DIV> </DIV>
<DIV> // result[12] =3D 0.0; result[13] =3D =
0.0; result[14] =3D 0.0;=20
result[15] =3D 1.0; <BR> mat =3D new=20
Matrix;<BR> mat->Mat =3D (double**)malloc(sizeof =
(double) *=20
16);<BR> mat->Rows =3D 4;<BR> =
mat->Columns =3D=20
4;<BR> memmove(mat->Mat, result, sizeof(double) *=20
16); </DIV>
<DIV> </DIV>
<DIV> return mat;</DIV>
<DIV> </DIV>
<DIV>//tranformation code</DIV>
<DIV> </DIV>
<DIV>if (pResult =3D=3D NULL || pResult->Mem =3D=3D=20
NULL){ <BR> AfxMessageBox("Out volume can not be=20
NULL");<BR> return =
NULL;<BR> }<BR> <BR> typedef=20
BufferToImageConversion <short> ConverterType; // Volume - Image =
converetr<BR> typedef ConverterType::ImageType UCharImage;</DIV>
<DIV> </DIV>
<DIV> typedef itk::AffineTransform<double, 3>=20
TransformType;<BR> itk::SmartPointer< TransformType> =20
transform;</DIV>
<DIV> </DIV>
<DIV> typedef itk::ResampleImageFilter<UCharImage, =
UCharImage>=20
ResampleFilter;<BR> typedef UCharImage::SizeType =20
ImageSizeType; </DIV>
<DIV> </DIV>
<DIV> typedef itk::LinearInterpolateImageFunction<UCharImage,=20
double> InterpolatorType;<BR> itk::SmartPointer<=20
InterpolatorType> =20
=
interpolator;<BR> //////////////////////////////////////////////////=
/////////////////////////////////////////<BR> //<BR> // in=
itializations<BR> //<BR> //////////////////////////////////////=
/////////////////////////////////////////////////////</DIV>
<DIV> </DIV>
<DIV> double result[12];<BR> memcpy(result, mat->Mat,=20
sizeof(double)*12);<BR> ConverterType converter =
;<BR> ConverterType=20
::ImagePointer fixedImage =3D=20
=
converter.GetImage(pFixed); &nbs=
p;=20
<BR> ConverterType ::ImagePointer movingImage =3D=20
=
converter.GetImage(pMoving); &nb=
sp;=20
<BR> <BR> <BR> ResampleFilter::Pointer =
resampleFilter =3D=20
ResampleFilter::New(); </DIV>
<DIV> </DIV>
<DIV> transform =3D=20
=
TransformType::New();<BR> interpolator=3DInterpolatorType::New();<BR=
> TransformType::ParametersType=20
params(transform->GetNumberOfParameters()) ;<BR> <BR> for =
(int=20
pi=3D0;pi<12;pi++)<BR> {<BR> params[pi] =3D=20
result[pi];<BR> <BR> }<BR> </DIV>
<DIV> </DIV>
<DIV><BR> UCharImage::Pointer outImage =3D=20
=
NULL; <BR> <BR> transform->SetParameters(=20
params);</DIV>
<DIV> </DIV>
<DIV> resampleFilter->SetInput(movingImage=20
); <BR> resampleFilter=20
->SetTransform(transform.GetPointer());<BR> resampleFilter=20
->SetInterpolator(interpolator.GetPointer()); =
<BR> =20
resampleFilter->SetSize(=20
=
fixedImage->GetLargestPossibleRegion().GetSize());<BR> &nbs=
p;=20
resampleFilter->SetOutputOrigin( fixedImage->GetOrigin()=20
);<BR> resampleFilter->SetOutputSpacing(=20
fixedImage->GetSpacing() );<BR> =20
resampleFilter->SetDefaultPixelValue( 0 );</DIV>
<DIV> </DIV>
=
<DIV> try<BR> {<BR> resampleFilter->Update();</DI=
V>
<DIV> </DIV>
<DIV> }catch(itk::ExceptionObject &Eo){</DIV>
<DIV> </DIV>
<DIV> AfxMessageBox(Eo.GetDescription());<BR> }</DIV>
<DIV> </DIV>
<DIV> outImage =3D =
resampleFilter->GetOutput();</DIV>
<DIV> </DIV>
<DIV> ConverterType::PixelType* pMem =3D=20
outImage->GetBufferPointer(); </DIV>
<DIV> </DIV>
<DIV><BR> if (pMem =3D=3D=20
=
NULL){<BR> AfxMessageBox("NUll");<BR> } <BR> </=
DIV>
<DIV> </DIV>
<DIV> memcpy(pResult->Mem, pMem, =
pFixed->bytesPerVolume);</DIV>
<DIV> </DIV>
<DIV><BR> return pResult; <BR>}</DIV>
<DIV> </DIV>
<DIV> </DIV>
<DIV>Regards,</DIV>
<DIV>CSPL</FONT></DIV></DIV>
<BLOCKQUOTE dir=3Dltr=20
style=3D"BORDER-LEFT: #000000 2px solid; MARGIN-LEFT: 5px; =
MARGIN-RIGHT: 0px; PADDING-LEFT: 5px; PADDING-RIGHT: 0px">
<DIV style=3D"FONT: 10pt arial">----- Original Message ----- </DIV>
<DIV=20
style=3D"BACKGROUND: #e4e4e4; FONT: 10pt arial; font-color: =
black"><B>From:</B>=20
<A href=3D"mailto:lng@insightful.com" =
title=3Dlng@insightful.com>Lydia Ng</A>=20
</DIV>
<DIV style=3D"FONT: 10pt arial"><B>To:</B> <A=20
href=3D"mailto:affable@hd2.dot.net.in" =
title=3Daffable@hd2.dot.net.in>cspl</A> ;=20
<A href=3D"mailto:insight-users@public.kitware.com"=20
=
title=3Dinsight-users@public.kitware.com>insight-users@public.kitware.com=
</A>=20
</DIV>
<DIV style=3D"FONT: 10pt arial"><B>Sent:</B> Monday, November 11, =
2002 11:17=20
PM</DIV>
<DIV style=3D"FONT: 10pt arial"><B>Subject:</B> RE: [Insight-users] =
Query on=20
Resampling</DIV>
<DIV><BR></DIV>
<DIV>
<DIV><FONT color=3D#0000ff face=3DArial size=3D2><SPAN=20
class=3D494524517-11112002>CSPL,</SPAN></FONT></DIV>
<DIV><FONT color=3D#0000ff face=3DArial size=3D2><SPAN=20
class=3D494524517-11112002></SPAN></FONT> </DIV>
<DIV><FONT color=3D#0000ff face=3DArial size=3D2><SPAN=20
class=3D494524517-11112002>Your code snippet looks ok - you might =
need to post=20
more</SPAN></FONT></DIV>
<DIV><FONT color=3D#0000ff face=3DArial size=3D2><SPAN =
class=3D494524517-11112002>of=20
your code to see if everything is setup =
correctly.</SPAN></FONT></DIV>
<DIV><FONT color=3D#0000ff face=3DArial size=3D2><SPAN =
class=3D494524517-11112002>I=20
wonder if it is a templating issue?</SPAN></FONT></DIV>
<DIV><FONT color=3D#0000ff face=3DArial size=3D2><SPAN=20
class=3D494524517-11112002>Are you using the ITK registration =
components?=20
</SPAN></FONT></DIV>
<DIV><FONT color=3D#0000ff face=3DArial size=3D2><SPAN=20
class=3D494524517-11112002>Have you set the transform and =
interpolation=20
template parameter</SPAN></FONT></DIV>
<DIV><FONT color=3D#0000ff face=3DArial size=3D2><SPAN =
class=3D494524517-11112002>of=20
the ResampleImageFilter to the appropriate =
classes?</SPAN></FONT></DIV>
<DIV><FONT color=3D#0000ff face=3DArial size=3D2><SPAN=20
class=3D494524517-11112002></SPAN></FONT> </DIV>
<DIV><FONT color=3D#0000ff face=3DArial size=3D2><SPAN=20
class=3D494524517-11112002>-Lydia</SPAN></FONT></DIV></DIV>
<BLOCKQUOTE dir=3Dltr style=3D"MARGIN-RIGHT: 0px">
<DIV align=3Dleft class=3DOutlookMessageHeader dir=3Dltr><FONT =
face=3DTahoma=20
size=3D2>-----Original Message-----<BR><B>From:</B> cspl=20
[mailto:affable@hd2.dot.net.in]<BR><B>Sent:</B> Sunday, November =
10, 2002=20
8:39 PM<BR><B>To:</B> =
insight-users@public.kitware.com<BR><B>Subject:</B>=20
[Insight-users] Query on Resampling<BR><BR></FONT></DIV>
<DIV>
<DIV><FONT face=3DArial size=3D2>Dear Friends,</FONT></DIV>
<DIV><FONT face=3DArial size=3D2><BR> I have some doubts =
regarding=20
Resample image filter.After </FONT></DIV>
<DIV> </DIV>
<DIV><FONT face=3DArial size=3D2>registring fixed and moving =
volumes, I=20
applied registration output matrix on </FONT></DIV>
<DIV> </DIV>
<DIV><FONT face=3DArial size=3D2>Moving volume.Output volume is =
not as=20
expected.In some other tool,I </FONT></DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV><FONT face=3DArial size=3D2>have applied the same matrix on =
same moving=20
volume.Output volume is </FONT></DIV>
<DIV> </DIV>
<DIV><FONT face=3DArial size=3D2>totally different from resample =
output.Can u=20
please tell what component </FONT></DIV>
<DIV> </DIV>
<DIV><FONT face=3DArial size=3D2>is missing in resampling.I am =
enclosing the=20
resample code. </FONT></DIV>
<DIV> </DIV>
<DIV><FONT face=3DArial size=3D2>Please click on this url to see =
the pictures=20
</FONT></DIV>
<DIV> </DIV>
<DIV><FONT face=3DArial size=3D2><A=20
=
href=3D"http://www.cspl.org/registration.jpg">http://www.cspl.org/registr=
ation.jpg</A></FONT></DIV>
<DIV><FONT face=3DArial size=3D2></FONT> </DIV>
<DIV><FONT face=3DArial size=3D2>//code<BR>ResampleFilter::Pointer =
resampleFilter =3D ResampleFilter::New(); </FONT></DIV>
<DIV> </DIV>
<DIV><FONT face=3DArial size=3D2> transform =3D=20
=
TransformType::New();<BR> interpolator=3DInterpolatorType::New();<BR=
> TransformType::ParametersType=20
</FONT></DIV>
<DIV> </DIV>
<DIV><FONT face=3DArial =
size=3D2>params(transform->GetNumberOfParameters())=20
;<BR> <BR> for (int=20
pi=3D0;pi<12;pi++)<BR> {<BR> params[pi] =
=3D=20
result[pi];<BR> =
<BR> }<BR> </FONT></DIV>
<DIV> </DIV><FONT face=3DArial size=3D2>
<DIV><BR> UCharImage::Pointer outImage =3D=20
=
NULL; <BR> <BR> transform->SetParameters(=20
params);</DIV>
<DIV> </DIV>
<DIV><BR> =
resampleFilter->SetInput(movingImage=20
); <BR> resampleFilter=20
->SetTransform(transform.GetPointer());<BR> resampleFilter =
->SetInterpolator(interpolator.GetPointer()); =
<BR> =20
resampleFilter->SetSize( </DIV>
<DIV> </DIV>
=
<DIV>fixedImage->GetLargestPossibleRegion().GetSize());<BR>  =
; =20
resampleFilter->SetOutputOrigin( fixedImage->GetOrigin()=20
);<BR> resampleFilter->SetOutputSpacing(=20
fixedImage->GetSpacing() );</DIV>
<DIV> </DIV>
=
<DIV><BR> try<BR> {<BR> resampleFilter->Update();=
</DIV>
<DIV> </DIV>
<DIV> }catch(itk::ExceptionObject &Eo){</DIV>
<DIV> </DIV>
=
<DIV> AfxMessageBox(Eo.GetDescription());<BR> }</DIV>
<DIV> </DIV>
<DIV> outImage =3D =
resampleFilter->GetOutput();</DIV>
<DIV> </DIV>
<DIV><BR>Regards,<BR>CSPL=20
<BR></FONT></DIV></DIV></BLOCKQUOTE></BLOCKQUOTE></BLOCKQUOTE></BODY></HT=
ML>
------_=_NextPart_001_01C28A8C.56D95676--