For the correct reconstruction result, there should be a rectangle in the middle, but the RTK reconstruction result looks like it is missing some data, becoming a quadrilateral. Is my projection data 845, covering 0-340 angles, and this reconstruction result exactly looks like it is missing 10 degrees of data? What is the problem? How should I modify it?
/*=========================================================================
*
* Copyright RTK Consortium
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* https://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#include "rtkfdk_ggo.h"
#include "rtkGgoFunctions.h"
#include "rtkConfiguration.h"
#include "rtkThreeDCircularProjectionGeometryXMLFile.h"
#include "rtkDisplacedDetectorForOffsetFieldOfViewImageFilter.h"
#include "rtkParkerShortScanImageFilter.h"
#include "rtkFDKConeBeamReconstructionFilter.h"
#ifdef RTK_USE_CUDA
# include "rtkCudaDisplacedDetectorImageFilter.h"
// TODO # include "rtkCudaDisplacedDetectorForOffsetFieldOfViewImageFilter.h"
# include "rtkCudaParkerShortScanImageFilter.h"
# include "rtkCudaFDKConeBeamReconstructionFilter.h"
#endif
#include "rtkFDKWarpBackProjectionImageFilter.h"
#include "rtkCyclicDeformationImageFilter.h"
#include "rtkProgressCommands.h"
#include <itkStreamingImageFilter.h>
#include <itkImageRegionSplitterDirection.h>
#include <itkImageFileWriter.h>
int
main(int argc, char* argv[])
{
GGO(rtkfdk, args_info);
using OutputPixelType = float;
constexpr unsigned int Dimension = 3;
using CPUOutputImageType = itk::Image<OutputPixelType, Dimension>;
#ifdef RTK_USE_CUDA
using OutputImageType = itk::CudaImage<OutputPixelType, Dimension>;
#else
using OutputImageType = CPUOutputImageType;
#endif
// Projections reader
using ReaderType = rtk::ProjectionsReader<OutputImageType>;
ReaderType::Pointer reader = ReaderType::New();
rtk::SetProjectionsReaderFromGgo<ReaderType, args_info_rtkfdk>(reader, args_info);
if (!args_info.lowmem_flag)
{
if (args_info.verbose_flag)
std::cout << "Reading... " << std::endl;
reader->Update();
TRY_AND_EXIT_ON_ITK_EXCEPTION(reader->Update())
}
//将读取到的投影图像另存一下
const char* proj = "projectiong.mha" ;
TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(reader->GetOutput(), proj))
std::string inputVolumeFile = "D:\\code\\CBCT\\threshold.mha";
std::string subprojection = "G:\\CBCT\\workflow\\CBCT_RTK\\subimage\\subprojection.mha";
std::string projection = "G:\\CBCT\\workflow\\CBCT_RTK\\CBCT_RTK\\projectiong.mha";
/* typedef itk::ImageFileReader<OutputImageType> ReaderType;
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName(projection);
reader->Update();*/
//// Geometry
/* rtk::ThreeDCircularProjectionGeometry::Pointer geometry;
TRY_AND_EXIT_ON_ITK_EXCEPTION(geometry = rtk::ReadGeometry(args_info.geometry_arg));
geometry->GetMagnificationMatrices();
for (auto x: geometry->GetSourcePosition(0))
{
std::cout <<"geometry->GetTiltAngles()"<< x << std::endl;;
}*/
// geometry->GetAngularGaps();
if (args_info.verbose_flag)
std::cout << "Reading geometry information from " << args_info.geometry_arg << "..." << std::endl;
// rtk::ThreeDCircularProjectionGeometry::Pointer geometry;
// TRY_AND_EXIT_ON_ITK_EXCEPTION(geometry = rtk::ReadGeometry("projection_matrix.xml"));
// geometry->GetMagnificationMatrices();
constexpr unsigned int NumberOfProjectionImages = 845;
using GeometryType = rtk::ThreeDCircularProjectionGeometry;
GeometryType::Pointer geometry = GeometryType::New();
// geometry->SetCollimationOfLastProjection();
double sid = 1000;
double sad =1.5* sid;
for (unsigned int noProj = 0; noProj < NumberOfProjectionImages; noProj++)
geometry->AddProjection(996.9, 1500.5, (noProj * 339.6 / NumberOfProjectionImages)-80,0,0,0 ,0 , 0, 0);
// geometry->AddProjection(1000, 1500, (noProj * 339.2 / NumberOfProjectionImages),3.5,10.5, 0, 0, 0, 0);
// geometry->AddProjectionInRadians();
// (noProj * 360.0 / NumberOfProjectionImages)
const char fileName[] = "rtkgeometryfiletest.out";
rtk::ThreeDCircularProjectionGeometryXMLFileWriter::Pointer xmlWriter =
rtk::ThreeDCircularProjectionGeometryXMLFileWriter::New();
xmlWriter->SetFilename(fileName);
xmlWriter->SetObject(geometry);
TRY_AND_EXIT_ON_ITK_EXCEPTION(xmlWriter->WriteFile())
// Check on hardware parameter
#ifndef RTK_USE_CUDA
if (!strcmp(args_info.hardware_arg, "cuda"))
{
std::cerr << "The program has not been compiled with cuda option" << std::endl;
return EXIT_FAILURE;
}
#endif
// Displaced detector weighting 暂时没用到
using DDFCPUType = rtk::DisplacedDetectorImageFilter<OutputImageType>;
using DDFOFFFOVType = rtk::DisplacedDetectorForOffsetFieldOfViewImageFilter<OutputImageType>;
#ifdef RTK_USE_CUDA
using DDFType = rtk::CudaDisplacedDetectorImageFilter;
#else
using DDFType = rtk::DisplacedDetectorForOffsetFieldOfViewImageFilter<OutputImageType>;
#endif
DDFCPUType::Pointer ddf;
if (!strcmp(args_info.hardware_arg, "cuda"))
ddf = DDFType::New();
else
ddf = DDFOFFFOVType::New();
ddf->SetInput(reader->GetOutput());
ddf->SetGeometry(geometry);
ddf->SetDisable(0);
// Short scan image filter 短扫描 阈值设置不是很明白
using PSSFCPUType = rtk::ParkerShortScanImageFilter<OutputImageType>;
#ifdef RTK_USE_CUDA
using PSSFType = rtk::CudaParkerShortScanImageFilter;
#else
using PSSFType = rtk::ParkerShortScanImageFilter<OutputImageType>;
#endif
PSSFCPUType::Pointer pssf;
if (!strcmp(args_info.hardware_arg, "cuda"))
pssf = PSSFType::New();
else
pssf = PSSFCPUType::New();
pssf->SetInput(reader->GetOutput());
pssf->SetGeometry(geometry);
pssf->InPlaceOn();
pssf->SetAngularGapThreshold(0.35);
pssf->Update();
// Create reconstructed image 重建图像源
using ConstantImageSourceType = rtk::ConstantImageSource<OutputImageType>;
ConstantImageSourceType::Pointer constantImageSource = ConstantImageSourceType::New();
rtk::SetConstantImageSourceFromGgo<ConstantImageSourceType, args_info_rtkfdk>(constantImageSource, args_info);
// Motion-compensated objects for the compensation of a cyclic deformation.
// Although these will only be used if the command line options for motion
// compensation are set, we still create the object before hand to avoid auto
// destruction.
using DVFPixelType = itk::Vector<float, 3>;
using DVFImageSequenceType = itk::Image<DVFPixelType, 4>;
using DVFImageType = itk::Image<DVFPixelType, 3>;
using DeformationType = rtk::CyclicDeformationImageFilter<DVFImageSequenceType, DVFImageType>;
using DVFReaderType = itk::ImageFileReader<DeformationType::InputImageType>;
DVFReaderType::Pointer dvfReader = DVFReaderType::New();
DeformationType::Pointer def = DeformationType::New();
def->SetInput(dvfReader->GetOutput());
using WarpBPType = rtk::FDKWarpBackProjectionImageFilter<OutputImageType, OutputImageType, DeformationType>;
WarpBPType::Pointer bp = WarpBPType::New();
bp->SetDeformation(def);
bp->SetGeometry(geometry);
// This macro sets options for fdk filter which I can not see how to do better
// because TFFTPrecision is not the same, e.g. for CPU and CUDA (SR)
#define SET_FELDKAMP_OPTIONS(f) \
f->SetInput(0, constantImageSource->GetOutput()); \
f->SetInput(1, pssf->GetOutput()); \
f->SetGeometry(geometry); \
f->GetRampFilter()->SetTruncationCorrection(0.0); \
f->GetRampFilter()->SetHannCutFrequency(0.5); \
f->GetRampFilter()->SetHannCutFrequencyY(0.8); \
f->SetProjectionSubsetSize(args_info.subsetsize_arg); \
if (args_info.verbose_flag) \
{ \
f->AddObserver(itk::AnyEvent(), progressCommand); \
}
// FDK reconstruction filtering FDK
using FDKCPUType = rtk::FDKConeBeamReconstructionFilter<OutputImageType>;
FDKCPUType::Pointer feldkamp;
#ifdef RTK_USE_CUDA
using FDKCUDAType = rtk::CudaFDKConeBeamReconstructionFilter;
FDKCUDAType::Pointer feldkampCUDA;
#endif
itk::Image<OutputPixelType, Dimension>* pfeldkamp = nullptr;
if (!strcmp(args_info.hardware_arg, "cpu"))
{
// Progress reporting
using PercentageProgressCommandType = rtk::PercentageProgressCommand<FDKCPUType>;
PercentageProgressCommandType::Pointer progressCommand = PercentageProgressCommandType::New();
feldkamp = FDKCPUType::New();
SET_FELDKAMP_OPTIONS(feldkamp);
// Motion compensated CBCT settings
if (args_info.signal_given && args_info.dvf_given)
{
dvfReader->SetFileName(args_info.dvf_arg);
def->SetSignalFilename(args_info.signal_arg);
feldkamp->SetBackProjectionFilter(bp.GetPointer());
}
pfeldkamp = feldkamp->GetOutput();
}
#ifdef RTK_USE_CUDA
else if (!strcmp(args_info.hardware_arg, "cuda"))
{
// Motion compensation not supported in cuda
if (args_info.signal_given && args_info.dvf_given)
{
std::cerr << "Motion compensation is not supported in CUDA. Aborting" << std::endl;
return EXIT_FAILURE;
}
// Progress reporting
using PercentageProgressCommandType = rtk::PercentageProgressCommand<FDKCUDAType>;
PercentageProgressCommandType::Pointer progressCommand = PercentageProgressCommandType::New();
feldkampCUDA = FDKCUDAType::New();
SET_FELDKAMP_OPTIONS(feldkampCUDA);
pfeldkamp = feldkampCUDA->GetOutput();
}
#endif
// Streaming depending on streaming capability of writer
using StreamerType = itk::StreamingImageFilter<CPUOutputImageType, CPUOutputImageType>;
StreamerType::Pointer streamerBP = StreamerType::New();
streamerBP->SetInput(pfeldkamp);
streamerBP->SetNumberOfStreamDivisions(args_info.divisions_arg);
itk::ImageRegionSplitterDirection::Pointer splitter = itk::ImageRegionSplitterDirection::New();
splitter->SetDirection(2); // Prevent splitting along z axis. As a result, splitting will be performed along y axis
streamerBP->SetRegionSplitter(splitter);
// Write
if (args_info.verbose_flag)
std::cout << "Reconstructing and writing... " << std::endl;
TRY_AND_EXIT_ON_ITK_EXCEPTION(itk::WriteImage(pfeldkamp, args_info.output_arg))
return EXIT_SUCCESS;
}