Discussion:
[ITK-users] Extract slice from ResampleFilter
sidharta
2017-04-06 12:34:14 UTC
Permalink
Dear all,

I am trying to extract density profile along a trajectory in CT. Using some
adjustments to contour-based segmentation and getting a mesh from the binary
mask, I computed certain trajectories tangential to certain faces/vertices
of the mesh.

Now, I have the coordinates of the target point and a tangent to the target.
So, depending on the spacing, I extract all the points along the tangent
which are 5mm away from the target. I am using following code to generate
densities after I have got the physical coordinates along this line.

void generateDensities(std::vector<Vector3d> points,
InternalImageType::Pointer image){

InternalImageType::IndexType pixelIndex;
PointType point;
itk::ContinuousIndex<double, Dimension> pixel;
linearInterpolationContinuousIndexType::Pointer interpolate =
linearInterpolationContinuousIndexType::New();

interpolate->SetInputImage(image);

std::vector<int> HU_values;
std::vector<double> HU_values_interpolate;
std::vector<InternalImageType::IndexType> Pixel_indices;
std::vector<itk::ContinuousIndex&lt;double, Dimension>>
Pixel_indices_interpolated;
// std::cout << "Generating with interpolation" << std::endl;
for (int i = 0; i < points.size(); i++)
{
Vector3d temp_point = points.at(i);
point[0] = temp_point[0];
point[1] = temp_point[1];
point[2] = temp_point[2];
// std::cout << "Getting continuous index " << std::endl;
bool isInside = image->TransformPhysicalPointToContinuousIndex(point,
pixel);
if (isInside)
{
InternalImageType::PixelType pixelValue =
interpolate->EvaluateAtContinuousIndex(pixel);
HU_values_interpolate.push_back(pixelValue);
Pixel_indices_interpolated.push_back(pixel);
}
else
std::cout << "Point " << point << " with index " << pixel << "is outside"
<< std::endl;
}
All_HU_values_interpolated.push_back(HU_values_interpolate);
All_pixel_indices_interpolated.push_back(Pixel_indices_interpolated);
HU_values_interpolate.clear();
HU_values.clear();
}

Now, for verification that I am actually getting in fact the correct values,
I would like to extract the oblique slice along this trajectory and I figure
I have to use the Resample filter to "transform" the CT so that rotation
puts the resampled CT coordinate system parallel to the trajectory and the
origin translated to the start of the trajectory. Now, along the trajectory,
with the set of collinear points (along the trajectory), I can get a lot of
planes (oblique slices essentially) which have the trajectory in them. So to
get around this, I figure if I give another point which is non-collinear to
these points, I can specify the exact oblique slice I want. However, this is
something I have to figure out for a lot of trajectories for the same CT and
then correspondingly for others. Can someone tell me how can I achieve this
"automatically" for all trajectories/CTs?
PS. All CTs have the same spacing and origin.
Kindly also mention if you think the function for generating the densities
is correct or not? Eitherways I would like to validate.

Thank you,
Sidharta




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

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/insight-users
Dženan Zukić
2017-04-08 20:51:00 UTC
Permalink
Hi Sidharta,

your function to calculate interpolated Hounsfield values seems correct.

Your problem of selecting oblique slice to resample can be reduced to
finding a plane which passed through your first and last point, and as an
arbitrary constraint passes through point [0,0,0]. If this point is
collinear with your line then pick point [1,0,0], [0,1,0] or [0,0,1]. One
of them will not be collinear.

In order to resample your oblique slice, you need to figure out a transform
matrix which transforms your first point to coordinates [0,0,0], your last
point to coordinates e.g. [1000,0,0] and transforms the normal of your
plane to be parallel with vector [0,0,1].

Then use resample filter, set your transform matrix and set as target
image's region index [0,0,0] and size [1001,1,1]. Alternatively index
[0,-10,0] and size [1001,21,1] if you want to have a real 2D slice to look
at.

Hope this helps,
DÅŸenan
Post by sidharta
Dear all,
I am trying to extract density profile along a trajectory in CT. Using some
adjustments to contour-based segmentation and getting a mesh from the binary
mask, I computed certain trajectories tangential to certain faces/vertices
of the mesh.
Now, I have the coordinates of the target point and a tangent to the target.
So, depending on the spacing, I extract all the points along the tangent
which are 5mm away from the target. I am using following code to generate
densities after I have got the physical coordinates along this line.
void generateDensities(std::vector<Vector3d> points,
InternalImageType::Pointer image){
InternalImageType::IndexType pixelIndex;
PointType point;
itk::ContinuousIndex<double, Dimension> pixel;
linearInterpolationContinuousIndexType::Pointer interpolate =
linearInterpolationContinuousIndexType::New();
interpolate->SetInputImage(image);
std::vector<int> HU_values;
std::vector<double> HU_values_interpolate;
std::vector<InternalImageType::IndexType> Pixel_indices;
std::vector<itk::ContinuousIndex&lt;double, Dimension>>
Pixel_indices_interpolated;
// std::cout << "Generating with interpolation" << std::endl;
for (int i = 0; i < points.size(); i++)
{
Vector3d temp_point = points.at(i);
point[0] = temp_point[0];
point[1] = temp_point[1];
point[2] = temp_point[2];
// std::cout << "Getting continuous index " << std::endl;
bool isInside = image->TransformPhysicalPointToContin
uousIndex(point,
pixel);
if (isInside)
{
InternalImageType::PixelType pixelValue =
interpolate->EvaluateAtContinuousIndex(pixel);
HU_values_interpolate.push_back(pixelValue);
Pixel_indices_interpolated.push_back(pixel);
}
else
std::cout << "Point " << point << " with index "
<< pixel << "is outside"
<< std::endl;
}
All_HU_values_interpolated.push_back(HU_values_interpolate);
All_pixel_indices_interpolated.push_back(Pixel_
indices_interpolated);
HU_values_interpolate.clear();
HU_values.clear();
}
Now, for verification that I am actually getting in fact the correct values,
I would like to extract the oblique slice along this trajectory and I figure
I have to use the Resample filter to "transform" the CT so that rotation
puts the resampled CT coordinate system parallel to the trajectory and the
origin translated to the start of the trajectory. Now, along the trajectory,
with the set of collinear points (along the trajectory), I can get a lot of
planes (oblique slices essentially) which have the trajectory in them. So to
get around this, I figure if I give another point which is non-collinear to
these points, I can specify the exact oblique slice I want. However, this is
something I have to figure out for a lot of trajectories for the same CT and
then correspondingly for others. Can someone tell me how can I achieve this
"automatically" for all trajectories/CTs?
PS. All CTs have the same spacing and origin.
Kindly also mention if you think the function for generating the densities
is correct or not? Eitherways I would like to validate.
Thank you,
Sidharta
--
View this message in context: http://itk-users.7.n7.nabble.
com/Extract-slice-from-ResampleFilter-tp38074.html
Sent from the ITK - Users mailing list archive at Nabble.com.
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users
sidharta
2017-04-10 09:29:55 UTC
Permalink
Dear Dženan,

Thank you for your reply. Actually, I found a third point which gives me the
exact (informative) slice I am looking for. Should I resample the image
twice, first rotate the image to the direction of the trajectory and move
the origin to the start of the trajectory and second rotate the image to the
angle between the trajectory and the third point? I will try what you said
and get back to you. I am asking because I don't think resampling twice
makes sense.

Thank you for your time.
Sidharta



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

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/list
Dženan Zukić
2017-04-10 13:24:50 UTC
Permalink
​Hi Sidharta,

you definitely should not resample twice. All you need is to compound the
two transforms (multiply their matrices)​. And apply a single resampling.
But for debugging purposes, you could also display the resampled output
after only the first transform.

Regards,
DÅŸenan
Dear DÅŸenan,
Thank you for your reply. Actually, I found a third point which gives me the
exact (informative) slice I am looking for. Should I resample the image
twice, first rotate the image to the direction of the trajectory and move
the origin to the start of the trajectory and second rotate the image to the
angle between the trajectory and the third point? I will try what you said
and get back to you. I am asking because I don't think resampling twice
makes sense.
Thank you for your time.
Sidharta
--
View this message in context: http://itk-users.7.n7.nabble.
com/Extract-slice-from-ResampleFilter-tp38074p38084.html
Sent from the ITK - Users mailing list archive at Nabble.com.
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users
_______________________________________________
Community mailing list
http://public.kitware.com/mailman/listinfo/community
sidharta
2017-04-11 16:00:13 UTC
Permalink
Hey Dženan,

So, I kind of got the equation of the plane and the normal to the plane
where all the three non-collinear points lie. Then, I got the angle between
this normal and the z-axis of the image. I resample the image by rotating in
3D about the z-axis. How do I extract the exact slice from the resampled
image? Is it actually required to resample the 3D Dicom, can I not resample
a 2D image (just for optimization purposes)? I read the documentation
(ITKSoftwareGuide) which shows one example for 3D resampling, however not
how to save a particular slice. Kindly correct me if I am wrong.


Thank you for your time.
Best Regards,
Sidharta



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

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http://public.k
Dženan Zukić
2017-04-11 18:03:19 UTC
Permalink
Hi Sidartha,

you can use extract filter
<https://itk.org/Doxygen/html/classitk_1_1ExtractImageFilter.html> to
extract a single slice from a 3D image. Take a look at the examples liked
from the documentation page.

Regards,
DÅŸenan
Hey DÅŸenan,
So, I kind of got the equation of the plane and the normal to the plane
where all the three non-collinear points lie. Then, I got the angle between
this normal and the z-axis of the image. I resample the image by rotating in
3D about the z-axis. How do I extract the exact slice from the resampled
image? Is it actually required to resample the 3D Dicom, can I not resample
a 2D image (just for optimization purposes)? I read the documentation
(ITKSoftwareGuide) which shows one example for 3D resampling, however not
how to save a particular slice. Kindly correct me if I am wrong.
Thank you for your time.
Best Regards,
Sidharta
--
View this message in context: http://itk-users.7.n7.nabble.
com/Extract-slice-from-ResampleFilter-tp38074p38091.html
Sent from the ITK - Users mailing list archive at Nabble.com.
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users
sidharta
2017-04-12 11:44:55 UTC
Permalink
Hi Dženan,

Of course I would use that (sorry for not mentioning it that I know this
already), but how do I get this particular plane (2D slice)? I don't want to
hard code anything.

Hope you understand
Thank you for your time
Sidharta



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

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http://public.kitware.com/
Francois Budin
2017-04-12 12:46:12 UTC
Permalink
Hello Sidharta,

You can use the filter called ExtractImageFilter [1]. You should be able to
achieve your goal with that.
Hope this helps,

Francois

[1] https://itk.org/Doxygen/html/classitk_1_1ExtractImageFilter.html
Hi DÅŸenan,
Of course I would use that (sorry for not mentioning it that I know this
already), but how do I get this particular plane (2D slice)? I don't want to
hard code anything.
Hope you understand
Thank you for your time
Sidharta
--
View this message in context: http://itk-users.7.n7.nabble.
com/Extract-slice-from-ResampleFilter-tp38074p38097.html
Sent from the ITK - Users mailing list archive at Nabble.com.
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users
Dženan Zukić
2017-04-12 13:03:52 UTC
Permalink
Hi Sidharta,

to get the particular plane you want, set up the transform in such a way
that that slice gets transformed to some know spatial position (e.g. your
first point to coordinates [0,0,0] and your last point to coordinates
[1000,0,0]). Then you only need to resample your volume to a 3D image of
thickness 1 slice which passes through points [0,0,0] and [1000,0,0]. You
don't even need extract filter with this approach.

Regards,
DÅŸenan Zukić, PhD, Senior R&D Engineer, Kitware (Carrboro, N.C.)
Hi DÅŸenan,
Of course I would use that (sorry for not mentioning it that I know this
already), but how do I get this particular plane (2D slice)? I don't want to
hard code anything.
Hope you understand
Thank you for your time
Sidharta
--
View this message in context: http://itk-users.7.n7.nabble.
com/Extract-slice-from-ResampleFilter-tp38074p38097.html
Sent from the ITK - Users mailing list archive at Nabble.com.
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users
sidharta
2017-04-13 13:04:04 UTC
Permalink
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.
_____________________________________
Powered by www.kitware.com

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http
Dženan Zukić
2017-04-13 14:37:12 UTC
Permalink
Hi Sidartha,

the code looks OK. The index is [0,0,0] by default, so just setting size is
enough. My idea was not to apply translation2, then you would not have to
set origin on the resample filter (it would start from [0,0,0]).

I think you are forgetting to apply second translation, and you are not
applying any transformation to the resample filter. Something like this is
needed:
filter->SetTransform(eulerTransform);

What results are you getting?

Regards,
DÅŸenan Zukić, PhD, Senior R&D Engineer, Kitware (Carrboro, N.C.)
Hi DÅŸenan,
// 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.
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users
sidharta
2017-04-13 15:01:24 UTC
Permalink
Hi Dženan,

Sorry about that, I missed the line in the post. I am, in fact, adding that
line filter->SetTransform().

Apologies, but I don't know how I can share .dcm/.mhd image. I am just going
to add the desirable which I did using Amira, which is where I got the idea
of using EulerTransform.
<Loading Image...> . The big sphere
is the third point and the two small ones are the start and end of the
trajectory. The slices on the bounding box are essentially on the origin.
So, with vector math I could get the equation of the plane and so forth. If
you want to see the output .dcm, kindly let me know how I can share it with
you ( from a cloud probably but I don't know what you prefer).

Thank you for your swift reply and patience.

Regards,
Sidharta



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

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailma
Dženan Zukić
2017-04-13 15:16:23 UTC
Permalink
​Hi Siddartha,

you can share using cloud tools, e.g. Dropbox or Google Drive.

But I think you forgot to ask the question: I don't know what you are
asking.

Regards,
DÅŸenan Zukić, PhD, Senior R&D Engineer, Kitware (Carrboro, N.C.)​
Hi DÅŸenan,
Sorry about that, I missed the line in the post. I am, in fact, adding that
line filter->SetTransform().
Apologies, but I don't know how I can share .dcm/.mhd image. I am just going
to add the desirable which I did using Amira, which is where I got the idea
of using EulerTransform.
<http://itk-users.7.n7.nabble.com/file/n38115/snapshot.jpg> . The big sphere
is the third point and the two small ones are the start and end of the
trajectory. The slices on the bounding box are essentially on the origin.
So, with vector math I could get the equation of the plane and so forth. If
you want to see the output .dcm, kindly let me know how I can share it with
you ( from a cloud probably but I don't know what you prefer).
Thank you for your swift reply and patience.
Regards,
Sidharta
--
View this message in context: http://itk-users.7.n7.nabble.
com/Extract-slice-from-ResampleFilter-tp38074p38115.html
Sent from the ITK - Users mailing list archive at Nabble.com.
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users
sidharta
2017-04-13 15:33:52 UTC
Permalink
So if you think the code is correct and is doing what you suggested (or what
I think it should do using the eulerTransform), you mentioned setting
targetRegionIndex in ResampleFilter to get one single slice. I don't see
this targetRegionIndex or do you mean set origin. I am going to send you the
link to the dropbox folder in a private mail.



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

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/insight-users
Dženan Zukić
2017-04-13 17:44:58 UTC
Permalink
The target region is the region (index and size) you set in the resample
filter. It is just a difference in terminology between you and me.

I took a look at the images you supplied. Because of transformations, this
is hard to debug.

I realized it would be easier to construct the output image to be in the
same physical space as the input image, and have identity transform applied
to resample filter (used by default). Pseudo-code:

ImageType::DirectionType dir;
ST.normalize();
dir[0]=ST; //vector from first point to second
vec3 second=normal.cross(ST); //make the second vector orthogonal to first
and the normal
second.normalize();
dir[1]=second;
normal.normalize();
dir[2]=normal;
filter->SetOutputDirection(dir);
filter->SetOrigin(traj_start);
//it will be more math to compute proper size (so that resampled image ends
at traj_end), and proper spacing (project original spacing onto coordinate
system established by new dir matrix or use isotropic spacing equal to
cubic root of sp[0]*sp[1]*sp[2]).

With this approach you can use e.g. Slicer <https://www.slicer.org/> to
look at both the original and resampled images, and also define fiducial
points for your 3 points of interest.

P.S. If you share any more images, please save them with compression
applied (writer->SetUseCompression(true);).

Regards,
DÅŸenan Zukić, PhD, Senior R&D Engineer, Kitware (Carrboro, N.C.)
Post by sidharta
So if you think the code is correct and is doing what you suggested (or what
I think it should do using the eulerTransform), you mentioned setting
targetRegionIndex in ResampleFilter to get one single slice. I don't see
this targetRegionIndex or do you mean set origin. I am going to send you the
link to the dropbox folder in a private mail.
--
View this message in context: http://itk-users.7.n7.nabble.
com/Extract-slice-from-ResampleFilter-tp38074p38118.html
Sent from the ITK - Users mailing list archive at Nabble.com.
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users
sidharta
2017-04-19 13:40:38 UTC
Permalink
Hi Dženan,

Thank you for your reply. I implemented this, and learned some cool stuff in
3DSlicer (in order to compare the two images) but as you said, the size is
difficult to compute. When I set the origin explicitly, to the start of the
trajectory, the resampled image looks good in the direction I want it to be
in, however, lot of information is lost. After putting my head into it for
the past two weeks, I realized there are other ways to achieve I want this
too. Does it make sense to you that I construct another 2D image (or 3D with
thickness 1) by essentially copying the voxel values at the trajectory
points and around it (probably the size equal to the dimensions of the
original image in the x,y direction)? The aim ultimately is in fact to get
one slice which best represents the HU values. I apologize for asking such
silly questions.

Just putting the code of what I did, let me know if it is in fact what you
suggested me to do:

Vector3d ST, SP, normal_to_plane;
ST = traj_end - traj_start;
SP = third_point - traj_start;
normal_to_plane = ST.cross(SP);
std::cout << "Vector ST : " << ST << " Vector SP : " << SP << " Direction
of normal : " << normal_to_plane << std::endl;

InternalImageType::DirectionType dir;
// vector from start to the end of the trajectory
Vector3d normalized_ST = ST.normalized();
dir[0][0] = normalized_ST[0];
dir[0][1] = normalized_ST[1];
dir[0][2] = normalized_ST[2];

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

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

// traj_start += normalized_second * 10;
double mpr_start[3] = { traj_start[0] , traj_start[1] , traj_start[2]};

std::cout << "Previous image direction " << image->GetDirection() <<
std::endl;
std::cout << "New Image Direction " << dir << std::endl;

getSliceMPR(image, mpr_start, dir);

void getSliceMPR(InternalImageType::Pointer &image, double
output_origin[Dimension], InternalImageType::DirectionType dir)
{

// Instantiate Resample Filter here
ResampleFilterType::Pointer filter = ResampleFilterType::New();
// Instantiate Transform for resampling
// Changing the Output direction instead of setting the transform
// filter->SetTransform(transform);

typedef itk::NearestNeighborInterpolateImageFunction<
InternalImageType, double> InterpolatorType;

// Instantiate the interpolator
InterpolatorType::Pointer interpolator = InterpolatorType::New();
filter->SetInterpolator(interpolator);

// Set pixels of the output image that end up being mapped outside the
extent of the input image to -1500
filter->SetDefaultPixelValue(-1500);

// Sampling grid of the output space - pixel spacing in mm
// Hard coded here
filter->SetOutputSpacing(image->GetSpacing());


// filter->SetOutputOrigin(image->GetOrigin());

filter->SetOutputOrigin(output_origin);
filter->SetOutputDirection(dir);
InternalImageType::SizeType size;
filter->SetSize(image->GetLargestPossibleRegion().GetSize());

filter->SetInput(image);

typedef itk::ImageFileWriter<InternalImageType> WriterType;

WriterType::Pointer writer = WriterType::New();
writer->SetFileName("C:/Users/api/Desktop/Resampled/resampled_2.mhd");
writer->SetUseCompression(true);
writer->SetInput(filter->GetOutput());
try{
writer->Update();

}
catch (itk::ExceptionObject & excep){

std::cerr << "Exception catched!.." << std::endl;
std::cerr << excep << std::endl;
}

}



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

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http://public.kitware
Dženan Zukić
2017-04-19 15:29:53 UTC
Permalink
Hi Sidharta,

you seem to have implemented what I suggested.

"essentially copying the voxel values at the trajectory points and around
it" can be achieved by using nearest neighbor interpolator instead of
linear interpolator. But it looks like you are already using nearest
neighbor, so you should try using linear interpolator, or even windowed sinc
<https://itk.org/Doxygen/html/classitk_1_1WindowedSincInterpolateImageFunction.html>
interpolator to get higher quality.

Regards,
DÅŸenan Zukić, PhD, Senior R&D Engineer, Kitware (Carrboro, N.C.)
Hi DÅŸenan,
Thank you for your reply. I implemented this, and learned some cool stuff in
3DSlicer (in order to compare the two images) but as you said, the size is
difficult to compute. When I set the origin explicitly, to the start of the
trajectory, the resampled image looks good in the direction I want it to be
in, however, lot of information is lost. After putting my head into it for
the past two weeks, I realized there are other ways to achieve I want this
too. Does it make sense to you that I construct another 2D image (or 3D with
thickness 1) by essentially copying the voxel values at the trajectory
points and around it (probably the size equal to the dimensions of the
original image in the x,y direction)? The aim ultimately is in fact to get
one slice which best represents the HU values. I apologize for asking such
silly questions.
Just putting the code of what I did, let me know if it is in fact what you
Vector3d ST, SP, normal_to_plane;
ST = traj_end - traj_start;
SP = third_point - traj_start;
normal_to_plane = ST.cross(SP);
std::cout << "Vector ST : " << ST << " Vector SP : " << SP << " Direction
of normal : " << normal_to_plane << std::endl;
InternalImageType::DirectionType dir;
// vector from start to the end of the trajectory
Vector3d normalized_ST = ST.normalized();
dir[0][0] = normalized_ST[0];
dir[0][1] = normalized_ST[1];
dir[0][2] = normalized_ST[2];
Vector3d second = normal_to_plane.cross(ST);
Vector3d normalized_second = second.normalized();
dir[1][0] = normalized_second[0];
dir[1][1] = normalized_second[1];
dir[1][2] = normalized_second[2];
Vector3d normalized_normal = normal_to_plane.normalized();
dir[2][0] = normalized_normal[0];
dir[2][1] = normalized_normal[1];
dir[2][2] = normalized_normal[2];
// traj_start += normalized_second * 10;
double mpr_start[3] = { traj_start[0] , traj_start[1] , traj_start[2]};
std::cout << "Previous image direction " << image->GetDirection() <<
std::endl;
std::cout << "New Image Direction " << dir << std::endl;
getSliceMPR(image, mpr_start, dir);
void getSliceMPR(InternalImageType::Pointer &image, double
output_origin[Dimension], InternalImageType::DirectionType dir)
{
// Instantiate Resample Filter here
ResampleFilterType::Pointer filter = ResampleFilterType::New();
// Instantiate Transform for resampling
// Changing the Output direction instead of setting the transform
// filter->SetTransform(transform);
typedef itk::NearestNeighborInterpolateImageFunction<
InternalImageType, double> InterpolatorType;
// Instantiate the interpolator
InterpolatorType::Pointer interpolator = InterpolatorType::New();
filter->SetInterpolator(interpolator);
// Set pixels of the output image that end up being mapped outside the
extent of the input image to -1500
filter->SetDefaultPixelValue(-1500);
// Sampling grid of the output space - pixel spacing in mm
// Hard coded here
filter->SetOutputSpacing(image->GetSpacing());
// filter->SetOutputOrigin(image->GetOrigin());
filter->SetOutputOrigin(output_origin);
filter->SetOutputDirection(dir);
InternalImageType::SizeType size;
filter->SetSize(image->GetLargestPossibleRegion().GetSize());
filter->SetInput(image);
typedef itk::ImageFileWriter<InternalImageType> WriterType;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName("C:/Users/api/Desktop/Resampled/
resampled_2.mhd");
writer->SetUseCompression(true);
writer->SetInput(filter->GetOutput());
try{
writer->Update();
}
catch (itk::ExceptionObject & excep){
std::cerr << "Exception catched!.." << std::endl;
std::cerr << excep << std::endl;
}
}
--
View this message in context: http://itk-users.7.n7.nabble.
com/Extract-slice-from-ResampleFilter-tp38074p38129.html
Sent from the ITK - Users mailing list archive at Nabble.com.
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users
sidharta
2017-05-11 15:20:49 UTC
Permalink
Dear Dženan,

First of all, thank you for all your help.

So, from the previous posts I was able to resample the CBCT image so that it
lies in the direction of the Trajectory. So far so good. Now, I am trying to
extract only a part of the volume in this resampled image. So, I figure that
since I already calculated the direction, I could just set the size of the
output image to be what I want in the "ResampleImageFilter" and it should
work fine. However, when I do this, the extracted volume is in a different
direction than I assumed it to be. For example,

<Loading Image...>

In the above image, I assumed the part of the volume that I extract is
actually in the direction shown in the "Green" line, however I instead get
the volume along the "Red" line. Why is this so? The direction I calculate
for the "ResampleImageFilter" was essentially along the "Green" line.

Additionally, what I am trying to do, is this not possible with the
"ResampleImageFilter"? Do I need to extract the volume using the
"ExtractImageFilter" instead?

Thank you for your time
Sidharta



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

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http
Dženan Zukić
2017-05-11 15:40:56 UTC
Permalink
Hi Sidharta,

just changing the size of the output image should not affect its
orientation. Can you share the code with comments which line you added to
try changing the size?

Regards,
DÅŸenan Zukić, PhD, Senior R&D Engineer, Kitware (Carrboro, N.C.)
Dear DÅŸenan,
First of all, thank you for all your help.
So, from the previous posts I was able to resample the CBCT image so that it
lies in the direction of the Trajectory. So far so good. Now, I am trying to
extract only a part of the volume in this resampled image. So, I figure that
since I already calculated the direction, I could just set the size of the
output image to be what I want in the "ResampleImageFilter" and it should
work fine. However, when I do this, the extracted volume is in a different
direction than I assumed it to be. For example,
<http://itk-users.7.n7.nabble.com/file/n38228/Capture.png>
In the above image, I assumed the part of the volume that I extract is
actually in the direction shown in the "Green" line, however I instead get
the volume along the "Red" line. Why is this so? The direction I calculate
for the "ResampleImageFilter" was essentially along the "Green" line.
Additionally, what I am trying to do, is this not possible with the
"ResampleImageFilter"? Do I need to extract the volume using the
"ExtractImageFilter" instead?
Thank you for your time
Sidharta
--
View this message in context: http://itk-users.7.n7.nabble.
com/Extract-slice-from-ResampleFilter-tp38074p38228.html
Sent from the ITK - Users mailing list archive at Nabble.com.
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users
sidharta
2017-05-11 16:00:00 UTC
Permalink
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.
_____________________________________
Powered by www.kitware.com

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/i
Dženan Zukić
2017-05-12 14:55:06 UTC
Permalink
Hi Sidharta,

the older versions of ITK-SNAP had index-based overlay, ignoring different
spatial position which stem from different origin or direction vectors. I
don't know whether the latest (3.6) breaks with that.

Slicer always uses spatial position and orientation, so for your
experiments you should use Slicer to verify whether your code does what you
want.

This difference in size setting should not matter, I think your calculation
of the direction vectors is more likely cause of error. You could check
that by setting different sizes in between the two extremes you trying
right now. In case there is some bug in ITK (unlikely), you will flush it
out with this. More likely this will make it more apparent to you what you
are doing wrong.

Regards,
DÅŸenan Zukić, PhD, Senior R&D Engineer, Kitware (Carrboro, N.C.)
Dear DÅŸenan,
Thank you for your swift reply.
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.
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users
sidharta
2017-05-12 15:28:37 UTC
Permalink
Hi Dženan,

You are right. I fixed it a while back. The problem was infact in the
direction vectors. I realized (not a good practice) how ITK deals with
coordinate systems. I actually had to do an inverse of the direction matrix
and then everything is great. Thank you for your help.

Regards
Sidharta



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

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http://public.kitw

Loading...