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: [email protected] [mailto:[email protected]] Im Auftrag von Simon Rit
Gesendet: Freitag, 13. Mai 2016 18:15
An: Tobias Stein <[email protected]>
Cc: [email protected]
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 <[email protected]
<mailto:[email protected]> > 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: <mailto:[email protected]> [email protected] [mailto:
<mailto:[email protected]> [email protected]] Im Auftrag von Simon Rit
Gesendet: Mittwoch, 4. Mai 2016 09:02
An: Tobias Stein < <mailto:[email protected]>
[email protected]>
Cc: <mailto:[email protected]> [email protected]
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 <[email protected]
<mailto:[email protected]> > 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
[email protected] <mailto:[email protected]>
http://public.kitware.com/mailman/listinfo/rtk-users
_______________________________________________
Rtk-users mailing list
[email protected]
http://public.kitware.com/mailman/listinfo/rtk-users