Click here to Skip to main content
15,898,010 members

Undefined result when Displaying Multiframe DICOM Image With Marching Cubes and OPENGL

AziSckii GVb asked:

Open original thread
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.
Tags: C++, 3D

Plain Text
ASM
ASP
ASP.NET
BASIC
BAT
C#
C++
COBOL
CoffeeScript
CSS
Dart
dbase
F#
FORTRAN
HTML
Java
Javascript
Kotlin
Lua
MIDL
MSIL
ObjectiveC
Pascal
PERL
PHP
PowerShell
Python
Razor
Ruby
Scala
Shell
SLN
SQL
Swift
T4
Terminal
TypeScript
VB
VBScript
XML
YAML

Preview



When answering a question please:
  1. Read the question carefully.
  2. Understand that English isn't everyone's first language so be lenient of bad spelling and grammar.
  3. If a question is poorly phrased then either ask for clarification, ignore it, or edit the question and fix the problem. Insults are not welcome.
  4. Don't tell someone to read the manual. Chances are they have and don't get it. Provide an answer or move on to the next question.
Let's work to help developers, not make them feel stupid.
Please note that all posts will be submitted under the http://www.codeproject.com/info/cpol10.aspx.



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