The result of calling the interpolation_projection shows an exception

[OHMYGODjack]

I am trying to call the 'interpolation_projection' function in the 'ray_interpolated_projection.cu' file to implement forward projection, but the projection results I get are abnormal. I set 100 angles, but in the resulting projections, only the first angle has values, while the projection values for all other angles are entirely 0. Additionally, I am using the 'head_phantom' data, but it appears that the orientation of the object being projected is different from that in the TIGRE examples. I would like to know what caused the abnormal projection results. Could you help me figure out if I made any mistakes? Could you provide an example of how to correctly call the 'interpolation_projection' function? Thank you.

the first angle projection：

the 34th angle projection:

the last angle projection:

the head_phantom data:

Code to reproduce the problem (If applicable)

#include <string.h>
#include "ray_interpolated_projection.hpp"
#include "GpuIds.hpp"
#include <cmath>
#include <fstream>
#include <iostream>
#include <string>

void Ax_projection_test()
{
    float* image_data;      
    Geometry geo;       
    float** result;     
    float* angles;      
    int nangles;        
    GpuIds gpuids;      

    std::string dataPath = R"(X:\Recon\Ax_test\OutputHead-float-64x64x64.raw)";
    int width = 64;
    int height = 64;
    int depth = 64;
    size_t length = width * height * depth * size_t(1);

//=========================================================================
    image_data = new float[length];
    bool ret = loadData<float>(image_data, length, 0, dataPath.c_str());

    nangles = 100;

    geo.nVoxelX = 64;    
    geo.nVoxelY = 64;
    geo.nVoxelZ = 64;

    geo.sVoxelX = 256;
    geo.sVoxelY = 256;
    geo.sVoxelZ = 256;

    geo.dVoxelX = geo.sVoxelX / geo.nVoxelX;
    geo.dVoxelY = geo.sVoxelY / geo.nVoxelY;
    geo.dVoxelZ = geo.sVoxelZ / geo.nVoxelZ;

    geo.DSD = new float[nangles] {1536.f};
    geo.DSO = new float[nangles] {1000.f};

    geo.nDetecU = 256;
    geo.nDetecV = 256;

    geo.sDetecU = 409.6;
    geo.sDetecV = 409.6;

    geo.dDetecU = geo.sDetecU / geo.nDetecU;
    geo.dDetecV = geo.sDetecV / geo.nDetecV;

    geo.offOrigX = new float[nangles] {0.f};
    geo.offOrigY = new float[nangles] {0.f};
    geo.offOrigZ = new float[nangles] {0.f};

    geo.offDetecU = new float[nangles] {0.f};
    geo.offDetecV = new float[nangles] {0.f};

    geo.accuracy = 0.5;

    geo.COR = new float[nangles] {0.f};

    geo.dRoll = new float[nangles] {0.f};
    geo.dPitch = new float[nangles] {0.f};
    geo.dYaw = new float[nangles] {0.f};

    geo.unitX = 1;
    geo.unitY = 1;
    geo.unitZ = 1;

    int deviceId = 0;
    gpuids.SetIds(1, &deviceId);
    gpuids.SetAllGpus(1);

    angles = new float[nangles * 3];
    float angel_delta = 3.1415926535 * 2 / nangles;
    for (int i = 0; i < nangles; i++)
    {
        angles[i * 3 + 0] = angel_delta * i;
        angles[i * 3 + 1] = 0.f;
        angles[i * 3 + 2] = 0.f;

        //std::cout << angles[i * 3 + 0] << std::endl;
    }

    unsigned long long nDetectorPixels = (unsigned long long)geo.nDetecU * geo.nDetecV;
    size_t total_projections = nangles; 
    float* outProjections = new float[nDetectorPixels * total_projections];
    result = new float* [total_projections];
    for (size_t i = 0; i < total_projections; ++i)
    {
        unsigned long long currProjIndex = nDetectorPixels * i;
        result[i] = &outProjections[currProjIndex];
    }

    interpolation_projection(image_data, geo, result, angles, nangles, gpuids);

    std::string save_name = "Projection-float-" + std::to_string(geo.nDetecU) + "x" + std::to_string(geo.nDetecV) + "x" + std::to_string(nangles) + ".raw";
    std::string save_path = "X:/Recon/Ax_test/" + save_name;
    saveData<float>(outProjections, nDetectorPixels * total_projections,
        save_path.c_str());

template<class T>
inline bool loadData(
    T* input, const size_t length, size_t offset ,const char* path)
{
    FILE* fp;
    fopen_s(&fp, path, "rb");
    if (nullptr != fp)
    {
        //fseek(fp, offset, SEEK_CUR);
        _fseeki64(fp, offset, SEEK_CUR);
        fread(input, 1, length * sizeof(T), fp);
        //_fread_nolock
        fclose(fp);
        return true;
    }
    else
    {
        std::cout<< "Open " << path << " failed!" << std::endl;
        return false;
    }

}

template<class T>
inline int saveData(
    T* input, const size_t length,const char* path
)
{
    FILE* fp;
    fopen_s(&fp, path, "wb");
    if (nullptr != fp)
    {
        fwrite(input, 1, length * sizeof(T), fp);
        fclose(fp);
        return true;
    }
    else
    {
        std::cout << "Save " << path << " failed!" << std::endl;
        return false;
    }
}

Specifications

• MATLAB version:
• OS:
• CUDA version:

[AnderBiguri]

hum, no idea. Is this the same for the Siddon one? That one should be faster too, I would recommend.

I can't see anything blatantly wrong. What I would suggest is to print the values of the variables (except image) before the call, and to compare that to the MATLAB/Python version, to try to see if anything is off, as the MATLAB/Python way of calling this should work.

For orientation, there is a compilation flag for IS_FOR_MATLAB_TIGRE that will control the direction. Ultimately we use the same source for MATLAB and Python and one has Fortran memory ordering and the other one C, so with this flag the order of operations are changed and the result may come in a different shape.

[OHMYGODjack]

When I call the siddon_ray_projection function to calculate projections, the first projection is correct, but the others are abnormal. During debugging, I compared the parameters in MATLAB and found that the geo parameters are the same. I'm unable to figure out where I made a mistake, which is frustrating.Notably, although the other angle projections are incorrect, they seem to exhibit structural textures.

the first angle projection：

the 41th angle projection:

the last angle projection:

[AnderBiguri]

Hum. It's quite hard to help without me setting up the project in the same way and start investigating the variables. It's almost certain that there is a big in your code, mostly because this works in python and MATLAB and I have myself called it directly from C++ before in some private projects, so I am certain it works.

Few things to try to figure out:

Are all the data types exactly the expected ones?

Are all the numerical values of the geoemtry exactly the ones that get passed into the function if you call it from standard Tigre?

Is it the first projection that always works, or the projection corresponding to that angle? Ie of you call it in a loop, do you get the right result?

These may help figure out what the issue is. Otherwise you may need to dig into the function and see what is being passed to individual calls. My guess is that some number is not what you think it is or that the memory is not being out in the right place (eg not allocated properly) . I say this because the code internslly computes several projections with the same GPU call, so it's strange that only 1 "computes". All of them must be properly being computed, but the geo data is wrong or the data is saved somewhere badly.

[AnderBiguri]

@OHMYGODjack did you find a solution?

[OHMYGODjack]

I've found the issue.just as you said, I had set the geometric parameters incorrectly, with only the first angle's geometric parameter being set correctly. Thank you very much for your help.