Ну, я работал с DICOM некоторое время, и теперь у меня есть некоторые проблемы с отображением изображений MULTIFRAME CT. Я использую марширующие кубы (конечно), чтобы построить поверхность, а затем opengl (под Qt 4.8.1), чтобы отобразить ее.
Этот фрагмент кода — то, что я использую, чтобы прочитать кадры один за другим. У меня есть 8 бит * ширина, высота, число_кадров *, я получил их из заголовка файла dicom
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 имеет следующую структуру
struct MCPoint {
float x, y, z;
float val;
};
Я также пытался сохранить извлеченные кадры один как .jpg изображение, и это работает, я получил их. но когда я пытаюсь отправить сгенерированную матрицу pixel_8 к функциям движущегося куба у меня появилось что-то странное (ниже я добавлю картинку)
После прочтения фреймов я вызываю марширующий куб со следующим кодом:
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();
а вот класс Marching cube Class (я не думаю, что в коде MC что-то не так, потому что я получил его из надежного места :)) я изменил только Шаг 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] -
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;
}
здесь также структура ТРЕУГОЛЬНИКА
typedef struct {
mpVector p[3];
mpVector norm;
} TRIANGLE;
По сути, я читаю матрицу из файла dicom, а затем отправляю ее в марширующий куб, как я уже сказал, есть что-то, что я должен изменить на матрице перед отправкой или что ???
После этого я использую эту функцию, чтобы отобразить то, что я получил (я разделил на 200, потому что с фактическими значениями у меня пустой черный экран, поэтому я попытался их масштабировать):
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();
Вот что я получил при попытке отобразить треугольники результата:
Ссылка на изображение:
Марширующий кубик результат
И это 16 кадров, которые у меня есть в файле DICOM после сохранения их в формате JPG:
Ссылка на изображение: 16 извлеченных изображений jpg из файла dicom
Задача ещё не решена.
Других решений пока нет …