How to allocate a stack of images in a P1 function in C++?

mesh
1

Hi all, I am using dolfin/C++.

What I have: a set of 2D images (pgm format), each with the same resolution, deriving from CT DICOM;

What I want to do: create a scalar 3D P1 function which allocates the stack of 2D images (a so-called 3D-image).

What I did:

std::shared_ptr<dolfin::Function> setupDICOM::amReadDicoms(std::string name, int n1, int n2, double dx, double dy, double dz)
{
    int imageWidth, imageHeight;
    std::ifstream infile(name+"/slice"+std::to_string(n1)+".pgm");
    std::stringstream ss;
    std::string inputLine = "";
    getline(infile,inputLine);
    getline(infile,inputLine);
    std::string spool = "";
    spool += inputLine[0]; spool += inputLine[1]; spool += inputLine[2];
    std::istringstream iss1(spool);
    iss1 >> imageWidth;
    spool = "";
    spool += inputLine[4]; spool += inputLine[5]; spool += inputLine[6];
    std::istringstream iss2(spool);
    iss2 >> imageHeight;
    ss << infile.rdbuf();

    auto p1 = std::make_shared<dolfin::Point>(0.0,0.0,0.0);
    auto p2 = std::make_shared<dolfin::Point>(dx*imageWidth,dy*imageHeight,dz*(n2-n1));

    auto Th = std::make_shared<dolfin::BoxMesh>(*p1,*p2,imageWidth-1,imageHeight-1,n2-n1-1);
    auto Vh = std::make_shared<varfs::CoefficientSpace_mu>(Th);

    auto image = std::make_shared<dolfin::Function>(Vh);

    std::vector<double> unrolledImge;

    for(int kkl=0; kkl<n2-n1; kkl++)
    {
        int index = n1+kkl;
        std::string str = name+"slice";
        str+=std::to_string(index);
        str+=".pgm";
        std::ifstream infileInner(str);
        if(infileInner.fail())
        {
            return nullptr;
        }
        std::stringstream ssInner;
        std::string inputLineInner = "";
        getline(infileInner,inputLineInner);
        getline(infileInner,inputLineInner);
        ssInner << infileInner.rdbuf();
        int spool2;
        ssInner >> spool2;
        for(int ik=0; ik<imageHeight; ik++)
        {
            for(int jk=0; jk<imageWidth; jk++)
            {
                unsigned char px;
                ssInner >> px;
                unrolledImge.push_back(px);
            }
        }
    }

    std::vector<std::size_t> dof2vert = dof_to_vertex_map(*Vh);
    int nnn = unrolledImge.size();

    for(int k=0; k<unrolledImge.size(); k++)
    {
        image->vector()->setitem(dolfin::la_index(k),unrolledImge[dof2vert[k]]);
    }

    return image;
}

What it produces: something, but the degrees of freedom are totally mixed up, showing all possible kind of patterns. I know how to do this with only one image (https://fenicsproject.org/qa/9339/fenics-image-processing-because-matrix-expression-poisson/), but even after 100 index permutations I can’t understand how how these degrees of freedom are ordered. Any suggestion?