[vtkusers] vtk dicom reslice axial, coronal and sagittal view

[p.l.]

I have a problem with calculating the correct axial, coronal and sagittal
images from a given DICOM series. I use the vtk-example in:
VTK\Examples\ImageProcessing\Cxx\ImageSlicing.cxx : 

     vtkSmartPointer<vtkDICOMReader> dicomReader =     
vtkSmartPointer<vtkDICOMReader>::New(); 

dicomReader->SetFileNames(sortedFileNames); 
dicomReader->SetDataScalarTypeToUnsignedChar(); 
dicomReader->AutoRescaleOn(); 

dicomReader->SetGlobalWarningDisplay(1); 
dicomReader->Update(); 

int* dicomExtent = new int[6]; 
dicomReader->GetDataExtent(dicomExtent); 

double* dicomOrigin = new double[3]; 
dicomReader->GetDataOrigin(dicomOrigin); 

double* dicomSpacing = new double[3]; 
dicomReader->GetDataSpacing(dicomSpacing); 

double dicomCenter[3];   
dicomCenter[0] = dicomOrigin[0] + dicomSpacing[0] * 0.5 * (dicomExtent[0] +
dicomExtent[1]); 
dicomCenter[1] = dicomOrigin[1] + dicomSpacing[1] * 0.5 * (dicomExtent[2] +
dicomExtent[3]); 
dicomCenter[2] = dicomOrigin[2] + dicomSpacing[2] * 0.5 * (dicomExtent[4] +
dicomExtent[5]); 

// Matrices for axial, coronal, sagittal, oblique view orientations 
static double axialElements[16] = { 
    1, 0, 0, 0, 
    0, 1, 0, 0, 
    0, 0, 1, 0, 
    0, 0, 0, 1 }; 

static double coronalElements[16] = { 
    1, 0, 0, 0, 
    0, 0, 1, 0, 
    0, -1, 0, 0, 
    0, 0, 0, 1 }; 

    static double sagittalElements[16] = { 
    0, 0, -1, 0, 
    1, 0, 0, 0, 
    0, -1, 0, 0, 
    0, 0, 0, 1 }; 

//static double obliqueElements[16] = { 
//         1, 0, 0, 0, 
//         0, 0.866025, -0.5, 0, 
//         0, 0.5, 0.866025, 0, 
//         0, 0, 0, 1 }; 

// Set the slice orientation 
vtkSmartPointer<vtkMatrix4x4> resliceAxes = 
    vtkSmartPointer<vtkMatrix4x4>::New(); 

resliceAxes->DeepCopy(axialElements); 

// Set the point through which to slice 
resliceAxes->SetElement(0, 3, dicomCenter[0]); 
resliceAxes->SetElement(1, 3, dicomCenter[1]); 
resliceAxes->SetElement(2, 3, dicomCenter[2]); 

// Extract a slice in the desired orientation 
vtkSmartPointer<vtkImageReslice> reslice = 
    vtkSmartPointer<vtkImageReslice>::New(); 
reslice->SetInputConnection(dicomReader->GetOutputPort()); 
reslice->SetOutputDimensionality(2); 
reslice->SetResliceAxes(resliceAxes); 
reslice->SetInterpolationModeToLinear(); 
reslice->Update(); 


The problem is, that the matrices for the axial, coronal and sagittal
orientation only work if the direction cosines form an identity matrix. But
I have lots of DICOM series with arbitrary direction cosines, like:
0.758921/ 0.498738/ 0.418688/ -0.00570692/ 0.648034/ -0.76159 So I thought
that I should use the direction cosines in the axial, coronal and sagittal
matrix; this is my code: 

    vtkSmartPointer<vtkDICOMReader> dicomReader =
vtkSmartPointer<vtkDICOMReader>::New(); 

dicomReader->SetFileNames(sortedFileNames); 
dicomReader->SetDataScalarTypeToUnsignedChar(); 
dicomReader->AutoRescaleOn(); 

dicomReader->SetGlobalWarningDisplay(1); 
dicomReader->Update(); 

int* dicomExtent = new int[6]; 
dicomReader->GetDataExtent(dicomExtent); 

double* dicomOrigin = new double[3]; 
dicomReader->GetDataOrigin(dicomOrigin); 

double* dicomSpacing = new double[3]; 
dicomReader->GetDataSpacing(dicomSpacing); 

double dicomCenter[3];   
dicomCenter[0] = dicomOrigin[0] + dicomSpacing[0] * 0.5 * (dicomExtent[0] +
dicomExtent[1]); 
dicomCenter[1] = dicomOrigin[1] + dicomSpacing[1] * 0.5 * (dicomExtent[2] +
dicomExtent[3]); 
dicomCenter[2] = dicomOrigin[2] + dicomSpacing[2] * 0.5 * (dicomExtent[4] +
dicomExtent[5]); 

double* directionCosine = new double[6]; 
dicomReader->GetMedicalImageProperties()->GetDirectionCosine(directionCosine); 
double* xAxis = new double[3]; 
xAxis[0] = directionCosine[0]; 
xAxis[1] = directionCosine[1]; 
xAxis[2] = directionCosine[2]; 

double* yAxis = new double[3]; 
yAxis[0] = directionCosine[3]; 
yAxis[1] = directionCosine[4]; 
yAxis[2] = directionCosine[5]; 

vtkSmartPointer<vtkMath> math = vtkSmartPointer<vtkMath>::New(); 

double* zAxis = new double[3]; 
math->Cross(xAxis, yAxis, zAxis); 

// Matrices for axial, coronal, sagittal, oblique view orientations 
//xy 
static double axialElements[16] = { 
    xAxis[0], yAxis[0], zAxis[0], 0, 
    xAxis[1], yAxis[1], zAxis[1], 0, 
    xAxis[2], yAxis[2], zAxis[2], 0, 
    0, 0, 0, 1 }; 


double* yCross = new double[3]; 
math->Cross(xAxis, zAxis, yCross); 
//xz 
static double coronalElements[16] = { 
    xAxis[0], zAxis[0], yCross[0], 0, 
    xAxis[1], zAxis[1], yCross[1], 0, 
    xAxis[2], zAxis[2], yCross[2], 0, 
    0, 0, 0, 1 }; 

double* xCross = new double[3]; 
math->Cross(yAxis, zAxis, xCross); 
//yz 
static double sagittalElements[16] = { 
    yAxis[0], zAxis[0], xCross[0], 0, 
    yAxis[1], zAxis[1], xCross[1], 0, 
    yAxis[2], zAxis[2], xCross[2], 0, 
    0, 0, 0, 1 }; 

//static double obliqueElements[16] = { 
//         1, 0, 0, 0, 
//         0, 0.866025, -0.5, 0, 
//         0, 0.5, 0.866025, 0, 
//         0, 0, 0, 1 }; 

// Set the slice orientation 
vtkSmartPointer<vtkMatrix4x4> resliceAxes = 
    vtkSmartPointer<vtkMatrix4x4>::New(); 

resliceAxes->DeepCopy(axialElements); 

// Set the point through which to slice 
resliceAxes->SetElement(0, 3, dicomCenter[0]); 
resliceAxes->SetElement(1, 3, dicomCenter[1]); 
resliceAxes->SetElement(2, 3, dicomCenter[2]); 

// Extract a slice in the desired orientation 
vtkSmartPointer<vtkImageReslice> reslice = 
    vtkSmartPointer<vtkImageReslice>::New(); 
reslice->SetInputConnection(dicomReader->GetOutputPort()); 
reslice->SetOutputDimensionality(2); 
reslice->SetResliceAxes(resliceAxes); 
reslice->SetInterpolationModeToLinear(); 
reslice->Update(); 

But it does not work. I also tried to use the inverse of the direction
cosines, but the result is not correct. Can someone help me with my problem?
Thank you in advance! 

Kind regards, Pia 



--
View this message in context: http://vtk.1045678.n5.nabble.com/vtk-dicom-reslice-axial-coronal-and-sagittal-view-tp5736244.html
Sent from the VTK - Users mailing list archive at Nabble.com.