[ITK-users] [ITK] Extract slice from ResampleFilter

[sidharta]

Hi Dženan,

So I am doing something as follows: 

// Vectors for calculating normal_to_plane - The plane/slice I want to
extract

Vector3d ST, SP, normal_to_plane; // This is a custom 3D Vector library
ST = traj_start - traj_end; // start and end of the trajectory - vector 1
SP = traj_start - third_point; // start and third point gives another vector
- vector 2
// Now with two vectors I can get the direction of the normal to the plane
normal_to_plane = ST.cross(SP); // cross product essentially

// Getting angles between normal to the plane and the coordinate axes
double angle_plane_z = std::acos((normal_to_plane[2] * 1) /
norm(normal_to_plane)); // angle between plane and z axis
double angle_plane_y = std::acos((normal_to_plane[1] * 1) /
norm(normal_to_plane)); // angle between plane and y axis
double angle_plane_x = std::acos((normal_to_plane[0] * 1) /
norm(normal_to_plane)); // angle between plane and x axis

// Instantiate an Euler Transform
typedef itk::Euler3DTransform<double> ETransformType;
ETransformType::Pointer eulerTransform = ETransformType::New();

// This because - Quote from ITKSoftwareGuide "Rotations are performed
around the origin of physical coordinates—not the image origin nor the image
center"
ETransformType::OutputVectorType translation1;
translation1[0] = -origin_dcm[0]; // origin_dcm is the origin of the dicom
translation1[1] = -origin_dcm[1]; 
translation1[2] = -origin_dcm[2];
eulerTransform->Translate(translation1);
eulerTransform->SetRotation(angle_plane_x, angle_plane_y, angle_plane_z);

// Get back to the origin of the dicom
ETransformType::OutputVectorType translation2;
translation2[0] = origin_dcm[0];
translation2[1] = origin_dcm[1];
translation2[2] = origin_dcm[2];

// Instantiate the ResampleFilter
// typedef double InternalPixelType;
// typedef itk::Image< InternalPixelType, Dimension > InternalImageType; //
For clarity
typedef itk::ResampleImageFilter<InternalImageType, InternalImageType>
ResampleFilterType;

ResampleFilterType::Pointer filter = ResampleFilterType::New();

typedef itk::NearestNeighborInterpolateImageFunction<
		InternalImageType, double> InterpolatorType;
InterpolatorType::Pointer interpolator = InterpolatorType::New();
filter->SetInterpolator(interpolator);

filter->SetDefaultPixelValue(0);
filter->SetOutputSpacing(image->GetSpacing());
filter->SetOutputOrigin(image->GetOrigin()); // If I change this to
traj_start, the dicom looks completely off, as in not intended.
filter->SetOutputDirection(image->GetDirection()); // I have my doubts about
this too
filter->SetSize(image->GetSize()); // Of course, as you suggested, I would
set the thickness here as 1, but I was doing this for testing purposes.
filter->SetInput(image);

So, what do you think about this code? Am I doing it correctly? I feel
really dumb asking you to check my code but unfortunately I can't find any
relevant example for 3D resampling. Additionally, (no offense) but I don't
get what you mean by "target region index" with respect to the resampling
filter? Do you mean the SetOrigin()? As you also see, I am essentially not
applying any translation to the origin. Just a rotation for testing
purposes.

Thank you for your time.
Best Regards,
Sidharta





--
View this message in context: http://itk-users.7.n7.nabble.com/Extract-slice-from-ResampleFilter-tp38074p38110.html
Sent from the ITK - Users mailing list archive at Nabble.com.