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

[sidharta]

Dear Dženan,


Thank you for your swift reply.

This is the entire function which I used:

int ResampleImage(std::vector<Vector3d> trajectory_points,
std::vector<Vector3d> coordinate_system, InputImageType::Pointer image){

	Vector3d traj_start, traj_end, third_point;
	traj_start = trajectory_points.at(0);
	traj_end = trajectory_points.back();

	third_point = traj_end + coordinate_system.at(3) * 2 /
image->GetSpacing()[0];

	Vector3d ST, SP, normal_to_plane;
	ST = traj_end - traj_start;
	SP = third_point - traj_start;
	normal_to_plane = ST.cross(SP);

	InputImageType::DirectionType imageDirection, trajectoryDirection;
	imageDirection = image->GetDirection();
	std::cout << "Image direction cosines : \n" << imageDirection << std::endl;

	trajectoryDirection[0][0] = ST.normalized()[0];
	trajectoryDirection[0][1] = ST.normalized()[1];
	trajectoryDirection[0][2] = ST.normalized()[2];

	Vector3d second = normal_to_plane.cross(ST);
	Vector3d normalized_second = second.normalized();
	trajectoryDirection[1][0] = normalized_second[0];
	trajectoryDirection[1][1] = normalized_second[1];
	trajectoryDirection[1][2] = normalized_second[2];

	Vector3d normalized_normal = normal_to_plane.normalized();
	trajectoryDirection[2][0] = normalized_normal[0];
	trajectoryDirection[2][1] = normalized_normal[1];
	trajectoryDirection[2][2] = normalized_normal[2];

	std::cout << "Trajectory direction cosines : \n" << trajectoryDirection <<
std::endl;

	InputImageType::IndexType voxelIndex;
	InputImageType::PointType point;

	point[0] = traj_start[0], point[1] = traj_start[1], point[2] =
traj_start[2];
	image->TransformPhysicalPointToIndex(point, voxelIndex);
	std::cout << "Index for point " << point << " is " << voxelIndex <<
std::endl;

	point[0] = traj_end[0], point[1] = traj_end[1], point[2] = traj_end[2];
	image->TransformPhysicalPointToIndex(point, voxelIndex);
	std::cout << "Index for point " << point << " is " << voxelIndex <<
std::endl;

	double mpr_start[3] = { traj_start[0], traj_start[1], traj_start[2] };

	typedef itk::ResampleImageFilter<InputImageType, InputImageType>
resampleFilterType;
	resampleFilterType::Pointer resampleFilter = resampleFilterType::New();

	typedef itk::LinearInterpolateImageFunction<InputImageType, InputPixelType>
interpolatorType;
	interpolatorType::Pointer linearInterpolator = interpolatorType::New();

	resampleFilter->SetInterpolator(linearInterpolator);
	resampleFilter->SetDefaultPixelValue(-2000);
	resampleFilter->SetOutputSpacing(image->GetSpacing());
	resampleFilter->SetOutputOrigin(mpr_start);
	resampleFilter->SetOutputDirection(trajectoryDirection);
	InputImageType::SizeType size;
	size[0] = trajectory_points.size();
	size[1] = 1.0 / image->GetSpacing()[1];
	size[2] = 1.0 / image->GetSpacing()[2];
	resampleFilter->SetSize(size); // This is for the trajectory, this
corresponds to a size of [137, 8, 8] which is essentially the volume I want
to extract.
	//resampleFilter->SetSize(image->GetLargestPossibleRegion().GetSize()); //
This is for the entire resampled image
	resampleFilter->SetInput(image);
	//resampleFilter->Update();

	//InputImageType::Pointer resampled_image = InputImageType::New();
	//resampled_image = resampleFilter->GetOutput();

	FileWriter(resampleFilter->GetOutput(),
"C:/Users/api/Desktop/resampled_volume_1.mhd");
	return EXIT_SUCCESS;

}

I should also add that when I open the original image with the resampled
image in ITK-SNAP (by adding another image), the overlay is exactly the
same. Additionally, when I opened this in Slicer, only the top part of the
volume seems to be rotated and the rest isn't. Kindly let me know if you
need more information. 
I apologize for asking too many questions. We don't have an experienced ITK
user in our lab, although vector understanding is sound :)

Thank you for your time. 



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