Click here to Skip to main content
13,091,920 members (58,903 online)
Rate this:
Please Sign up or sign in to 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
if (samplesPerPixel == 1 &&  bitsAllocated == 8)
    pixel_8=( byte*)malloc(numPixels);

    for(int f=0;f<number_of_frames;f++)
        for(int i=0;i<height;i++)
            for(int j=0;j<width;j++)
                int il=(i*width)+j;
                MCPoint *vert=&mcPoints[f][i][j]; 
                vert->val = unsig;  

MCPoint has the follow stucture
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:
MarchingCubeObject = new CMarchingCubes(width, height, number_of_frames, 200, mcPoints);
gl_widget=new CGlWidget(NULL,MarchingCubeObject->getNumTraingle(),MCTriangles,0);

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
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] - 
        }   //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;


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);
        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
typedef struct {
mpVector p[3];
mpVector norm;

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 ) :

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++)                            

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 19-Jan-13 0:48am

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

  Print Answers RSS
Top Experts
Last 24hrsThis month

Advertise | Privacy |
Web03 | 2.8.170813.1 | Last Updated 19 Jan 2013
Copyright © CodeProject, 1999-2017
All Rights Reserved. Terms of Service
Layout: fixed | fluid

CodeProject, 503-250 Ferrand Drive Toronto Ontario, M3C 3G8 Canada +1 416-849-8900 x 100