[Rtk-users] forward and back projection - MITK
Simon Rit
simon.rit at creatis.insa-lyon.fr
Sat May 14 06:17:18 EDT 2016
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/20160514/5aaf2cfb/attachment-0010.html>
More information about the Rtk-users
mailing list