I tried and I could not understand what you've done in your code. Enclosed
is what I would use which gives a more reasonable result but which is full
of cone-beam artefacts (I think).
On Wed, Jun 3, 2020 at 7:59 PM Benjamin W. Maloney <
benjamin.w.maloney...@dartmouth.edu> wrote:
> Thanks.
>
>
> https://drive.google.com/drive/folders/1DTREuI80fWrIcBa5LIUfnAzsEPSRvhyZ?usp=sharing
>
> At that link are a .tif file for my projections and the .cxx
> reconstruction file I have been using
>
> Some notes about my system:
> SID: 306.6 mm
> SDD: 561.3 mm
>
> The detector panel is:
> 230mm x 70mm with a pixel matrix of 3072 x 864
> The projections are binned by 4 resulting in a pixel size of ~0.03mm-1 and
> a pixel array of 768 x 216
> There are 720 projection angles with an angular step of 0.5 degrees
>
> The system also has the sample on a turntable with a fixed source and
> detector but I believe that should be mathematically equivalent
>
> Thanks again!
> Ben
> ------------------------------
> *From:* Simon Rit <simon....@creatis.insa-lyon.fr>
> *Sent:* Wednesday, June 3, 2020 1:04 PM
> *To:* Benjamin W. Maloney <benjamin.w.maloney...@dartmouth.edu>
> *Cc:* Chao Wu <wucha...@gmail.com>; Andreas Andersen <
> andreasg...@gmail.com>; rtk-users@public.kitware.com <
> rtk-users@public.kitware.com>
> *Subject:* Re: [Rtk-users] Reconstruction Artifact
>
> Yes, that's allowed, other people did it before.
>
> On Wed, Jun 3, 2020 at 6:28 PM Benjamin W. Maloney <
> benjamin.w.maloney...@dartmouth.edu> wrote:
>
> Hi all,
>
> Thanks for the help!
> I'll be sure to make sure my preprocessed data is handled correctly.
>
> The reconstruction is meant to be a square, it is of a CIRS ACR digital
> mammography phantom
>
> I'm trying to send an example of my data but I can't get it under 147 MB
> data limit, is there a preferred method for sending this? I could send a
> google drive link if that is allowed
>
> Thanks again!
> Ben
>
> ------------------------------
> *From:* Chao Wu <wucha...@gmail.com>
> *Sent:* Wednesday, June 3, 2020 11:20 AM
> *To:* Simon Rit <simon....@creatis.insa-lyon.fr>
> *Cc:* Andreas Andersen <andreasg...@gmail.com>;
> rtk-users@public.kitware.com <rtk-users@public.kitware.com>; Benjamin W.
> Maloney <benjamin.w.maloney...@dartmouth.edu>
> *Subject:* Re: [Rtk-users] Reconstruction Artifact
>
> If your projections record the number of X-ray photon counts instead of
> attenuation, an I0 must be set correctly and a logarithmic operation is
> needed before the projection data can be fed to the reconstruction loop.
> Both the I0 and the logarithmic operation can be handled by the RTK
> projection reader or manually, depending on your implementation. Zero and
> negative numbers must be coerced for the logarithmic operation by you if
> this is not the case in the RTK code you use.
>
> Regards,
> Chao
>
> Simon Rit <simon....@creatis.insa-lyon.fr> 于2020年6月3日周三 下午2:35写道：
>
> Yes, there is no such limitation as far as I know, you can use negative
> numbers and value above 1.
> Your result is really strange, it it supposed to be a square? I don't know
> what is the problem but that's clearly a geometry issue. We can always have
> a look if you're able to share some data.
> Simon
>
> On Wed, Jun 3, 2020 at 1:56 PM Andreas Andersen <andreasg...@gmail.com>
> wrote:
>
> I don't think there are any restrictions technically.
> You should be able to use negative values, as the main arithmetic is just
> a sum
> <https://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2FSimonRit%2FRTK%2Fblob%2Fmaster%2Finclude%2FrtkBackProjectionImageFilter.hxx%23L106&data=02%7C01%7CBenjamin.W.Maloney.TH%40dartmouth.edu%7C7090376e438d462c76cf08d807dff9f4%7C995b093648d640e5a31ebf689ec9446f%7C0%7C0%7C637268005857264400&sdata=ta2wjTGRMl8MRkv5qeqUZDYED%2FD%2BNtzkFEMdIHp7YfA%3D&reserved=0>
> and a multiplication
> <https://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2FSimonRit%2FRTK%2Fblob%2Fmaster%2Finclude%2FrtkBackProjectionImageFilter.hxx%23L115&data=02%7C01%7CBenjamin.W.Maloney.TH%40dartmouth.edu%7C7090376e438d462c76cf08d807dff9f4%7C995b093648d640e5a31ebf689ec9446f%7C0%7C0%7C637268005857269384&sdata=pGY1uH5hjzUepFxR1AiDFKUFAxwNQM8mPIVDxwe%2B7N8%3D&reserved=0>
> .
>
> The only restriction I can see is that this sum and multiplication should
> not overflow (or underflow for negative values) the underlying type of the
> output image, as that would be undefined behaviour. Over- and underflow is
> unlikely for float, unless you have extremely high values (see wikipedia
> for floating point range
> <https://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FFloating-point_arithmetic%23Range_of_floating-point_numbers&data=02%7C01%7CBenjamin.W.Maloney.TH%40dartmouth.edu%7C7090376e438d462c76cf08d807dff9f4%7C995b093648d640e5a31ebf689ec9446f%7C0%7C0%7C637268005857274380&sdata=0wwsrWW8ELu0WjD5Ed9pJcKPFNI%2FNrSrIcGOV45jdgA%3D&reserved=0>
> ).
>
> /Andreas
>
> __________________________________
>
> Andreas Gravgaard Andersen
>
> Danish Center for Particle Therapy,
>
> Aarhus University Hospital
>
> Palle Juul-Jensens Blvd. 99,
>
> 8200, Aarhus
>
> Mail: agravga...@protonmail.com
>
> Cell: +45 3165 8140
>
>
> On Wed, 3 Jun 2020 at 03:34, Benjamin W. Maloney <
> benjamin.w.maloney...@dartmouth.edu> wrote:
>
> Hi,
>
> I thought the same in regards to trying to rotating in the other
> direction. Unfortunately, that has a similar artifact but with the
> reconstruction flipped.
> Interestingly the overlap happens in a similar place but the internal
> structures are flipped
>
> 1. Thanks!
> 2. I should have worded that better. My projection images will be
> preprocessed in a float format. I wanted to check if there were
> restrictions on these float values. Can input image data have negative
> values or high values? Or are they expected to have values between 0 and 1?
> I ask because some of the tools I have used to do this preprocessing
> (outside of RTK) have given negative values or 'stretched' the data from 0
> to 255 before saving.
>
> Ben
>
> ------------------------------
> *From:* Simon Rit <simon....@creatis.insa-lyon.fr>
> *Sent:* Tuesday, June 2, 2020 5:28 PM
> *To:* Benjamin W. Maloney <benjamin.w.maloney...@dartmouth.edu>
> *Cc:* rtk-users@public.kitware.com <rtk-users@public.kitware.com>
> *Subject:* Re: [Rtk-users] Reconstruction Artifact
>
> Hi,
> Sometimes rotating in the wrong direction gives this kind of artefacts.
> It's quite possible that we don't use the same convention as other toolkits
> regarding this.
> For other questions:
> 1. Yes, you can use pixel as the unit. Then the image spacing should be 1
> obviously and indeed, sdd and sid should be in pixels.
> 2. I don't fully understand. If your data is the output of a count
> detector, then you either rely on RTK to guess the counts without object to
> compute the line integral, or your preprocess your projections to pass line
> integrals in a float format.
> Simon
>
> On Tue, Jun 2, 2020 at 4:51 PM Benjamin W. Maloney <
> benjamin.w.maloney...@dartmouth.edu> wrote:
>
> Hi all,
>
> Not sure if this is the right group to post to but here's my question:
>
> I have a code that pulls in my own projection images and uses
> the FDKConeBeamReconstructionFilter for reconstruction
>
> The reconstruction I am getting has an artifact where it looks like there
> are two overlapping objects rotated.
> I have used other reconstruction toolboxes (mainly TIGRE in MATLAB) with
> the similar geometry inputs and not had this issue. I suspect the
> difference is in the parts of the geometry that are set by default. My
> question is if anyone has seen this before and what input I should look
> into?
>
> I have a few more questions I may or may not be related:
> 1. I assume that since I set the origin etc of my images in pixel, sid and
> sdd should be in pixels as well?
> 2. Are there restrictions related the scalar values of the projection
> data? My data will be in detector counts rather than linear attenuation
> coefficients, is that okay?
> I have attached an image to show this issue. It is supposed to be a
> rectangular mammography phantom. It is a slice in XZ plane
>
> Thanks for the help!
> Ben
>
>
>
> _______________________________________________
> Rtk-users mailing list
> Rtk-users@public.kitware.com
> https://public.kitware.com/mailman/listinfo/rtk-users
> <https://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fpublic.kitware.com%2Fmailman%2Flistinfo%2Frtk-users&data=02%7C01%7CBenjamin.W.Maloney.TH%40dartmouth.edu%7C7090376e438d462c76cf08d807dff9f4%7C995b093648d640e5a31ebf689ec9446f%7C0%7C0%7C637268005857279367&sdata=rN0WtXgPU48YJ%2BDkxnnx2StLGljd1AGBRX6IaXgRaY0%3D&reserved=0>
>
> _______________________________________________
> Rtk-users mailing list
> Rtk-users@public.kitware.com
> https://public.kitware.com/mailman/listinfo/rtk-users
> <https://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fpublic.kitware.com%2Fmailman%2Flistinfo%2Frtk-users&data=02%7C01%7CBenjamin.W.Maloney.TH%40dartmouth.edu%7C7090376e438d462c76cf08d807dff9f4%7C995b093648d640e5a31ebf689ec9446f%7C0%7C0%7C637268005857284358&sdata=cbP3Anssi1CsT65ReQWUQACc4nprtjW%2Bx%2F279Q0GH1k%3D&reserved=0>
>
> _______________________________________________
> Rtk-users mailing list
> Rtk-users@public.kitware.com
> https://public.kitware.com/mailman/listinfo/rtk-users
> <https://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fpublic.kitware.com%2Fmailman%2Flistinfo%2Frtk-users&data=02%7C01%7CBenjamin.W.Maloney.TH%40dartmouth.edu%7C7090376e438d462c76cf08d807dff9f4%7C995b093648d640e5a31ebf689ec9446f%7C0%7C0%7C637268005857289350&sdata=L%2FfRy4pSoo7gTUgVUPupSI3gwHIV%2FAlDyTSTqVfPzIg%3D&reserved=0>
>
>
// RTK includes
#include <rtkConstantImageSource.h>
#include <rtkThreeDCircularProjectionGeometryXMLFileWriter.h>
#include <rtkRayEllipsoidIntersectionImageFilter.h>
#include <rtkFDKConeBeamReconstructionFilter.h>
// ITK includes
#include <itkImageFileWriter.h>
#include <itkImageFileReader.h>
#include <itkTIFFImageIO.h>
#include <itkImage.h>
int main()
{
using GeometryType = rtk::ThreeDCircularProjectionGeometry;
GeometryType::Pointer geometry = GeometryType::New();
unsigned int numberOfProjections = 720;
double firstAngle = 0;
double angularArc = 360;
const double sid = 306.6; // source to isocenter distance, 306.6mm, 0.3 pixel size, image binned by 4 in x and y
const double sdd = 561.3; // source to detector distance
for ( unsigned int noProj = 0; noProj < numberOfProjections; noProj++ )
{
double angle = firstAngle + noProj * angularArc / numberOfProjections;
geometry->AddProjection( sid, sdd, angle );
}
// Write the geometry to disk
rtk::ThreeDCircularProjectionGeometryXMLFileWriter::Pointer xmlWriter;
xmlWriter = rtk::ThreeDCircularProjectionGeometryXMLFileWriter::New();
xmlWriter->SetFilename( "C:\\RTK_Example\\FirstReconstruction\\data\\geometry.xml" );
xmlWriter->SetObject( geometry );
xmlWriter->WriteFile();
// Defines the image type
using ImageType = itk::Image< float, 3 >;
// Projection Images
using ConstantImageSourceType = rtk::ConstantImageSource< ImageType >;
ConstantImageSourceType::PointType origin;
ConstantImageSourceType::SpacingType spacing;
ConstantImageSourceType::SizeType sizeOutput;
// Binned x4
origin[0] = -0.5 * (768-1) * 0.3;
origin[1] = -0.5 * (216-1) * 0.3;
origin[2] = 0.;
sizeOutput[0] = 768;
sizeOutput[1] = 216;
sizeOutput[2] = numberOfProjections;
spacing.Fill( 0.3);
// Load in my data instead
typedef float PixelType; // Pixel type double rather than float
const unsigned char Dimension = 3;
typedef itk::Image< PixelType, Dimension > ImageType;
ImageType::Pointer image;
typedef itk::ImageFileReader< ImageType > ReaderType;
ReaderType::Pointer reader;
reader = ReaderType::New();
reader->SetImageIO( itk::TIFFImageIO::New() );
reader->SetFileName( "../data/Projections_tif_full.tif" ); // bin x4
reader->Update();
reader->GetOutput()->SetOrigin( origin );
reader->GetOutput()->SetSpacing( spacing );
// Create reconstructed image
ConstantImageSourceType::Pointer constantImageSource2 = ConstantImageSourceType::New();
// x4
sizeOutput.Fill( 768 );
sizeOutput[1] = 216;
origin.Fill( -0.5 * ( 768 - 1 ) * 0.3*sid/sdd );
origin[1] = ( -0.5 * ( 216 - 1 ) * 0.3*sid/sdd );
spacing.Fill( 0.3*sid/sdd );
constantImageSource2->SetOrigin( origin );
constantImageSource2->SetSpacing( spacing );
constantImageSource2->SetSize( sizeOutput );
constantImageSource2->SetConstant( 0. );
// FDK reconstruction
using FDKCPUType = rtk::FDKConeBeamReconstructionFilter< ImageType >;
FDKCPUType::Pointer feldkamp = FDKCPUType::New();
feldkamp->SetInput( 0, constantImageSource2->GetOutput() );
feldkamp->SetInput( 1, reader->GetOutput() );
feldkamp->SetGeometry( geometry );
feldkamp->GetRampFilter()->SetTruncationCorrection( 0. );
feldkamp->GetRampFilter()->SetHannCutFrequency( 0.0 );
// Writer
using WriterType = itk::ImageFileWriter< ImageType >;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName( "C:\\RTK_Example\\FirstReconstruction\\data\\reconstructionRTK.tif" );
writer->SetInput( feldkamp->GetOutput() );
writer->Update();
std::cout << "Done!" << std::endl;
}
_______________________________________________
Rtk-users mailing list
Rtk-users@public.kitware.com
https://public.kitware.com/mailman/listinfo/rtk-users
