[Rtk-users] forward and back projection - MITK
Simon Rit
simon.rit at creatis.insa-lyon.fr
Mon May 16 05:50:41 EDT 2016
Answers below. I think you need to understand better the coordinate system
of an ITK image, maybe check ITK's doc.
Simon
On Sun, May 15, 2016 at 2:38 PM, Tobias Stein <tstein at stud.hs-heilbronn.de>
wrote:
> If I understood you correct, the input of the Joseph-Filter needs to be
> oriented in the y-direction when I am dealing with axial oriented slices.
>
Yes, or in other words, if you want to take projections around the axial
axis, be aware that y is the default axis of rotation in RTK.
> Now I changed the direction of the 3D Slice by rotating it 90 degrees
> around the x-axis with the itk::ChangeInformationImageFilter to have my
> slice in the xz-layer in the y-direction:
>
> double rotationX = 90;
>
> itk::Versor< double > rotation;
>
> const double angleInRadians = rotationX * vnl_math::pi / 180.0;
>
> rotation.SetRotationAroundX(angleInRadians);
>
> const ImageType::DirectionType direction =
> inputImage->GetDirection();
>
> const ImageType::DirectionType newDirection = direction *
> rotation.GetMatrix();
>
> changeFilter->SetOutputDirection(newDirection);
>
> changeFilter->ChangeAll();
>
> changeFilter->SetInput(inputImage);
>
> changeFilter->UpdateOutputInformation();
>
>
>
> Then I added a zero-border as you told me (This changed the size of my
> region from [512, 512, 1] to [514, 514, 3]).
>
> While executing the Joseph Filter I got an error because a vector
> subscript was out of range. Maybe that’s because the index of the region is
> [-1,-1,-1] caused by the padding.
>
That should not have been a problem but that might be a bug...
> I tried the padding after the Joseph Filter, and got a result which has a
> size [514, 3, 45], while 45 was the number of projection images but the
> image was still empty.
>
I guess you should check the origin of your projected volume and your
projections. A simple test is to center your volume (check using ITK doc
how you can define the origin so that (0,0,0) is at the center) and your
projection (for the projection stack, only the first 2 coordinates count).
>
>
> I also tried the Jospeh Filter with a real 3D volume stack of slices,
> rotated it as described below and got a result.
>
> Now I am confused how I should create the geometry of the Joseph Filter to
> get a good result.
>
>
>
> typedef rtk::ThreeDCircularProjectionGeometry GeometryType;
>
> GeometryType::Pointer geometry = GeometryType::New();
>
> for (unsigned int i = 0; i < NumberOfProjectionImages; i++)
>
> geometry->AddProjection(300., 500., i * 180 /
> NumberOfProjectionImages);
>
>
>
> I am trying to use a 180-degree rotation, but i don’t have the feeling for
> the distance values.
>
> I this code example I used an 3D Image with dimensions = [256, 256, 49]
>
This doesn't tell you anything about the size of your volume, we also need
the origin and the spacing. Look at the ITK doc to understand the ITK
coordinate systems.
>
>
> My idea was, since there is no way to get a result with a single 3D slice,
> I would use a stack of slices as input and slice the output volume in
> y-direction to get my sinogram as you said.
>
> Tobias
>
>
>
> *Von:* simon.rit at gmail.com [mailto:simon.rit at gmail.com] *Im Auftrag von *Simon
> Rit
> *Gesendet:* Samstag, 14. Mai 2016 12:17
>
> *An:* Tobias Stein <tstein at stud.hs-heilbronn.de>
> *Cc:* rtk-users at public.kitware.com
> *Betreff:* Re: [Rtk-users] forward and back projection - MITK
>
>
>
> The output of the Joseph filter (of any forward projector) is a stack of
> projection images, i.e., a 3D image. If you slice this 3D image in the z
> direction, you get one projection images. If you slice this 3D image in the
> y direction, you get a sinogram, e.g., the 2D sinogram of the central slice
> (plane of the source trajectory) if y=0. The best way to have this sinogram
> is to set size[1]=1 and there is no better way to do it in RTK.
>
> For zero-padding, you need to set extendRegion[?]=1, with (for all
> dimensions, i.e., ?=0,1,2) and SetConstant(0). I'm surprised though that
> your image is a black square, something else must be going on. One thing to
> remember, the typical RTK rotation is around the y axis so your 3D slice
> must be 1 in the y direction. A simple way to change the orientation of one
> slice is to change its direction with ChangeImageFilter
>
> Simon
>
>
>
>
>
> On Sat, May 14, 2016 at 11:48 AM, Tobias Stein <
> tstein at stud.hs-heilbronn.de> wrote:
>
> Hello Simon,
>
>
>
> Thanks for your advice so far. I’m getting to know this filter better. At
> least I can now use it and save the result into a file.
>
>
>
> I updated my code which now uses a volume.
>
> The values of the input and output volume are equal except of the
> z-dimension of the output which I changed to the number of projection
> images.
>
> const unsigned int Dimension = 3;
>
> …Reading image…
>
> ImageType::Pointer inputImage = reader->GetOutput();
>
> ImageType::Pointer outputImage = ImageType::New();
>
>
>
> RegionType inputRegion = inputImage->GetLargestPossibleRegion();
>
>
>
> const unsigned int NumberOfProjectionImages = 45;
>
>
>
> ImageType::IndexType start;
>
> start[0] = inputRegion.GetIndex()[0];
>
> start[1] = inputRegion.GetIndex()[1];
>
> start[2] = inputRegion.GetIndex()[2];
>
>
>
> RegionType::SizeType size;
>
> size[0] = inputRegion.GetSize()[0];
>
> size[1] = inputRegion.GetSize()[1];
>
> size[2] = NumberOfProjectionImages;
>
>
>
> RegionType region;
>
> region.SetSize(size);
>
> region.SetIndex(start);
>
>
>
> ImageType::SpacingType spacing;
>
> spacing[0] = inputImage->GetSpacing()[0];
>
> spacing[1] = inputImage->GetSpacing()[1];
>
> spacing[2] = inputImage->GetSpacing()[2];
>
>
>
> ImageType::PointType origin;
>
> origin[0] = inputImage->GetOrigin()[0];
>
> origin[1] = inputImage->GetOrigin()[1];
>
> origin[2] = inputImage->GetOrigin()[2];
>
>
>
> outputImage->SetRegions(region);
>
> outputImage->SetSpacing(spacing);
>
> outputImage->SetOrigin(origin);
>
> outputImage->Allocate();
>
> outputImage->FillBuffer(size[0] * size[1] * size[2]);
>
>
>
> // Geometry
>
> typedef rtk::ThreeDCircularProjectionGeometry GeometryType;
>
> GeometryType::Pointer geometry = GeometryType::New();
>
> for (unsigned int i = 0; i < NumberOfProjectionImages; i++)
>
> geometry->AddProjection(510., 510., i*8.);
>
>
>
> typedef rtk::JosephForwardProjectionImageFilter<ImageType,
> ImageType> ForwardType;
>
> ForwardType::Pointer forward = ForwardType::New();
>
> forward->SetInput(1, inputImage);
>
> forward->SetInput(0, outputImage);
>
> forward->SetGeometry(geometry);
>
> forward->Update();
>
>
>
> I’m not sure how I have to add projection for the geometry to get a result
> which I’m familiar with.
>
>
>
> Can you describe how do I have to interpret the result of the
> Joseph-Filter? I am familiar with sinograms and how they are computed, but
> one sinogram is created by the projections of one image slice as far as I
> understood it. What I want is to manipulate the image in the radon space
> and then perform a filtered back projection. Typically the articles
> describe techniques to reduce metal artifacts by using one 2D Slice. Is it
> possible to achieve something like that with the JosephFilter? I’m guessing
> since this filter won’t work well for 2D I have to use another filter. Can
> you give me a hint if I can do something like that with RTK or another
> ITK-based framework? Here <https://i.imgur.com/1yBz3o2.png> a diagram
> which shows an example what I want to try.
>
>
>
> I tried also to use my code with one slice and added a zero padding like
> that:
>
> typedef itk::ConstantPadImageFilter <ImageType, ImageType>
> ConstantPadImageFilterType;
>
> ConstantPadImageFilterType::Pointer padFilter =
> ConstantPadImageFilterType::New();
>
> padFilter->SetInput(inputImage);
>
>
>
> ImageType::SizeType extendRegion;
>
> extendRegion[0] = 0;
>
> extendRegion[1] = 0;
>
> extendRegion[2] = 0;
>
>
>
> padFilter->SetPadLowerBound(extendRegion);
>
> padFilter->SetPadUpperBound(extendRegion);
>
> padFilter->SetConstant(1);
>
> padFilter->Update();
>
> But the result of my png is a black square.
>
> Tobias
>
>
>
>
>
> *Von:* simon.rit at gmail.com [mailto:simon.rit at gmail.com] *Im Auftrag von *Simon
> Rit
> *Gesendet:* Freitag, 13. Mai 2016 18:15
>
>
> *An:* Tobias Stein <tstein at stud.hs-heilbronn.de>
> *Cc:* rtk-users at public.kitware.com
> *Betreff:* Re: [Rtk-users] forward and back projection - MITK
>
>
>
> Dear Tobias,
>
> If you intialize the data of your itk::Image with FillBuffer with the
> corresponding constant value, it is quite similar. The constant image value
> allow us to stream the computation whereas you have to allocate the whole
> image if you pass directly itk::Image.
>
> We only work in 3D and probably our projector don't work with 2D, that is
> probably the problem here. You would have to correct the code for 2D but
> that's a lot of work, we generally prefer pseudo 2D, see next point.
>
> I doubt that our forward projector will work well with one slice only if
> you use 3D, we assume that the volume starts at the first pixel and ends at
> the last pixel. You should probably add a zero border around with
> itk::PadImage to see something.
>
> To avoid aliasing, you can cut frequencies in the ramp filter (look for
> CutFrequency in the doxygen of rtk::FFTRampImageFilter).
>
> Simon
>
>
>
> On Fri, May 13, 2016 at 5:35 PM, Tobias Stein <tstein at stud.hs-heilbronn.de>
> wrote:
>
> Thanks for the information.
>
>
>
> To keep it simple I want to project forward a single slice which I load as
> a DICOM file.
>
> Now I want to compute the sinogram of that slice and write it in a new
> file.
>
> I got a little confused about the ConstantImageSource which is used in the
> test class. Do I need it to compute a sinogram or is just a replacement of
> a itk::Image?
>
>
>
> Here is my code so far:
>
>
>
> const unsigned int Dimension = 2;
>
> typedef float PixelType;
>
> typedef itk::Image< PixelType, Dimension > ImageType;
>
> …
>
> (reading image by itk::ImageFileReader which works)
>
> …
>
> ImageType::Pointer inputImage = reader->GetOutput();
>
>
>
> const unsigned int NumberOfProjectionImages = 1;
>
>
>
> ImageType::Pointer outputImage = ImageType::New();
>
>
>
> typedef rtk::JosephForwardProjectionImageFilter<ImageType,
> ImageType> ForwardType;
>
> ForwardType::Pointer forward = ForwardType::New();
>
> forward->SetInput(1, inputImage);
>
> forward->SetInput(0, outputImage);
>
> forward->Update();
>
>
>
> …
>
> (writing output image into a file)
>
> …
>
>
>
> The problem is that I get some errors when I instantiate the
> forward-Filter.
>
> The errors are like that one:
>
> error C2784: "vnl_matrix<T> operator *(const vnl_diag_matrix<T> &,const
> vnl_matrix<T> &)": template-Argument für "const vnl_diag_matrix<T> &"
> konnte nicht von "vnl_matrix_fixed<T,4,4>" hergeleitet werden.
>
>
>
> Maybe there is something wrong with the PixelType that i used?
>
>
>
> If the forward projection works, which filter should i use to achieve a
> filtered back projection without aliasing?
>
>
>
> Many thanks in advance.
>
> Tobias
>
>
>
> *Von:* simon.rit at gmail.com [mailto:simon.rit at gmail.com] *Im Auftrag von *Simon
> Rit
> *Gesendet:* Mittwoch, 4. Mai 2016 09:02
> *An:* Tobias Stein <tstein at stud.hs-heilbronn.de>
> *Cc:* rtk-users at public.kitware.com
> *Betreff:* Re: [Rtk-users] forward and back projection - MITK
>
>
>
> Dear Tobias,
>
> The forward projection has 2 inputs:
>
> - input 0: the stack of projections in which you wish to forward project,
>
> - input 1: the volume you wish to forward project.
>
> For the backprojection, it's exactly the same:
>
> - input 0: the volume in which you wish to backproject,
>
> - input 1: the stack of projections you wish to backproject.
>
> Be aware that JosephBackProjectionImageFilter is the transpose of the
> forward projection and will have some aliasing (see, e.g., De Man and Basu
> <http://iopscience.iop.org/article/10.1088/0031-9155/49/11/024/meta;jsessionid=87B598ABFDC2AA07D1DED2DEE607F0E7.c2.iopscience.cld.iop.org>
> to see what I mean).
>
> I have personally never used MITK but don't hesitate to share your
> experience on the mailing list.
>
> Simon
>
>
>
> On Wed, May 4, 2016 at 12:59 AM, Tobias Stein <tstein at stud.hs-heilbronn.de>
> wrote:
>
> Hi all,
>
>
>
> i want to use the forward and backward projection to reduce metal
> artefacts in ct images. I’ve seen so far that I may use the
> JosephForwardProjectionImageFilter to perform the forward projection. I’ve
> also seen the test for this class but I don’t get it, where should i put my
> 2D slice as input to execute the transformation. About the back
> transformation with the JosephBackProjectionImageFilter I also need more
> information how I can use it to transform a sinogram back to a ct slice.
> There is a test at the documentation, but the link is missing there.
>
>
>
> I’ve got another independent question.
>
> Are some of you familiar with MITK and know how to get the superbuild up
> and running with RTK? I am writing a MITK-Plugin and want to reduce the
> metal artifacts before a segmentation. So if some of you have experience
> with the combination of RTK and MITK it would be nice if you will share it
> ;)
>
>
>
> Best regards,
>
> Tobias
>
>
>
>
> _______________________________________________
> Rtk-users mailing list
> Rtk-users at public.kitware.com
> http://public.kitware.com/mailman/listinfo/rtk-users
>
>
>
>
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/rtk-users/attachments/20160516/bfd7330e/attachment-0010.html>
More information about the Rtk-users
mailing list