Nifti, Java, and VTK [updated]

We settled on the Nifti dataformat for our project. Nifti is a data format for MRI volume data dreamed up at NIMH (see here for more information on Nifti) and probably the next unifying thing. Most respectable programs in the field support Nifti -- so it's probably a good way to go.

In order to integrate a Nifti-Reader into VTK using Java I downloaded the Java library from Nifti's Sourceforge project. Then I added the niftijlib.jar to my CLASSPATH and the following code converts Nifti to VTK with the trick used in the ImageJ vtk plugin. Essentially we read the Nifti-File and create a temporary file in VTK format. Then we apply the transformations described in nifti1.h (I currently support only method 1 and 2 -- the code for method 3 is untested)This file then gets read with the standard VTK library function:

import java.io.BufferedOutputStream;
import java.io.File;
import java.io.FileOutputStream;
import java.util.Calendar;
import java.util.Date;

import vtk.vtkImageData;

import vtk.vtkImageChangeInformation;
import vtk.vtkImageFlip;
import vtk.vtkImageResample;
import vtk.vtkImageReslice;
import vtk.vtkImageShiftScale;
import vtk.vtkImageShrink3D;
import vtk.vtkImageSource;
import vtk.vtkImageTranslateExtent;
import vtk.vtkMatrix4x4;
import vtk.vtkMatrixToLinearTransform;
import vtk.vtkStructuredPointsReader;

public  class NiftiVTKReader extends vtkImageSource {
 private final int OLD_WAY = 0;
 private final int NORMAL_CASE = 1;
 private final int METHOD_3 = 2;
 
 private vtkImageData image = null;
 private String filename = null;
 
 
 
 private static String uniqueName = null;

 private static String generateUniqueName() {
  if (uniqueName==null) {
   Calendar calendar = Calendar.getInstance();
   uniqueName = "dflImageData"+(new Long(calendar.getTimeInMillis())).toString();
  }
  return uniqueName;
 }
 
 public String GetFileName() {
  // TODO Auto-generated method stub
  return filename;
 }

 
 

 @Override
 public vtkImageData GetOutput() {
  if (image==null) convertNiftiData();
  return image;
 }





 public void SetFileName(String arg0) {
  filename = arg0;
  image = null; //reset image
 }
 
 private static String  CreateVTKHeader(int maxx, int maxy, int maxz, double voxelX, double voxelY, double voxelZ) { 
  String header = "";
  header = "# vtk DataFile Version 2.0\n"
      + " generated with NiftiVTKReader on " + (new Date()).toString() + "\n"
      + "BINARY\n" 
      + "DATASET STRUCTURED_POINTS\n"
      + "DIMENSIONS " + maxx + " " + maxy + " " + maxz + "\n"
      + "ASPECT_RATIO " + voxelX + " " + voxelY + " " + (voxelZ<0.0001?1.5:voxelZ) + "\n"
      + "ORIGIN 0 0 0 \n"
      + "POINT_DATA " + maxx*maxy*maxz + "\n"
      + "SCALARS volume_scalars short 1 \n"
      + "LOOKUP_TABLE default\n"; 
  System.out.println(header);
  return header;
 }
 
 /*---------------------------------------------------------------------------*/
 /*! Given the quaternion parameters (etc.), compute a transformation matrix. 
  * (taken from nifti1_io.c -- Robert W Cox)

    See comments in nifti1.h for details.
      - qb,qc,qd = quaternion parameters
      - qx,qy,qz = offset parameters
      - dx,dy,dz = grid stepsizes (non-negative inputs are set to 1.0)
      - qfac     = sign of dz step (< 0 is negative; >= 0 is positive)

    <pre>
    If qx=qy=qz=0, dx=dy=dz=1, then the output is a rotation matrix.
    For qfac >= 0, the rotation is proper.
    For qfac <  0, the rotation is improper.
    </pre>

    \see "QUATERNION REPRESENTATION OF ROTATION MATRIX" in nifti1.h
    \see nifti_mat44_to_quatern, nifti_make_orthog_mat44,
        nifti_mat44_to_orientation

 *//*-------------------------------------------------------------------------*/
 private vtkMatrix4x4 nifti_quatern_to_mat44( float qb, float qc, float qd,
                               float qx, float qy, float qz,
                               float dx, float dy, float dz, float qfac )
 {
    vtkMatrix4x4 R = new vtkMatrix4x4();
    double a,b=qb,c=qc,d=qd , xd,yd,zd ;

    /* last row is always [ 0 0 0 1 ] */

    R.SetElement(3,0,0);
    R.SetElement(3,1,0);
    R.SetElement(3,2,0);
    R.SetElement(3,3,1.0);

    /* compute a parameter from b,c,d */

    a = 1.0 - (b*b + c*c + d*d) ;
    if( a < 1.e-7 ){                   /* special case */
      a = 1.0 / Math.sqrt(b*b+c*c+d*d) ;
      b *= a ; c *= a ; d *= a ;        /* normalize (b,c,d) vector */
      a = 0.0 ;                        /* a = 0 ==> 180 degree rotation */
    } else{
      a = Math.sqrt(a) ;                     /* angle = 2*arccos(a) */
    }

    /* load rotation matrix, including scaling factors for voxel sizes */

    xd = (dx > 0.0) ? dx : 1.0 ;       /* make sure are positive */
    yd = (dy > 0.0) ? dy : 1.0 ;
    zd = (dz > 0.0) ? dz : 1.0 ;
    
    System.out.println("xd: " + xd + "yd: " + yd + "zd: " + zd);

    if( qfac < 0.0 ) zd = -zd ;         /* left handedness? */

    R.SetElement(0,0, (a*a+b*b-c*c-d*d) * xd);
    R.SetElement(0,1, 2.0 * (b*c-a*d        ) * yd);
    R.SetElement(0,2, 2.0 * (b*d+a*c        ) * zd);
    R.SetElement(1,0, 2.0 * (b*c+a*d        ) * xd);
    R.SetElement(1,1, (a*a+c*c-b*b-d*d) * yd);
    R.SetElement(1,2, 2.0 * (c*d-a*b        ) * zd);
    R.SetElement(2,0, 2.0 * (b*d-a*c        ) * xd);
    R.SetElement(2,1, 2.0 * (c*d+a*b        ) * yd );
    R.SetElement(2,2,  (a*a+d*d-c*c-b*b) * zd);
    

    /* load offsets */

    R.SetElement(0,3, qx);
    R.SetElement(1,3,qy);
    R.SetElement(2,3,qz);

    return R ;
 }

 private double sqr(double a) {
  return a*a;
 }
 /*
  * Conversion of a quaternion into a direction cosine matrix   
  */
 private double[] QuaternionToDirectionCosines(double q2, double q3, double q4) {
  double mat[] = new double[9];
  double q1 = Math.sqrt(1-sqr(q2) + sqr(q3) + sqr(q4));
  
  //x
  mat[0] = sqr(q1) - sqr(q2) - sqr(q3) + sqr(q4);
  mat[1] = 2*(q1*q2 - q3*q4);
  mat[2] = 2*(q1*q3+q2*q4);
  
  //y
  mat[3] = 2*(q1*q2 + q3*q4);
  mat[4] = -sqr(q1) +sqr(q2) - sqr(q3) + sqr(q4);
  mat[5] = 2*(q2*q3 - q1*q4);
  
  //z
  mat[6] = 2*(q1*q3 - q2*q4);
  mat[7] = 2*(q2*q3 + q1*q4);
  mat[8] = -sqr(q1) - sqr(q2) + sqr(q3) + sqr(q4);
  
  System.out.println("Direction Cosines: " + SegementationPane.printDouble(mat));
  
  return mat;
 }

 private double MatrixVectorNorm(double[] matrix, int pos, double value) {
  double result = 0;
  for (int i = pos*3;i<pos*3+3;i++) {
   result = result + sqr(matrix[i]*value);
  }
  return Math.sqrt(result);
  
  /*double vector[] = new double[3];
  for (int i=0; i<3; i++) vector[i] = 0;
  vector[pos] = value;
  
  double result[] = {0,0,0};
  
  for (int i=0; i<3;i++) {
   for (int j=0; j<3; j++) {
    result[j] = result[j] + vector[i]*matrix[j+i];
   }
  }
  double res = 0;
  for (int i=0; i<3; res=res+sqr(result[i++]));
  return Math.sqrt(res);*/
  
 }
 
 /** do everything
  */
 private void convertNiftiData( )
    {
   Nifti1Dataset nds = new Nifti1Dataset(filename);
   try {
    nds.readHeader();
    nds.printHeader();
    
    //Determine 3D Image (Volume) Orientation and Location in Space
    int orientation = NORMAL_CASE;
    if (nds.qform_code == 0) { //The "old" way (ANALYZE 7.5 way)
     orientation = OLD_WAY;
    } else if (nds.qform_code > 0) { //"normal" case
     orientation = NORMAL_CASE;
    } else if (nds.sform_code > 0) { // affine transformation case
        orientation = METHOD_3;
        
       }
    //generate temporary file
    File tmpFile = File.createTempFile(generateUniqueName(), ".vtk");
       String tmpFileName = tmpFile.getAbsolutePath();
    
       //write vtk Data from Nifti
       BufferedOutputStream bos = new BufferedOutputStream(new FileOutputStream(tmpFileName));
       byte [] volume = nds.readData();
       
       String header = "";

       /* Determione header */
       switch (orientation) {
        case OLD_WAY:
         header = CreateVTKHeader(nds.dim[1], nds.dim[2], 
           nds.dim[3], nds.pixdim[1], nds.pixdim[2], nds.pixdim[3]);
         break;
        case NORMAL_CASE: 
          /* For some Voodoo reason you can't chnage the spacing
           * in VTK after you load the image -- so we do it in the
           * header
           */
          double direction_cosines[] = this.QuaternionToDirectionCosines
           (nds.quatern[0], nds.quatern[1], nds.quatern[2]);
         
          header = CreateVTKHeader(nds.dim[1], nds.dim[2], nds.dim[3], 
            MatrixVectorNorm(direction_cosines, 0, nds.pixdim[1]),
            MatrixVectorNorm(direction_cosines, 1, nds.pixdim[2]),  
            MatrixVectorNorm(direction_cosines, 1, nds.qfac*nds.pixdim[3]));
         break;
        case METHOD_3: 
         header = CreateVTKHeader(nds.dim[1], nds.dim[2], 
           nds.dim[3], nds.srow_x[3], nds.srow_y[3], nds.srow_z[3]);
         break;
       }
      
       bos.write(header.getBytes());
       
       //write volume (check endianess!!)
       if (!nds.big_endian) {
        for (int i=1; i<volume.length; i+=2) {
         //switch endianess
         bos.write(volume[i]);
         bos.write(volume[i-1]);
        }
       } else {
        bos.write(volume); //no need to switch
       }
       
       bos.close();
       
       System.out.println("Wrote VTK file");
    
       // Read the temporary file using VTK.
       vtkStructuredPointsReader reader = new vtkStructuredPointsReader();
       reader.SetFileName(tmpFileName);
       reader.Update();
       reader.CloseVTKFile();

      image = (vtkImageData) reader.GetOutput();

       // Remove the temporary file.
       tmpFile.delete();
       
       System.out.println("VTK file deleted");
       
       switch (orientation) {
        case OLD_WAY: break; //Do nothing
        case NORMAL_CASE: { 
         
         vtkMatrix4x4 mat = this.nifti_quatern_to_mat44(
           nds.quatern[0],nds.quatern[1], nds.quatern[2], 
           nds.qoffset[0], nds.qoffset[1], nds.qoffset[2],
           //0,0,0,
           nds.pixdim[1], nds.pixdim[2], nds.pixdim[3],
           //1,1,1,
           nds.qfac);
         //vtkMatrixToLinearTransform trans = new vtkMatrixToLinearTransform();
         //trans.SetInput(mat); 
         
         //bring it into the Cox coordinate system
         vtkImageFlip flip = new vtkImageFlip();
      flip.SetFilteredAxes(2); //flip z-axes
      flip.SetInput(reader.GetOutput());
         
         vtkImageChangeInformation info = new vtkImageChangeInformation();
         info.SetInput(reader.GetOutput());
         info.SetOutputOrigin(0,0,0 );
         //info.SetExtentTranslation(-image.GetDimensions()[0]/2,0, 0);
         //info.CenterImageOn();
         System.out.println("Extend 1: " + SegementationPane.printExtend(image.GetExtent()));
         
         vtkImageReslice reslice = new vtkImageReslice();
         reslice.SetInput(info.GetOutput());
         //reslice.SetResliceAxes(mat);
         reslice.SetResliceAxesOrigin(QuaternionToDirectionCosines
           (nds.quatern[0], nds.quatern[1], nds.quatern[2]));
         //reslice.SetResliceTransform(trans);
         reslice.SetOutputSpacing(nds.pixdim[1],nds.pixdim[2], nds.pixdim[3]);
         
         vtkImageChangeInformation info2 = new vtkImageChangeInformation();
         info2.SetInput(reslice.GetOutput());
         info2.SetOutputSpacing(nds.pixdim[1], nds.pixdim[2], nds.pixdim[3]);
         //reslice.WrapOn();
         //reslice.MirrorOn();
         //reslice.Update();
      
        /* vtkImageShrink3D shrink = new vtkImageShrink3D();
         shrink.SetShrinkFactors(Math.round(nds.pixdim[1]*10),Math.round(nds.pixdim[2]*10), Math.round(nds.pixdim[3])*10);
         shrink.SetInput(reslice.GetOutput());
         
         vtkImageShiftScale scale = new vtkImageShiftScale();
         scale.SetScale(1);
         scale.SetInput(shrink.GetOutput());*/
         
         
         
         //transle
         /*vtkImageTranslateExtent translate = new vtkImageTranslateExtent();
         translate.SetInput(resample.GetOutput());
         translate.SetTranslation(Math.round(nds.qoffset[0]), Math.round(nds.qoffset[1]), Math.round(nds.qoffset[2]));
         */
        
         //reslice.DebugOn();
         System.out.println("Before Reslice");
        
         //scale.Update();
         info.Update();
         
         image = info.GetOutput();
         //image = scale.GetOutput(); 
         System.out.println("Extend 2: " + SegementationPane.printExtend(image.GetDimensions()));
         System.out.println("spacing: " + SegementationPane.printDouble(image.GetSpacing()));
         System.out.println("After Reslice");          
         break;
        }
        case METHOD_3: //do the affine transformation
         System.out.println("Starting reslice!");
         vtkMatrix4x4 mat = new vtkMatrix4x4();
         for (int j=0;j<3; j++) mat.SetElement(0,j, nds.srow_x[j]);
         for (int j=0;j<3; j++) mat.SetElement(1,j, nds.srow_y[j]);
         for (int j=0;j<3; j++) mat.SetElement(2,j, nds.srow_z[j]);
         
         vtkMatrixToLinearTransform trans = new vtkMatrixToLinearTransform();
         trans.SetInput(mat);
         trans.Print();
         
         vtkImageReslice reslice = new vtkImageReslice();
         reslice.SetInput(image);
         reslice.SetResliceTransform(trans);
         reslice.Update();
         image = reslice.GetOutput();
         System.out.println("Reslice done!"); 
         break;
       }
            } catch (Exception e) { 
             System.err.println("Error!!");
             System.err.print(e.getMessage());
             image = null;
             };
         
    } // end method
  

}