[Insight-users] Query on Resampling

cspl affable@hd2.dot.net.in
Tue, 12 Nov 2002 11:56:16 +0530


This is a multi-part message in MIME format.

------=_NextPart_000_0088_01C28A42.81582F70
Content-Type: text/plain;
	charset="iso-8859-1"
Content-Transfer-Encoding: quoted-printable

Hi!

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.


Registration:=20

//Converting volume data to ITK image data
  typedef BufferToImageConversion<short> ConverterType;=20
  typedef ConverterType::ImageType FixedImageType;
  typedef itk::StatisticsImageFilter<FixedImageType> =
StatisticsImageFilterType;


  //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


  //Get Fixed and Moving images
     ConverterType converter;
     FixedImageType::Pointer fixedImage =3D converter.GetImage(pFixed);
     FixedImageType::Pointer movingImage =3D =
converter.GetImage(pMoving);
=20


    ////////////////////////////////////////////////////////////////////
    //Set up the statistic filter=20
    ////////////////////////////////////////////////////////////////////
   =20
    StatisticsImageFilterType::Pointer fixedImageStatisticsFilter =3D=20
        StatisticsImageFilterType::New();
    fixedImageStatisticsFilter->SetInput( fixedImage );
    fixedImageStatisticsFilter->Update();=20

    StatisticsImageFilterType::Pointer movingImageStatisticsFilter =3D=20
        StatisticsImageFilterType::New();
    movingImageStatisticsFilter->SetInput( movingImage );
    movingImageStatisticsFilter->Update();=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

  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
  Matrix *mat =3D NULL;
  double result[12];

    TransformType::Pointer      transform     =3D TransformType::New();
    OptimizerType::Pointer      optimizer     =3D OptimizerType::New();
    InterpolatorType::Pointer   interpolator  =3D =
InterpolatorType::New();
  =20

 ////////////////////////////////////////////////////////////////////
 // Set up initial transform parameters
 ////////////////////////////////////////////////////////////////////    =


   RegistrationType1::ParametersType initialParameters(=20
     transform->GetNumberOfParameters() );

   //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 );

   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 );

 //  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);*/



   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
   }
   */

    ////////////////////////////////////////////////////////////////////
    // 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++ )
    {

     parametersScales[j] =3D scale;
    }
   =20
    optimizer->SetScales( parametersScales );
    =20
    optimizer->SetNumberOfIterations(m_iterations);
  =20
  =20
   ////////////////////////////////////////////////////////////////////
   // Set up the registrator.
   ////////////////////////////////////////////////////////////////////

  // connect up the components
 =20
   registration->SetOptimizer( optimizer );
   registration->SetTransform( transform );
   registration->SetFixedImage(fixedImage );
   registration->SetMovingImage( movingImage  );
   registration->SetInterpolator( interpolator );


   ProgressUpdate(75, " Start registration", MRIVolumeName);
   try
   {
    registration->StartRegistration();    =20

   }catch(itk::ExceptionObject& Eo)
   {  =20
   AfxMessageBox(Eo.GetDescription());
   return NULL;
   }

  =20
   TransformType::ParametersType lastParams =3D =
registration->GetLastTransformParameters();    // Get the results


   for (int pi=3D0;pi<12;pi++)
    result[pi] =3D lastParams[pi];

  // 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

   return mat;

//tranformation code

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;

 typedef itk::AffineTransform<double, 3> TransformType;
 itk::SmartPointer< TransformType>  transform;

 typedef itk::ResampleImageFilter<UCharImage, UCharImage>  =
ResampleFilter;
 typedef UCharImage::SizeType  ImageSizeType; =20

 typedef itk::LinearInterpolateImageFunction<UCharImage, double> =
InterpolatorType;
 itk::SmartPointer< InterpolatorType>  interpolator;
 =
/////////////////////////////////////////////////////////////////////////=
//////////////////
 //
 // initializations
 //
 =
/////////////////////////////////////////////////////////////////////////=
//////////////////

 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

 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


 UCharImage::Pointer outImage  =3D NULL; =20
 =20
 transform->SetParameters( params);

    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 );

 try
 {
  resampleFilter->Update();

 }catch(itk::ExceptionObject &Eo){

  AfxMessageBox(Eo.GetDescription());
 }

 outImage   =3D resampleFilter->GetOutput();

 ConverterType::PixelType* pMem =3D outImage->GetBufferPointer();=20


 if (pMem =3D=3D NULL){
  AfxMessageBox("NUll");
 }=20
=20

 memcpy(pResult->Mem, pMem, pFixed->bytesPerVolume);


 return pResult;=20
}


Regards,
CSPL
  ----- Original Message -----=20
  From: Lydia Ng=20
  To: cspl ; insight-users@public.kitware.com=20
  Sent: Monday, November 11, 2002 11:17 PM
  Subject: RE: [Insight-users] Query on Resampling


  CSPL,

  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?

  -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

    registring fixed and moving volumes, I applied registration output =
matrix on=20

    Moving volume.Output volume is not as expected.In some other tool,I=20

    have applied the same matrix on same moving volume.Output volume is=20

    totally different from resample output.Can u please tell what =
component=20

    is missing in resampling.I am enclosing the resample code. =20

    Please click on this url to see the pictures=20

    http://www.cspl.org/registration.jpg

    //code
    ResampleFilter::Pointer resampleFilter =3D ResampleFilter::New(); =20

     transform =3D TransformType::New();
     interpolator=3DInterpolatorType::New();
     TransformType::ParametersType=20

    params(transform->GetNumberOfParameters()) ;
    =20
     for (int pi=3D0;pi<12;pi++)
     {
        params[pi] =3D result[pi];
       =20
     }
    =20


     UCharImage::Pointer outImage  =3D NULL; =20
     =20
     transform->SetParameters( params);


        resampleFilter->SetInput(movingImage );   =20
     resampleFilter ->SetTransform(transform.GetPointer());
     resampleFilter ->SetInterpolator(interpolator.GetPointer());=20
        resampleFilter->SetSize(=20

    fixedImage->GetLargestPossibleRegion().GetSize());
        resampleFilter->SetOutputOrigin( fixedImage->GetOrigin() );
        resampleFilter->SetOutputSpacing( fixedImage->GetSpacing() );


     try
     {
      resampleFilter->Update();

     }catch(itk::ExceptionObject &Eo){

      AfxMessageBox(Eo.GetDescription());
     }

     outImage   =3D resampleFilter->GetOutput();


    Regards,
    CSPL=20


------=_NextPart_000_0088_01C28A42.81582F70
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=3DContent-Type content=3D"text/html; =
charset=3Diso-8859-1">
<META content=3D"MSHTML 5.50.4134.600" name=3DGENERATOR>
<STYLE></STYLE>
</HEAD>
<BODY bgColor=3D#ffffff>
<DIV>
<DIV><FONT face=3DArial size=3D2>Hi!</FONT></DIV>
<DIV><FONT face=3DArial size=3D2></FONT>&nbsp;</DIV>
<DIV><FONT face=3DArial size=3D2>Thankyou for your =
response.</FONT></DIV>
<DIV><FONT face=3DArial size=3D2>Please find herewith the =
code&nbsp;snippet of=20
Registration and&nbsp;Resample volume.</FONT></DIV>
<DIV><FONT face=3DArial size=3D2>We&nbsp;set all the reigstration =
components. For=20
registration and resampling we set the transformation componentto Affine =

Transform and interpolation component to Linear =
interpolation.</FONT></DIV>
<DIV><FONT face=3DArial size=3D2></FONT>&nbsp;</DIV>
<DIV><FONT face=3DArial size=3D2></FONT>&nbsp;</DIV>
<DIV><FONT face=3DArial size=3D2>Registration:</FONT>&nbsp;</DIV>
<DIV><FONT face=3DArial size=3D2></FONT>&nbsp;</DIV>
<DIV><FONT face=3DArial size=3D2>//Converting volume data to ITK image=20
data<BR>&nbsp;&nbsp;typedef BufferToImageConversion&lt;short&gt;=20
ConverterType;&nbsp;<BR>&nbsp;&nbsp;typedef ConverterType::ImageType=20
FixedImageType;<BR>&nbsp;&nbsp;typedef=20
itk::StatisticsImageFilter&lt;FixedImageType&gt;=20
StatisticsImageFilterType;</FONT></DIV>
<DIV>&nbsp;</DIV><FONT face=3DArial size=3D2>
<DIV><BR>&nbsp;&nbsp;//Resgitartion components<BR>&nbsp;&nbsp;typedef=20
itk::AffineTransform&lt;double&gt; TransformType;<BR>&nbsp;&nbsp;typedef =

itk::RegularStepGradientDescentOptimizer OptimizerType; =
<BR>&nbsp;&nbsp;typedef=20
itk::LinearInterpolateImageFunction &lt;FixedImageType, double&gt;=20
InterpolatorType;<BR>&nbsp;&nbsp;typedef itk::ImageRegistrationMethod=20
&lt;FixedImageType , FixedImageType &gt; RegistrationType1;&nbsp; </DIV>
<DIV>&nbsp;</DIV>
<DIV><BR>&nbsp;&nbsp;//Get Fixed and Moving =
images<BR>&nbsp;&nbsp;&nbsp;&nbsp;=20
ConverterType converter;<BR>&nbsp;&nbsp;&nbsp;&nbsp; =
FixedImageType::Pointer=20
fixedImage =3D converter.GetImage(pFixed);<BR>&nbsp;&nbsp;&nbsp;&nbsp;=20
FixedImageType::Pointer movingImage =3D=20
converter.GetImage(pMoving);<BR>&nbsp;</DIV>
<DIV>&nbsp;</DIV>
<DIV><BR>&nbsp;&nbsp;&nbsp;=20
////////////////////////////////////////////////////////////////////<BR>&=
nbsp;&nbsp;&nbsp;=20
//Set up the statistic filter <BR>&nbsp;&nbsp;&nbsp;=20
////////////////////////////////////////////////////////////////////<BR>&=
nbsp;&nbsp;&nbsp;=20
<BR>&nbsp;&nbsp;&nbsp; StatisticsImageFilterType::Pointer=20
fixedImageStatisticsFilter =3D=20
<BR>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;StatisticsImageFilter=
Type::New();<BR>&nbsp;&nbsp;&nbsp;=20
fixedImageStatisticsFilter-&gt;SetInput( fixedImage =
);<BR>&nbsp;&nbsp;&nbsp;=20
fixedImageStatisticsFilter-&gt;Update(); </DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;&nbsp;&nbsp; StatisticsImageFilterType::Pointer=20
movingImageStatisticsFilter =3D=20
<BR>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;StatisticsImageFilter=
Type::New();<BR>&nbsp;&nbsp;&nbsp;=20
movingImageStatisticsFilter-&gt;SetInput( movingImage =
);<BR>&nbsp;&nbsp;&nbsp;=20
movingImageStatisticsFilter-&gt;Update(); </DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;&nbsp;////////////////////////////////////////////////////////=
////////////<BR>&nbsp;&nbsp;&nbsp;=20
//Set up the metric. <BR>&nbsp;&nbsp;&nbsp;=20
////////////////////////////////////////////////////////////////////<BR>&=
nbsp;&nbsp;&nbsp;=20
<BR>&nbsp;&nbsp;&nbsp; =
/*metric-&gt;SetMovingImageStandardDeviation(&nbsp;=20
movingImageStatisticsFilter-&gt;GetSigma() * 0.4 =
);<BR>&nbsp;&nbsp;&nbsp;=20
metric-&gt;SetFixedImageStandardDeviation(=20
fixedImageStatisticsFilter-&gt;GetSigma() * 0.4 );<BR>&nbsp;&nbsp;&nbsp; =

metric-&gt;SetNumberOfSpatialSamples(m_spatialsamples);<BR>&nbsp;&nbsp;&n=
bsp;=20
metric-&gt;SetFixedImageRegion( fixedImage-&gt;GetBufferedRegion()=20
);*/<BR>&nbsp;&nbsp;&nbsp; =
<BR>&nbsp;&nbsp;RegistrationType1::Pointer&nbsp;=20
registration&nbsp; =3D RegistrationType1::New();&nbsp;&nbsp;</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;&nbsp;if (m_metrictype=3D=3D1)&nbsp;=20
<BR>&nbsp;&nbsp;{<BR>&nbsp;&nbsp;&nbsp;typedef=20
itk::MutualInformationImageToImageMetric&lt;FixedImageType, =
FixedImageType&gt;=20
MetricType;<BR>&nbsp;&nbsp;&nbsp;MetricType::Pointer&nbsp;&nbsp;&nbsp;&nb=
sp;&nbsp;&nbsp;&nbsp;&nbsp;=20
metric&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; =3D=20
MetricType::New();<BR>&nbsp;&nbsp;&nbsp;metric-&gt;SetMovingImageStandard=
Deviation(&nbsp;=20
movingImageStatisticsFilter-&gt;GetSigma() * 0.4=20
);<BR>&nbsp;&nbsp;&nbsp;metric-&gt;SetFixedImageStandardDeviation(=20
fixedImageStatisticsFilter-&gt;GetSigma() * 0.4=20
);<BR>&nbsp;&nbsp;&nbsp;metric-&gt;SetNumberOfSpatialSamples(m_spatialsam=
ples);<BR>&nbsp;&nbsp;&nbsp;metric-&gt;SetFixedImageRegion(=20
fixedImage-&gt;GetBufferedRegion()=20
);<BR>&nbsp;&nbsp;&nbsp;registration-&gt;SetMetric( metric=20
);<BR>&nbsp;&nbsp;}<BR>&nbsp;&nbsp;else if (m_metrictype=3D=3D2)=20
<BR>&nbsp;&nbsp;{<BR>&nbsp;&nbsp;&nbsp;typedef=20
itk::NormalizedCorrelationImageToImageMetric&lt;FixedImageType,=20
FixedImageType&gt;=20
MetricType;<BR>&nbsp;&nbsp;&nbsp;MetricType::Pointer&nbsp;&nbsp;&nbsp;&nb=
sp;&nbsp;&nbsp;&nbsp;&nbsp;=20
metric&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; =3D=20
MetricType::New();<BR>&nbsp;&nbsp;&nbsp;metric-&gt;SetFixedImageRegion(=20
fixedImage-&gt;GetBufferedRegion()=20
);<BR>&nbsp;&nbsp;&nbsp;registration-&gt;SetMetric( metric=20
);<BR>&nbsp;&nbsp;}<BR>&nbsp;&nbsp;else&nbsp;=20
<BR>&nbsp;&nbsp;{<BR>&nbsp;&nbsp;&nbsp;typedef=20
itk::MeanSquaresImageToImageMetric&lt;FixedImageType, FixedImageType&gt; =

MetricType;<BR>&nbsp;&nbsp;&nbsp;MetricType::Pointer&nbsp;&nbsp;&nbsp;&nb=
sp;&nbsp;&nbsp;&nbsp;&nbsp;=20
metric&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; =3D=20
MetricType::New();<BR>&nbsp;&nbsp;&nbsp;metric-&gt;SetFixedImageRegion(=20
fixedImage-&gt;GetBufferedRegion()=20
);<BR>&nbsp;&nbsp;&nbsp;registration-&gt;SetMetric( metric=20
);<BR>&nbsp;&nbsp;}</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<BR>&nbsp;&nbsp;Matrix *mat =3D=20
NULL;<BR>&nbsp;&nbsp;double result[12];</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;&nbsp;&nbsp; =
TransformType::Pointer&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;=20
transform&nbsp;&nbsp;&nbsp;&nbsp; =3D =
TransformType::New();<BR>&nbsp;&nbsp;&nbsp;=20
OptimizerType::Pointer&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;=20
optimizer&nbsp;&nbsp;&nbsp;&nbsp; =3D =
OptimizerType::New();<BR>&nbsp;&nbsp;&nbsp;=20
InterpolatorType::Pointer&nbsp;&nbsp; interpolator&nbsp; =3D=20
InterpolatorType::New();<BR>&nbsp;&nbsp; </DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;//////////////////////////////////////////////////////////////=
//////<BR>&nbsp;//=20
Set up initial transform=20
parameters<BR>&nbsp;/////////////////////////////////////////////////////=
///////////////&nbsp;&nbsp;&nbsp;=20
</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;&nbsp; RegistrationType1::ParametersType initialParameters(=20
<BR>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;transform-&gt;GetNumberOfParameters() =
);</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;&nbsp; //Set up transform parameters<BR>&nbsp;&nbsp; //&nbsp; =

transform-&gt;SetIdentity();<BR>&nbsp; <BR>&nbsp;&nbsp;=20
TransformType::OutputVectorType shift;<BR>&nbsp;&nbsp; shift[0]=20
=3D1;<BR>&nbsp;&nbsp; shift[1] =3D0;<BR>&nbsp;&nbsp; shift[2] =
=3D0;&nbsp; //=20
translation of (10,15,20)<BR>//&nbsp;&nbsp; transform-&gt;Translate( =
shift=20
);</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;&nbsp; TransformType::OutputVectorType axis;<BR>&nbsp;&nbsp; =
axis[0]=20
=3D 1.0;<BR>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; axis[1] =3D=20
0.0;<BR>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; axis[2] =3D 0.0;&nbsp; // axis=3D =
vector=20
parallel to the Z axis<BR>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; const double =
angle =3D=20
100;<BR>&nbsp; //&nbsp;&nbsp;&nbsp; transform-&gt;Rotate3D( axis, =
angle,0=20
);</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;//&nbsp; TransformType::OutputVectorType =
shift;<BR>&nbsp;&nbsp;=20
shift[0] =3D1;<BR>&nbsp;&nbsp; shift[1] =3D 0;<BR>&nbsp;&nbsp; shift[2] =
=3D 0;&nbsp;=20
// translation of (10,15,20)<BR>//&nbsp;&nbsp; transform-&gt;Translate( =
shift );=20
<BR>&nbsp;/*&nbsp; TransformType::OutputVectorType shr;<BR>&nbsp;&nbsp;=20
shr[0]=3D0.5;<BR>&nbsp;&nbsp; shr[1]=3D0.0;<BR>&nbsp;&nbsp;=20
shr[2]=3D0.0;<BR>&nbsp;&nbsp; transform-&gt;Scale(shr,0);*/</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;&nbsp; registration-&gt;SetInitialTransformParameters(=20
transform-&gt;GetParameters());<BR>&nbsp;<BR>&nbsp;&nbsp; =
<BR>/*&nbsp;&nbsp;=20
TransformType::ParametersType parameter; <BR>&nbsp;&nbsp;=20
parameter=3Dtransform-&gt;GetParameters();<BR>&nbsp;&nbsp; CString=20
str;<BR>&nbsp;&nbsp; for(int i=3D0;i&lt;12;i++)<BR>&nbsp;&nbsp;=20
{&nbsp;&nbsp;&nbsp;&nbsp; <BR>&nbsp;&nbsp;&nbsp; double=20
parameters[12];<BR>&nbsp;&nbsp;&nbsp; parameters[i]=3D=20
parameter[i];<BR>&nbsp;&nbsp;&nbsp; str.Format("%lf",=20
parameters[i]);<BR>&nbsp;&nbsp;&nbsp; =
AfxMessageBox(str);<BR>&nbsp;&nbsp;&nbsp;=20
<BR>&nbsp;&nbsp; }<BR>&nbsp;&nbsp; */</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;&nbsp;&nbsp;=20
////////////////////////////////////////////////////////////////////<BR>&=
nbsp;&nbsp;&nbsp;=20
// Set up the optimizer.<BR>&nbsp;&nbsp;&nbsp;=20
////////////////////////////////////////////////////////////////////&nbsp=
;&nbsp;&nbsp;=20
<BR>&nbsp;&nbsp;&nbsp; <BR>&nbsp;&nbsp;&nbsp; typedef =
OptimizerType::ScalesType=20
ScalesType;<BR>&nbsp;&nbsp;&nbsp; ScalesType parametersScales(=20
transform-&gt;GetNumberOfParameters() );<BR>&nbsp;&nbsp;&nbsp;=20
<BR>&nbsp;&nbsp;&nbsp; parametersScales.Fill( 1.0 =
);<BR>&nbsp;&nbsp;&nbsp;=20
<BR>&nbsp;&nbsp;&nbsp; double scale =3D 1.0 /=20
vnl_math_sqr(m_translatescale);<BR>&nbsp;&nbsp; <BR>&nbsp;&nbsp;&nbsp; =
for=20
(int&nbsp; j =3D 9; j &lt; 12; j++ )<BR>&nbsp;&nbsp;&nbsp; {</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;&nbsp;&nbsp;&nbsp; parametersScales[j] =3D =
scale;<BR>&nbsp;&nbsp;&nbsp;=20
}<BR>&nbsp;&nbsp;&nbsp; <BR>&nbsp;&nbsp;&nbsp; optimizer-&gt;SetScales(=20
parametersScales );<BR>&nbsp;&nbsp;&nbsp;&nbsp; <BR>&nbsp;&nbsp;&nbsp;=20
optimizer-&gt;SetNumberOfIterations(m_iterations);<BR>&nbsp;&nbsp;=20
<BR>&nbsp;&nbsp; <BR>&nbsp;&nbsp;=20
////////////////////////////////////////////////////////////////////<BR>&=
nbsp;&nbsp;=20
// Set up the registrator.<BR>&nbsp;&nbsp;=20
////////////////////////////////////////////////////////////////////</DIV=
>
<DIV>&nbsp;</DIV>
<DIV>&nbsp; // connect up the components<BR>&nbsp; <BR>&nbsp;&nbsp;=20
registration-&gt;SetOptimizer( optimizer );<BR>&nbsp;&nbsp;=20
registration-&gt;SetTransform( transform );<BR>&nbsp;&nbsp;=20
registration-&gt;SetFixedImage(fixedImage );<BR>&nbsp;&nbsp;=20
registration-&gt;SetMovingImage( movingImage&nbsp; );<BR>&nbsp;&nbsp;=20
registration-&gt;SetInterpolator( interpolator );</DIV>
<DIV>&nbsp;</DIV>
<DIV><BR>&nbsp;&nbsp; ProgressUpdate(75, " Start registration",=20
MRIVolumeName);<BR>&nbsp;&nbsp; try<BR>&nbsp;&nbsp; =
{<BR>&nbsp;&nbsp;&nbsp;=20
registration-&gt;StartRegistration();&nbsp;&nbsp;&nbsp;&nbsp; </DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;&nbsp; }catch(itk::ExceptionObject&amp; Eo)<BR>&nbsp;&nbsp;=20
{&nbsp;&nbsp;&nbsp;<BR>&nbsp;&nbsp;=20
AfxMessageBox(Eo.GetDescription());<BR>&nbsp;&nbsp; return =
NULL;<BR>&nbsp;&nbsp;=20
}</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;&nbsp;&nbsp;<BR>&nbsp;&nbsp; TransformType::ParametersType =
lastParams=20
=3D registration-&gt;GetLastTransformParameters();&nbsp;&nbsp;&nbsp; // =
Get the=20
results</DIV>
<DIV>&nbsp;</DIV>
<DIV><BR>&nbsp;&nbsp; for (int =
pi=3D0;pi&lt;12;pi++)<BR>&nbsp;&nbsp;&nbsp;=20
result[pi] =3D lastParams[pi];</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp; // result[12] =3D 0.0;&nbsp;result[13] =3D =
0.0;&nbsp;result[14] =3D 0.0;=20
result[15] =3D 1.0;&nbsp;&nbsp;&nbsp;<BR>&nbsp;&nbsp; mat =3D new=20
Matrix;<BR>&nbsp;&nbsp; mat-&gt;Mat =3D (double**)malloc(sizeof (double) =
*=20
16);<BR>&nbsp;&nbsp; mat-&gt;Rows =3D 4;<BR>&nbsp;&nbsp; mat-&gt;Columns =
=3D=20
4;<BR>&nbsp;&nbsp; memmove(mat-&gt;Mat, result, sizeof(double) *=20
16);&nbsp;&nbsp;&nbsp;</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;&nbsp; return mat;</DIV>
<DIV>&nbsp;</DIV>
<DIV>//tranformation code</DIV>
<DIV>&nbsp;</DIV>
<DIV>if (pResult =3D=3D NULL || pResult-&gt;Mem =3D=3D=20
NULL){&nbsp;<BR>&nbsp;&nbsp;AfxMessageBox("Out volume can not be=20
NULL");<BR>&nbsp;&nbsp;return =
NULL;<BR>&nbsp;}<BR>&nbsp;<BR>&nbsp;typedef=20
BufferToImageConversion &lt;short&gt; ConverterType; // Volume - Image=20
converetr<BR>&nbsp;typedef ConverterType::ImageType UCharImage;</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;typedef itk::AffineTransform&lt;double, 3&gt;=20
TransformType;<BR>&nbsp;itk::SmartPointer&lt; TransformType&gt;&nbsp;=20
transform;</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;typedef itk::ResampleImageFilter&lt;UCharImage, =
UCharImage&gt;=20
&nbsp;ResampleFilter;<BR>&nbsp;typedef UCharImage::SizeType&nbsp;=20
ImageSizeType;&nbsp;&nbsp;</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;typedef itk::LinearInterpolateImageFunction&lt;UCharImage, =
double&gt;=20
InterpolatorType;<BR>&nbsp;itk::SmartPointer&lt; =
InterpolatorType&gt;&nbsp;=20
interpolator;<BR>&nbsp;//////////////////////////////////////////////////=
/////////////////////////////////////////<BR>&nbsp;//<BR>&nbsp;//&nbsp;in=
itializations<BR>&nbsp;//<BR>&nbsp;//////////////////////////////////////=
/////////////////////////////////////////////////////</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;double result[12];<BR>&nbsp;memcpy(result, mat-&gt;Mat,=20
sizeof(double)*12);<BR>&nbsp;ConverterType converter =
;<BR>&nbsp;ConverterType=20
::ImagePointer fixedImage =3D=20
converter.GetImage(pFixed);&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbs=
p;=20
<BR>&nbsp;ConverterType ::ImagePointer movingImage =3D=20
converter.GetImage(pMoving);&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nb=
sp;=20
&nbsp;<BR>&nbsp; <BR>&nbsp;<BR>&nbsp;ResampleFilter::Pointer =
resampleFilter =3D=20
ResampleFilter::New();&nbsp;&nbsp;</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;transform =3D=20
TransformType::New();<BR>&nbsp;interpolator=3DInterpolatorType::New();<BR=
>&nbsp;TransformType::ParametersType=20
params(transform-&gt;GetNumberOfParameters()) ;<BR>&nbsp;<BR>&nbsp;for =
(int=20
pi=3D0;pi&lt;12;pi++)<BR>&nbsp;{<BR>&nbsp;&nbsp;&nbsp; params[pi] =3D=20
result[pi];<BR>&nbsp;&nbsp;&nbsp; <BR>&nbsp;}<BR>&nbsp;</DIV>
<DIV>&nbsp;</DIV>
<DIV><BR>&nbsp;UCharImage::Pointer outImage&nbsp; =3D=20
NULL;&nbsp;&nbsp;<BR>&nbsp;&nbsp;<BR>&nbsp;transform-&gt;SetParameters(=20
params);</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;&nbsp;&nbsp; resampleFilter-&gt;SetInput(movingImage=20
);&nbsp;&nbsp;&nbsp; <BR>&nbsp;resampleFilter=20
-&gt;SetTransform(transform.GetPointer());<BR>&nbsp;resampleFilter=20
-&gt;SetInterpolator(interpolator.GetPointer()); <BR>&nbsp;&nbsp;&nbsp;=20
resampleFilter-&gt;SetSize(=20
fixedImage-&gt;GetLargestPossibleRegion().GetSize());<BR>&nbsp;&nbsp;&nbs=
p;=20
resampleFilter-&gt;SetOutputOrigin( fixedImage-&gt;GetOrigin()=20
);<BR>&nbsp;&nbsp;&nbsp; resampleFilter-&gt;SetOutputSpacing(=20
fixedImage-&gt;GetSpacing() );<BR>&nbsp;&nbsp;&nbsp;=20
resampleFilter-&gt;SetDefaultPixelValue( 0 );</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;try<BR>&nbsp;{<BR>&nbsp;&nbsp;resampleFilter-&gt;Update();</DI=
V>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;}catch(itk::ExceptionObject &amp;Eo){</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;&nbsp;AfxMessageBox(Eo.GetDescription());<BR>&nbsp;}</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;outImage&nbsp;&nbsp; =3D =
resampleFilter-&gt;GetOutput();</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;ConverterType::PixelType* pMem =3D=20
outImage-&gt;GetBufferPointer();&nbsp;</DIV>
<DIV>&nbsp;</DIV>
<DIV><BR>&nbsp;if (pMem =3D=3D=20
NULL){<BR>&nbsp;&nbsp;AfxMessageBox("NUll");<BR>&nbsp;}&nbsp;<BR>&nbsp;</=
DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;memcpy(pResult-&gt;Mem, pMem, =
pFixed-&gt;bytesPerVolume);</DIV>
<DIV>&nbsp;</DIV>
<DIV><BR>&nbsp;return pResult; <BR>}</DIV>
<DIV>&nbsp;</DIV>
<DIV>&nbsp;</DIV>
<DIV>Regards,</DIV>
<DIV>CSPL</FONT></DIV></DIV>
<BLOCKQUOTE dir=3Dltr=20
style=3D"PADDING-RIGHT: 0px; PADDING-LEFT: 5px; MARGIN-LEFT: 5px; =
BORDER-LEFT: #000000 2px solid; MARGIN-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 title=3Dlng@insightful.com href=3D"mailto:lng@insightful.com">Lydia =
Ng</A>=20
  </DIV>
  <DIV style=3D"FONT: 10pt arial"><B>To:</B> <A =
title=3Daffable@hd2.dot.net.in=20
  href=3D"mailto:affable@hd2.dot.net.in">cspl</A> ; <A=20
  title=3Dinsight-users@public.kitware.com=20
  =
href=3D"mailto:insight-users@public.kitware.com">insight-users@public.kit=
ware.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 face=3DArial color=3D#0000ff size=3D2><SPAN=20
  class=3D494524517-11112002>CSPL,</SPAN></FONT></DIV>
  <DIV><FONT face=3DArial color=3D#0000ff size=3D2><SPAN=20
  class=3D494524517-11112002></SPAN></FONT>&nbsp;</DIV>
  <DIV><FONT face=3DArial color=3D#0000ff size=3D2><SPAN =
class=3D494524517-11112002>Your=20
  code snippet looks ok - you might need to post =
more</SPAN></FONT></DIV>
  <DIV><FONT face=3DArial color=3D#0000ff size=3D2><SPAN =
class=3D494524517-11112002>of=20
  your code to see if everything is setup correctly.</SPAN></FONT></DIV>
  <DIV><FONT face=3DArial color=3D#0000ff size=3D2><SPAN =
class=3D494524517-11112002>I=20
  wonder if it is a templating issue?</SPAN></FONT></DIV>
  <DIV><FONT face=3DArial color=3D#0000ff size=3D2><SPAN =
class=3D494524517-11112002>Are=20
  you using the ITK registration components? </SPAN></FONT></DIV>
  <DIV><FONT face=3DArial color=3D#0000ff size=3D2><SPAN =
class=3D494524517-11112002>Have=20
  you set the transform and interpolation template =
parameter</SPAN></FONT></DIV>
  <DIV><FONT face=3DArial color=3D#0000ff size=3D2><SPAN =
class=3D494524517-11112002>of=20
  the ResampleImageFilter to the appropriate =
classes?</SPAN></FONT></DIV>
  <DIV><FONT face=3DArial color=3D#0000ff size=3D2><SPAN=20
  class=3D494524517-11112002></SPAN></FONT>&nbsp;</DIV>
  <DIV><FONT face=3DArial color=3D#0000ff size=3D2><SPAN=20
  class=3D494524517-11112002>-Lydia</SPAN></FONT></DIV></DIV>
  <BLOCKQUOTE dir=3Dltr style=3D"MARGIN-RIGHT: 0px">
    <DIV class=3DOutlookMessageHeader dir=3Dltr align=3Dleft><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>&nbsp;I have some doubts =
regarding Resample=20
    image filter.After </FONT></DIV>
    <DIV>&nbsp;</DIV>
    <DIV><FONT face=3DArial size=3D2>registring fixed and moving =
volumes, I applied=20
    registration output matrix on </FONT></DIV>
    <DIV>&nbsp;</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>&nbsp;</DIV>
    <DIV><FONT face=3DArial size=3D2>have applied the same matrix on =
same moving=20
    volume.Output volume is </FONT></DIV>
    <DIV>&nbsp;</DIV>
    <DIV><FONT face=3DArial size=3D2>totally different from resample =
output.Can u=20
    please tell what component </FONT></DIV>
    <DIV>&nbsp;</DIV>
    <DIV><FONT face=3DArial size=3D2>is missing in resampling.I am =
enclosing the=20
    resample code.&nbsp; </FONT></DIV>
    <DIV>&nbsp;</DIV>
    <DIV><FONT face=3DArial size=3D2>Please click on this url to see the =
pictures=20
    </FONT></DIV>
    <DIV>&nbsp;</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>&nbsp;</DIV>
    <DIV><FONT face=3DArial size=3D2>//code<BR>ResampleFilter::Pointer=20
    resampleFilter =3D ResampleFilter::New();&nbsp;&nbsp;</FONT></DIV>
    <DIV>&nbsp;</DIV>
    <DIV><FONT face=3DArial size=3D2>&nbsp;transform =3D=20
    =
TransformType::New();<BR>&nbsp;interpolator=3DInterpolatorType::New();<BR=
>&nbsp;TransformType::ParametersType=20
    </FONT></DIV>
    <DIV>&nbsp;</DIV>
    <DIV><FONT face=3DArial =
size=3D2>params(transform-&gt;GetNumberOfParameters())=20
    ;<BR>&nbsp;<BR>&nbsp;for (int=20
    pi=3D0;pi&lt;12;pi++)<BR>&nbsp;{<BR>&nbsp;&nbsp;&nbsp; params[pi] =
=3D=20
    result[pi];<BR>&nbsp;&nbsp;&nbsp; <BR>&nbsp;}<BR>&nbsp;</FONT></DIV>
    <DIV>&nbsp;</DIV><FONT face=3DArial size=3D2>
    <DIV><BR>&nbsp;UCharImage::Pointer outImage&nbsp; =3D=20
    =
NULL;&nbsp;&nbsp;<BR>&nbsp;&nbsp;<BR>&nbsp;transform-&gt;SetParameters(=20
    params);</DIV>
    <DIV>&nbsp;</DIV>
    <DIV><BR>&nbsp;&nbsp;&nbsp; resampleFilter-&gt;SetInput(movingImage=20
    );&nbsp;&nbsp;&nbsp; <BR>&nbsp;resampleFilter=20
    -&gt;SetTransform(transform.GetPointer());<BR>&nbsp;resampleFilter=20
    -&gt;SetInterpolator(interpolator.GetPointer()); =
<BR>&nbsp;&nbsp;&nbsp;=20
    resampleFilter-&gt;SetSize( </DIV>
    <DIV>&nbsp;</DIV>
    =
<DIV>fixedImage-&gt;GetLargestPossibleRegion().GetSize());<BR>&nbsp;&nbsp=
;&nbsp;=20
    resampleFilter-&gt;SetOutputOrigin( fixedImage-&gt;GetOrigin()=20
    );<BR>&nbsp;&nbsp;&nbsp; resampleFilter-&gt;SetOutputSpacing(=20
    fixedImage-&gt;GetSpacing() );</DIV>
    <DIV>&nbsp;</DIV>
    =
<DIV><BR>&nbsp;try<BR>&nbsp;{<BR>&nbsp;&nbsp;resampleFilter-&gt;Update();=
</DIV>
    <DIV>&nbsp;</DIV>
    <DIV>&nbsp;}catch(itk::ExceptionObject &amp;Eo){</DIV>
    <DIV>&nbsp;</DIV>
    =
<DIV>&nbsp;&nbsp;AfxMessageBox(Eo.GetDescription());<BR>&nbsp;}</DIV>
    <DIV>&nbsp;</DIV>
    <DIV>&nbsp;outImage&nbsp;&nbsp; =3D =
resampleFilter-&gt;GetOutput();</DIV>
    <DIV>&nbsp;</DIV>
    <DIV><BR>Regards,<BR>CSPL=20
<BR></FONT></DIV></DIV></BLOCKQUOTE></BLOCKQUOTE></BODY></HTML>

------=_NextPart_000_0088_01C28A42.81582F70--