Tuesday, November 22, 2005

ITK Geodesic Active Contour 3D Example (Java)

So like promised in the last posting the sample source code for a three dimensional Geodesic Active Contour. (Please note that you need the patches from the previous posting for this to work)

But first I have added links to some good books covering the Segmentation Algorithm used in this example and the whole field of either PDE Segementation, Level Sets, or Visualization -- after all Christmas is around the corner. I just got Shapiro's book in the mail (left) and I might post a thorough review later;-)







And here comes the code. Look at p. 537 of the ITK User Guide for more information on the parameters and the algorithms:


private void myITKCalls ( String inputFile, String OutputFile, float propagationScaling,
            double sigma, double alpha, double beta, int seedX, int seedY, int seedZ,
            float initialDistance ) {
    /* # We're going to build the following pipelines:
    # 1. reader -> smoothing -> gradientMagnitude -> sigmoid -> FI
    # 2. fastMarching -> geodesicActiveContour ( FI ) -> thresholder -> writer
    # The output of pipeline 1 is a feature image that is used by the
    # geodesicActiveContour object.  Also see figure 9.18 in the ITK
    # Software Guide.*/



        itkImageFileReaderUS3_Pointer reader = itkImageFileReaderUS3.itkImageFileReaderUS3_New ( ) ;
            itkImageFileWriterUC3_Pointer writer = itkImageFileWriterUC3.itkImageFileWriterUC3_New ( ) ;

            itkCastImageFilterUS3F3_Pointer inputCast = itkCastImageFilterUS3F3.itkCastImageFilterUS3F3_New ( ) ;

            itkCurvatureAnisotropicDiffusionImageFilterF3F3_Pointer smoothing = itkCurvatureAnisotropicDiffusionImageFilterF3F3.itkCurvatureAnisotropicDiffusionImageFilterF3F3_New ( ) ;

            itkRescaleIntensityImageFilterUS3UC3_Pointer outputCast = itkRescaleIntensityImageFilterUS3UC3.itkRescaleIntensityImageFilterUS3UC3_New ( ) ;

            itkGradientMagnitudeRecursiveGaussianImageFilterF3F3_Pointer gradientMagnitude = itkGradientMagnitudeRecursiveGaussianImageFilterF3F3.itkGradientMagnitudeRecursiveGaussianImageFilterF3F3_New ( ) ;

            itkSigmoidImageFilterF3F3_Pointer sigmoid = itkSigmoidImageFilterF3F3.itkSigmoidImageFilterF3F3_New ( ) ;
            
            itkFastMarchingImageFilterF3F3_Pointer fastMarching = itkFastMarchingImageFilterF3F3.itkFastMarchingImageFilterF3F3_New ( ) ;

            itkGeodesicActiveContourLevelSetImageFilterF3F3_Pointer geodesicActiveContour = itkGeodesicActiveContourLevelSetImageFilterF3F3.itkGeodesicActiveContourLevelSetImageFilterF3F3_New ( ) ;

            itkBinaryThresholdImageFilterF3US3_Pointer thresholder = itkBinaryThresholdImageFilterF3US3.itkBinaryThresholdImageFilterF3US3_New ( ) ;

            reader.SetFileName ( inputFile ) ;
            writer.SetFileName ( OutputFile ) ;

            inputCast.SetInput ( reader.GetOutput ( ) ) ;
            smoothing.SetInput ( inputCast.GetOutput ( ) ) ;

            smoothing.SetNumberOfIterations (   5 ) ;
            smoothing.SetTimeStep ( 0.0625 ) ;
            smoothing.SetConductanceParameter ( 3.0 ) ;
            smoothing.Update ( ) ;
            System.out.println ( "Smoothing done" ) ;

            gradientMagnitude.SetInput ( smoothing.GetOutput ( ) ) ;
            gradientMagnitude.SetSigma ( sigma ) ;
            gradientMagnitude.Update ( ) ;
            System.out.println ( "Gradient Magnitude done" ) ;

            sigmoid.SetInput ( gradientMagnitude.GetOutput ( ) ) ;
            sigmoid.SetOutputMinimum ( ( float ) 0.0 ) ;
            sigmoid.SetOutputMaximum ( ( float ) 1.0 ) ;
            sigmoid.SetAlpha ( alpha ) ;
            sigmoid.SetBeta ( beta ) ;
            sigmoid.Update ( ) ;
            System.out.println ( "Sigmoid done" ) ;
            
            //Fast Marching
            itkIndex3 seedPosition = new itkIndex3 ( ) ; //seed
            seedPosition.SetElement ( 0, seedX ) ;
            seedPosition.SetElement ( 1, seedY ) ;
            seedPosition.SetElement ( 2, seedZ ) ;

            float seedValue = -initialDistance;
            itkLevelSetNodeF3 node = new itkLevelSetNodeF3 ( ) ;
            node.SetValue ( seedValue ) ;
            node.SetIndex ( seedPosition ) ;

            itkNodeContainerF3_Pointer seeds = itkNodeContainerF3.itkNodeContainerF3_New ( ) ;
            seeds.Initialize ( ) ;
            seeds.InsertElement ( 0, node ) ;

            fastMarching.SetTrialPoints ( seeds.GetPointer ( ) ) ;
            fastMarching.SetSpeedConstant ( 1.0 ) ;
            fastMarching.SetOutputSize ( reader.GetOutput ( ) .GetBufferedRegion ( ) .GetSize ( ) ) ;
            fastMarching.Update ( ) ;
            System.out.println ( "Fast Marching done" ) ;

            geodesicActiveContour.SetInput ( fastMarching.GetOutput ( ) ) ;
            geodesicActiveContour.SetFeatureImage ( sigmoid.GetOutput ( ) ) ;
            
            geodesicActiveContour.SetPropagationScaling ( propagationScaling ) ;
            geodesicActiveContour.SetCurvatureScaling ( ( float ) 1.0 ) ;
            geodesicActiveContour.SetAdvectionScaling ( ( float ) 1.0 ) ;
            geodesicActiveContour.SetMaximumRMSError ( 0.02 ) ;
            geodesicActiveContour.SetNumberOfIterations ( 800 ) ;
            geodesicActiveContour.Update ( ) ;
            System.out.println ( "Geodesic done" ) ;

            thresholder.SetLowerThreshold ( ( float ) -1000.0 ) ;
            thresholder.SetUpperThreshold ( ( float ) 0.0 ) ;
            thresholder.SetOutsideValue ( 0 ) ;
            thresholder.SetInsideValue ( 65535 ) ;
            
            thresholder.SetInput ( geodesicActiveContour.GetOutput ( ) ) ;
            System.out.println ( "Threshold done" ) ;

            outputCast.SetInput ( thresholder.GetOutput ( ) ) ;
            writer.SetInput ( outputCast.GetOutput ( ) ) ;

            outputCast.SetOutputMinimum (    ( short ) 0   ) ;
            outputCast.SetOutputMaximum ( ( short ) 255 ) ;

            writer.Update ( ) ;

     }

Friday, November 18, 2005

Java/Python ITK LevelSet Problems

As you might have noticed I like programming in Java and especially in ITK's case using ccmake for my own code doesn't really encourages me to use C. Anyway the Wrappings for other languages are really poor -- especially when it comes to level sets.

For the 2D case I ended up writing shell scripts which called the corresponding ITK Examples. That worked really well and I was even able to replace some of the C stuff with the Python examples -- but for some reasons which probably need further investigation the Python code was probably 5 times slower than the C code.

In the 3D case I was essentially forced to write my own code. So I started out with Java and - boom - the ITK people didn't ptovide the proper binding. Their seems to be a hack as describe in the Mailinglist to fix that -- and their is also a distribution which contains patches for better Wrapping support at inra.

Wednesday, November 16, 2005

Google Analytics

Just as an FYI: I have added Google's Analytics code to this blog so I can have some statistics on how you guys use the links on the right. Somehow I have to pay bills;-)

I played around with it in the last few days and it looks really neat and seems to provide some value. However it has like most things Google issues with https pages. For some reasons the lock looks crossed through on Mozilla after I add anything Google. But I guess that's the price to pay;-)

Friday, November 04, 2005

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


}


Tuesday, November 01, 2005

Put nx user's in jail

Because in our peoject we want to give access to our servers almost to anybody who manages to register on our webpage security is a big issue. There are some virtualization solutions like Xen, qemu, or virtual machine/pc but since our machines are only used for this task we implemented a more light weight solution like a chroot environment. We used jail from JMC Research for our little project:

I installed dozens of program with addjailsw to make nx and java (we want people to use java and vtk;-) run. Below you will find the programs from my directory structure:

bin/

awk cat egrep ln more rm sh uname basename cp grep ls mv rmdir sleep vi bash cut hostname mkdir pwd sed touch


/usr/bin

dirname java nxagent nxloadconfig nxprint nxviewer tr expect javac nxclient nxnode nxproxy passwd which head md5sum nxdesktop nxnode-login nxserver ssh xterm id nslookup nxkeygen nxpasswd nxsetup tail

I also copied the whole java directory in /usr/local (including links) and the whole X11R6 directory to /usr/, too. You might need to check if you need all the binaries in /usr/X11R6/bin

I also copied all the NX files wherever they were into to chroot envionment but then deleted in the chroot-nx-home-directory .ssh -- you only need that if you are going to run nx in the jail, too. I am not sure if that is a good idea because then the keys would be theoretically accessible by the users in the jail environment -- unless you would use a second jail for nx...

I have printed the content of the /etc directory below. I just copied the files from the real environment except passwd and groups -- they are maintained by the jail environment's programs. You might also want to edit the sshd directory to not give so many keys away or use more restrictive settings.


/etc

ld.so.cache mtab nxserver resolv.conf shadow termcap X11 hosts localtime nsswitch.conf passwd services ssh vimrc

nxserver and X11 are directories.

I also implemented some hack (see here on the mailinglist) into jail to allow essentially the usage of ssh in the jail. You might be able to get by without the patch if you make nx su to the other users.

This is a rough sketch from what I did -- drop me a note if you have questions.

Note: I had trouble opening an xterm in the jail -- there seem to be issues with pty. You can find some information here but that didn't help me. So I ended up linking the Java application I wanted to publish on the web to xterm, e.g. ln -sf /var/chroot/home/vtk/myJava /var/chroot/usr/bin/xterm -- now moznx (see other entry) instead of starting some console window starts my application, too. Of course you can't give users now no access to an xterm anymore -- so thsi won't work for everybody.