[Rtk-users] inline Recon help
Simon Rit
simon.rit at creatis.insa-lyon.fr
Wed Dec 3 11:27:28 CET 2025
Please refrain to address personal messages, use the mailing list (cc).
This code is incomplete, there is no information on extractor. How is it
initialized? You need to check that the third coordinate of the index is
incremented from frame to frame.
Simon
On Wed, Dec 3, 2025 at 8:13 AM 贾雪辉 <jxhzk1314 at 163.com> wrote:
> In src code,inline test demo code run ok, and first with last frame the
> parker is 0,and the other projections is non-zero.
> But I transfer this demo to a interface project,make the inlinerecon demo
> as a thread when user call me, they only send me data. after first frame,
> the parker output is still 0!!!
> below is my code,can you help me figure out the wrong usage? thank you!
> void reconloop()
> {
> itk::MetaImageIOFactory::RegisterOneFactory();
> itk::GDCMImageIOFactory::RegisterOneFactory();
> itk::FFTWFFTImageFilterInitFactory::RegisterOneFactory();
> #ifdef RTK_USE_CUDA
> // Use CUDA for Parker short scan image filter
> auto ddf = rtk::CudaDisplacedDetectorImageFilter::New();
> auto parker = rtk::CudaParkerShortScanImageFilter::New();
> ddf->SetGeometry(m_pItkPointers->m_pCorrectGeometry);
> ddf->SetDisable(0);
> parker->SetGeometry(m_pItkPointers->m_pCorrectGeometry);
> parker->InPlaceOff();
> auto reconstructionSource =
> rtk::ConstantImageSource<FCudaImageType>::New();
> reconstructionSource->SetOrigin(itk::MakePoint(m_pReconInfo.imageOrigin[0],
> m_pReconInfo.imageOrigin[1], m_pReconInfo.imageOrigin[2]));
> reconstructionSource->SetSpacing(itk::MakeVector(m_pReconInfo.imageSpace[0],
> m_pReconInfo.imageSpace[1], m_pReconInfo.imageSpace[2]));
> reconstructionSource->SetSize(itk::MakeSize(m_pReconInfo.imageDimension[0],
> m_pReconInfo.imageDimension[1], m_pReconInfo.imageDimension[2]));
> itk::Matrix<double, 3, 3> TDimageDirection;
> TDimageDirection[0][0] = 1;
> TDimageDirection[1][2] = 1;
> TDimageDirection[2][1] = 1;
> reconstructionSource->SetDirection(TDimageDirection);
> TRY_AND_EXIT_ON_ITK_EXCEPTION(reconstructionSource->UpdateOutputInformation());
>
> auto fdk = rtk::CudaFDKConeBeamReconstructionFilter::New();
> #else
> // Use CPU for Parker short scan image filter
> auto parker = rtk::ParkerShortScanImageFilter<ImageType>::New();
> parker->SetInput(extractor->GetOutput());
> parker->SetGeometry(geometryRec);
> auto reconstructionSource = rtk::ConstantImageSource<ImageType>::New();
> reconstructionSource->SetOrigin(itk::MakePoint(originValue, originValue,
> originValue));
> reconstructionSource->SetSpacing(itk::MakeVector(spacingValue,
> spacingValue, spacingValue));
> reconstructionSource->SetSize(itk::MakeSize(imageSize, imageSize,
> imageSize));
> auto fdk = rtk::FDKConeBeamReconstructionFilter<ImageType>::New();
> #endif
> fdk->SetGeometry(m_pItkPointers->m_pCorrectGeometry);
> fdk->SetInput(0, reconstructionSource->GetOutput());
> fdk->SetInput(1, parker->GetOutput());
> //process in filter
> fdk->GetRampFilter()->SetHannCutFrequency(0.46);
> fdk->GetRampFilter()->SetHannCutFrequencyY(0.46);
> fdk->GetRampFilter()->SetTruncationCorrection(m_pCorrectionInfo.padding);
> // Online reconstruction loop
> int reconstructedProjectionsCount = 0;
> while (reconstructedProjectionsCount != totalFrameN)
> {
> //auto start = std::chrono::high_resolution_clock::now();
> if (reconstructedProjectionsCount <= currentFrameN)
> {
> std::cout << "Processing projection #" << reconstructedProjectionsCount
> << "\r";
> if (reconstructedProjectionsCount == 0)
> {
> reconReady = false;
> }
> else
> {
> #ifdef RTK_USE_CUDA
> FCudaImageType::Pointer reconstructedImage = fdk->GetOutput();
> #else
> ImageType::Pointer reconstructedImage = fdk->GetOutput();
> #endif
> reconstructedImage->DisconnectPipeline();
> fdk->SetInput(reconstructedImage);
> }
> #ifdef RTK_USE_CUDA
> m_pItkPointers->m_pInputImageFilter = itk::ImportImageFilter<
> InputPixelType, DIMENSION_3D >::New();
> m_pItkPointers->m_pInputImageTypeCuda = InputImageTypeCuda::New();
> m_pItkPointers->m_pTempImage = InputImage3DType::New();
> InputImageTypeCuda::SizeType size;
> size[0] = m_pDicomImageInfo->nImageRowCount; //width! so,input must
> right.
> size[1] = m_pDicomImageInfo->nImageColumnCount;//height
> size[2] = 1;// for inline case nFrame is 1.
> InputImageTypeCuda::IndexType start;
> start.Fill(0);
> InputImageTypeCuda::RegionType region;
> region.SetSize(size);
> region.SetIndex(start);
> m_pItkPointers->m_pInputImageFilter->SetRegion(region);
> m_pItkPointers->m_pInputImageTypeCuda->SetRegions(region);
> float origin[3] = { 0.0f };
> memcpy(origin, m_pDicomImageInfo->fgImagePositionPatient,
> sizeof(origin));
> m_pItkPointers->m_pInputImageFilter->SetOrigin(origin);
> m_pItkPointers->m_pInputImageTypeCuda->SetOrigin(origin);
> float spacing[3] = { 0.0f };
> if (memcmp(spacing, m_pDicomImageInfo->fgImagePixel,
> sizeof(m_pDicomImageInfo->fgImagePixel)) != 0)
> {
> memcpy(spacing, m_pDicomImageInfo->fgImagePixel,
> sizeof(m_pDicomImageInfo->fgImagePixel));
> }
> else if (memcmp(spacing, m_pDicomImageInfo->fgPixelSpacing,
> sizeof(m_pDicomImageInfo->fgPixelSpacing)) != 0)
> {
> memcpy(spacing, m_pDicomImageInfo->fgPixelSpacing,
> sizeof(m_pDicomImageInfo->fgPixelSpacing));
> }
> spacing[2] = 1;
> m_pItkPointers->m_pInputImageFilter->SetSpacing(spacing);
> m_pItkPointers->m_pInputImageTypeCuda->SetSpacing(spacing);
> m_pItkPointers->m_pInputImageTypeCuda->Allocate(true);
> m_pItkPointers->m_pInputImageTypeCuda->Update();
> float* tempData = new float[size[0] * size[1]]();
> int pointCount = size[0] * size[1];
> // here I use the data user send to me and I push back them as a list.
> unsigned short* t = (unsigned
> short*)frameDataList[reconstructedProjectionsCount];// based on
> reconstructedProjectionsCount
> for (size_t i = 0; i < pointCount; i++)
> {
> tempData[i] = (float)(((unsigned short*)t)[i]);
> }
> m_pItkPointers->m_pInputImageFilter->SetImportPointer(tempData,
> size[0] * size[1], true);
> m_pItkPointers->m_pInputImageFilter->Update();
> m_pItkPointers->m_pTempImage =
> m_pItkPointers->m_pInputImageFilter->GetOutput();
> m_pItkPointers->m_pTempImage->Update();
>
> m_pItkPointers->m_pInputImageTypeCuda->SetPixelContainer(m_pItkPointers->m_pTempImage->GetPixelContainer());
>
>
> //m_pItkPointers->m_pInputImageTypeCuda->CopyInformation(m_pItkPointers->m_pTempImage);
>
>
> //m_pItkPointers->m_pInputImageTypeCuda->SetBufferedRegion(m_pItkPointers->m_pTempImage->GetRequestedRegion());
>
>
> //m_pItkPointers->m_pInputImageTypeCuda->SetRequestedRegion(m_pItkPointers->m_pTempImage->GetRequestedRegion());
>
> m_pItkPointers->m_pInputImageTypeCuda->Update();
> // Prepare projection for CUDA processing
> using DuplicatorType = itk::ImageDuplicator<FCudaImageType>;
> DuplicatorType::Pointer duplicator = DuplicatorType::New();
> duplicator->SetInputImage(m_pItkPointers->m_pInputImageTypeCuda);
> duplicator->Update();
> FCudaImageType::Pointer projection = duplicator->GetOutput();//try
> using a deep copy but problem still there
> itk::WriteImage(projection,
> "before_minuslog-like-in-test.mha");//always ok
> calcMinusLogFrame(projection);
> itk::WriteImage(projection,
> "after_minuslog-like-in-test.mha");//always ok
> ddf->SetInput(projection);
> itk::WriteImage(ddf->GetOutput(),
> "after_ddf-like-in-test.mha");//always ok
> parker->SetInput(ddf->GetOutput());
> itk::WriteImage(parker->GetOutput(),
> "after_parker-like-in-test.mha");//why here always zero?
> int tmp = 1;
> #endif
> TRY_AND_EXIT_ON_ITK_EXCEPTION(fdk->Update());
> reconstructedProjectionsCount++;
> }
> else
> {
> // Sleep briefly to avoid busy waiting
> std::this_thread::sleep_for(std::chrono::milliseconds(10));
> }
> }
> huData = fdk->GetOutput();
> reconReady = true;
> }
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20251203/39a76cad/attachment-0001.htm>
More information about the Rtk-users
mailing list