i build ITK5.2.1&RTK2.3 in C++,RTK’s demo FirstReconstruction can run normally，but when i modified it for my project ,it failed,The detailed error message is the range of reconstructed data is [0.0,0.0],and display it in VTK shows nothing. my original data type is raw, which is named test.raw, the size of which is 30082496360 ,here are my mainly code,please do me a favor and give some advices,Thank you!
using PixelType = unsigned short int;
constexpr unsigned int Dimension = 3;
using ImageType = itk::Image<PixelType, Dimension>;
using GPUImageType = itk::CudaImage<float, 3>;
int width = 3008;
int height = 2496;
typedef itk::RawImageIO<unsigned short , 3> ImageIOType;
typedef itk::ImageFileReader<ImageType> ImageRType;
ImageIOType::Pointer rawImgIO = ImageIOType::New();
rawImgIO->SetFileTypeToBinary(); //设置二进制形式打开
rawImgIO->SetFileDimensionality(3); //设置图像维度
rawImgIO->SetByteOrderToLittleEndian(); //设置小端序
rawImgIO->SetDimensions(0, width);
rawImgIO->SetDimensions(1, height);
rawImgIO->SetDimensions(2, numOfProj);
ImageType::SpacingType Spacing;
Spacing[0] = spacing_origion; // spacing along X
Spacing[1] = spacing_origion; // spacing along Y
Spacing[2] = spacing_origion; // spacing along Z
rawImgIO->SetSpacing(0, Spacing[0]);
rawImgIO->SetSpacing(1, Spacing[1]);
rawImgIO->SetSpacing(2, Spacing[2]);
ImageType::PointType newOrigin;
newOrigin[0] = spacing_origion * (width - 1) * (0.5);
newOrigin[1] = spacing_origion * (height - 1) * 0.5;
newOrigin[2] = 0.0;
rawImgIO->SetOrigin(0, newOrigin[0]);
rawImgIO->SetOrigin(1, newOrigin[1]);
rawImgIO->SetOrigin(2, newOrigin[2]);
rawImgIO->SetHeaderSize(0); //文件头大小
rawImgIO->SetNumberOfComponents(1); //像素通道数(灰度为1，RGB为3)
rawImgIO->SetPixelType(itk::ImageIOBase::SCALAR); //设置标量类型，即灰度图
rawImgIO->SetComponentType(itk::ImageIOBase::IOComponentEnum::USHORT);//设置像素数值类型为ushort
ImageRType::Pointer reader = ImageRType::New();
reader->SetFileName(rawimages);
std::cout << "Read Data!" << std::endl;
reader->SetImageIO(rawImgIO);
//reader->Update();
using ShrinkImageFilterType = itk::ShrinkImageFilter<ImageType, ImageType>;
ShrinkImageFilterType::Pointer shrinkFilter = ShrinkImageFilterType::New();
shrinkFilter->SetInput(reader->GetOutput());
shrinkFilter->SetShrinkFactor(0, 4); // shrink the first dimension by a factor of 2 (i.e. 100 gets changed to 50)752*624
shrinkFilter->SetShrinkFactor(1, 4); // shrink the second dimension by a factor of 3 (i.e. 100 gets changed to 33)
shrinkFilter->SetShrinkFactor(2, 1);
std::cout << "shrinkFilter!" << std::endl;
/* shrinkFilter->Update();
GPUInputImageType::SpacingType Spacing1 = shrinkFilter->GetOutput()->GetSpacing();
*/
using GeometryType = rtk::ThreeDCircularProjectionGeometry;
GeometryType::Pointer geometry = GeometryType::New();
double firstAng = 0.0;
double angleRac = 360.0;
#pragma omp parallel for
for (int i = 0; i < numOfProj; i++)
{
double angle = firstAng + i * angleRac / numOfProj;
geometry->AddProjection(sid, sdd, angle, shift, 0.0, 0.0, deangle, 0.0, 0.0);
}
std::cout << "Create Geometry!" << std::endl;
using CastFilter = itk::CastImageFilter<ImageType, GPUImageType>;
CastFilter::Pointer Castfilter = CastFilter::New();
Castfilter->SetInput(shrinkFilter->GetOutput());
std::cout << "Castfilter!" << std::endl;
//filter->Update();
using ConstantImageType = rtk::ConstantImageSource<GPUImageType>;
ConstantImageType::PointType origion;
ConstantImageType::SpacingType spacing;
ConstantImageType::SizeType outSize;
float m = width / OutX;
float scale[] = { width / OutX,width / OutY,height / OutZ };
spacing[0] = 0.4 / (sdd / sid) * scale[0];
spacing[1] = 0.4 / (sdd / sid) * scale[1];
spacing[2] = 0.4 / (sdd / sid) * scale[2];
origion[0] = -spacing[0] * (OutX - 1) * 0.5;
origion[1] = -spacing[1] * (OutY - 1) * 0.5;
origion[2] = -spacing[2] * (OutZ - 1) * 0.5;
outSize[0] = OutX;
outSize[1] = OutY;
outSize[2] = OutZ;
ConstantImageType::Pointer constanttype = ConstantImageType::New();
constanttype->SetOrigin(origion);
constanttype->SetSpacing(spacing);
constanttype->SetSize(outSize);
constanttype->SetConstant(0.);
//constanttype->Update();
std::cout << "Create ConstantImage!" << std::endl;
using FDKFilter = rtk::CudaFDKConeBeamReconstructionFilter;
FDKFilter::Pointer fdkfilter = FDKFilter::New();
fdkfilter->SetInput(0, constanttype->GetOutput());
fdkfilter->SetInput(1, Castfilter->GetOutput());
fdkfilter->SetGeometry(geometry);
fdkfilter->GetRampFilter()->SetTruncationCorrection(0.);
fdkfilter->GetRampFilter()->SetHannCutFrequency(0.0);
//fdkfilter->Update();
std::cout << "fdkfilter!" << std::endl;
using FOVFilter = rtk::FieldOfViewImageFilter<GPUImageType, GPUImageType>;
FOVFilter::Pointer fovfilter = FOVFilter::New();
fovfilter->SetInput(0, fdkfilter->GetOutput());
fovfilter->SetProjectionsStack(Castfilter->GetOutput());
fovfilter->SetGeometry(geometry);
//fovfilter->Update();
std::cout << "fovfilter!" << std::endl;
/* vtkfilter->SetInput(fovfilter->GetOutput());*/
using WriterType = itk::ImageFileWriter<GPUImageType>;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName("res.mha");
writer->SetInput(fovfilter->GetOutput());
try
{
writer->Update();
std::cout << "Finished!" << std::endl;
}
catch (const itk::ExceptionObject& err)
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
}
