[vtkusers] new post

[Pia Langoth]

vtk dicom reslice axial, coronal and sagittal view<http://stackoverflow.com/questions/35083162/vtk-dicom-reslice-axial-coronal-and-sagittal-view>

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




-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/vtkusers/attachments/20160129/548cf886/attachment-0001.html>