Click here to Skip to main content
15,885,216 members
Please Sign up or sign in to vote.
4.00/5 (1 vote)
See more:
Hey Guys. Well, i've been working with DICOM for a while and now i've some problems ,with displaying MULTIFRAME CT images .it's my first post here on SO ,so i'm sorry if it's not that good :'( i'm using marching cubes (of course) to construct the surface and then opengl (under Qt 4.8.1) to display it . This piece of code is what i'm using to read the frames one by one .i have a 8bits width,height,number_of_frames i've got them from the header of the dicom file
C++
if (samplesPerPixel == 1 &&  bitsAllocated == 8)
{

    pixel_8=( byte*)malloc(numPixels);


    for(int f=0;f<number_of_frames;f++)
    {
        fread(pixel_8,1,numPixels,dicom_file);
        for(int i=0;i<height;i++)
        {
            for(int j=0;j<width;j++)
            {
                int il=(i*width)+j;
                unsig=pixel_8[il];
                buff[j*3]=unsig;
                buff[(j*3)+1]=unsig;
                buff[(j*3)+2]=unsig;
                MCPoint *vert=&mcPoints[f][i][j]; 

                vert->x=f;
                vert->y=i;
                vert->z=j;
                vert->val = unsig;  
            }
        }
    }
}

MCPoint has the follow stucture
C++
struct MCPoint { 
float x, y, z; 
float val; 
 }; 


i've also tried to save the extracted frames one as a .jpg image and it works i got them . but when i try to send the generated matrice pixel_8 to the marching cube functions i got something weird (i will add a picture bellow)

After reading the frames i'm calling marching cube with the following code:
C++
MarchingCubeObject = new CMarchingCubes(width, height, number_of_frames, 200, mcPoints);
MCTriangles=MarchingCubeObject->getTriangles();
gl_widget=new CGlWidget(NULL,MarchingCubeObject->getNumTraingle(),MCTriangles,0);
gl_widget->show();

and here is the Marching cube Class ( i don't think that there is something wrong on the code of MC because i've got it from a trusted place :) ) i have changed only the Step 3
C++
CMarchingCubes::CMarchingCubes(int ncellsX, int ncellsY, int ncellsZ, float minValue, MCPoint *** points)
{
TRIANGLE * triangles = new TRIANGLE[3*ncellsX*ncellsY*ncellsZ];
numTriangles = int(0);

int YtimeZ = (ncellsY+1)*(ncellsZ+1);
//go through all the points
for(int k=0; k < ncellsZ; k++)          //x axis
    for(int j=0; j < ncellsY; j++)      //y axis
        for(int i=0; i < ncellsX; i++)  //z axis
        {
            //initialize vertices
            mp4Vector verts[8];

            int ind = i*YtimeZ + j*(ncellsZ+1) + k;
            /*(step 3)*/ 
            verts[0]=construct_vert_from_pt(points[k  ][j  ][i  ]); 
            verts[1] = construct_vert_from_pt(points[k][j  ][i+1  ]);
            verts[2] = construct_vert_from_pt(points[k+1][j  ][i+1]); 
            verts[3] =  construct_vert_from_pt(points[k+1  ][j  ][i+1]); 
            verts[4] = construct_vert_from_pt(points[k  ][j+1][i  ]); 
            verts[5] = construct_vert_from_pt(points[k][j+1][i+1  ]); 
            verts[6] =construct_vert_from_pt(points[k+1][j+1][i+1]);   
            verts[7] =construct_vert_from_pt(points[k+1  ][j+1][i]);


            //get the index
            int cubeIndex = int(0);
            for(int n=0; n < 8; n++) /*(step 4)*/       
                if(verts[n].val >= minValue) cubeIndex |= (1 << n);

            //check if its completely inside or outside
            /*(step 5)*/
            if(!edgeTable[cubeIndex]) continue;

            //get intersection vertices on edges and save into the array
            mpVector intVerts[12];
            /*(step 6)*/ 
            if(edgeTable[cubeIndex] & 1) intVerts[0] = LinearInterp(verts[0], verts[1], minValue);
            if(edgeTable[cubeIndex] & 2) intVerts[1] = LinearInterp(verts[1], verts[2], minValue);
            if(edgeTable[cubeIndex] & 4) intVerts[2] = LinearInterp(verts[2], verts[3], minValue);
            if(edgeTable[cubeIndex] & 8) intVerts[3] = LinearInterp(verts[3], verts[0], minValue);
            if(edgeTable[cubeIndex] & 16) intVerts[4] = LinearInterp(verts[4], verts[5], minValue);
            if(edgeTable[cubeIndex] & 32) intVerts[5] = LinearInterp(verts[5], verts[6], minValue);
            if(edgeTable[cubeIndex] & 64) intVerts[6] = LinearInterp(verts[6], verts[7], minValue);
            if(edgeTable[cubeIndex] & 128) intVerts[7] = LinearInterp(verts[7], verts[4], minValue);
            if(edgeTable[cubeIndex] & 256) intVerts[8] = LinearInterp(verts[0], verts[4], minValue);
            if(edgeTable[cubeIndex] & 512) intVerts[9] = LinearInterp(verts[1], verts[5], minValue);
            if(edgeTable[cubeIndex] & 1024) intVerts[10] = LinearInterp(verts[2], verts[6], minValue);
            if(edgeTable[cubeIndex] & 2048) intVerts[11] = LinearInterp(verts[3], verts[7], minValue);

            //now build the triangles using triTable
            for (int n=0; triTable[cubeIndex][n] != -1; n+=3) {
                /*(step 7)*/    
                triangles[numTriangles].p[0] = intVerts[triTable[cubeIndex][n+2]];
                triangles[numTriangles].p[1] = intVerts[triTable[cubeIndex][n+1]];
                triangles[numTriangles].p[2] = intVerts[triTable[cubeIndex][n]];
                /*(step 8)*/    
                    triangles[numTriangles].norm = ((triangles[numTriangles].p[1] - 
                    triangles[numTriangles].p[0]).Cross(triangles[numTriangles].p[2] - 
                    triangles[numTriangles].p[0])).Normalize();
                numTriangles++;
            }

        }   //END OF FOR LOOP

        //free all the wasted space
        retTriangles = new TRIANGLE[numTriangles];
        for(int i=0; i < numTriangles; i++)
            retTriangles[i] = triangles[i];
        delete [] triangles;


}


CMarchingCubes::~CMarchingCubes(void)
{
}

mpVector CMarchingCubes::LinearInterp( mp4Vector p1, mp4Vector p2, float value )
{
    mpVector p;
    if(p1.val != p2.val)
        p = (mpVector)p1 + ((mpVector)p2 - (mpVector)p1)/(p2.val - p1.val)*(value - p1.val);
    else 
        p = (mpVector)p1;
    return p;
}

TRIANGLE * CMarchingCubes::getTriangles()
{
    return retTriangles;
}

mp4Vector CMarchingCubes::construct_vert_from_pt( MCPoint p )
{
    return mp4Vector(p.x,p.y,p.z,p.val );
}

int CMarchingCubes::getNumTraingle()
{
    return numTriangles;
}

here is also the TRIANGLE structure
C++
typedef struct {
mpVector p[3];
mpVector norm;
} TRIANGLE;


Basically i'm reading the matrice from dicom file and then send it to Marching cube as i said , there is anything the i have to change on the matrice before send it or what ?

After that i use this function to display what i got (i've divided by 200 because with actual values i've got an empty black screen ,so i tried to scale them ) :

C++
 glViewport(0,0,700,700);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

glLoadIdentity();
glBegin(GL_TRIANGLES);

for(int i=0; i < numOfTriangles; i++){
    glNormal3f(MCTriangles[i].norm.x, MCTriangles[i].norm.y, MCTriangles[i].norm.z);
    for(int j=0; j < 3; j++)                            
        glVertex3f(MCTriangles[i].p[j].x/200,MCTriangles[i].p[j].y/200,MCTriangles[i].p[j].z/200);
}

glEnd();


This is what i got when trying to display the result triangles [ Marching Cube result Picture]
And these are the 16 frames that i have on the DICOM file after saving them as JPG :
[16 Extracted JPG from the DICOM file]
Thank you again for your help (waiting for your commments guys :'( ).
Best Regards.
Posted

This content, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)



CodeProject, 20 Bay Street, 11th Floor Toronto, Ontario, Canada M5J 2N8 +1 (416) 849-8900