levelset.java

LevelSet algorithme written in java language

import java.io.*;
import java.util.*;
import java.lang.*;

/***
 * This class defines a Point (pixel) class , and some basic operations on a Point object  
 ***/

class Point {
    private int x;  // column
    private int y;  // row
    
    public Point(int x, int y) {
		this.x = x;
		this.y = y;
    }
    
    public int getX() {
		return x;
    }

    public int getY() {
		return y;
    }
    
    public boolean equals(Point p) {
		return (x == p.getX() && y == p.getY());
    }
    
    public double distance(Point p) {
		return Math.sqrt((x - p.getX()) * (x - p.getX()) + (y - p.getY()) * (y - p.getY()));
    }
}



/***
 * This class is the main class. It implements the fast marching and level set methods  
 ***/
public class LevelSet {
    
    // auxilarry constants used for selecting max and min values of certain variables
    private static final double BIG_NUMBER = Double.MAX_VALUE;
    private static final double SMALL_NUMBER = Integer.MIN_VALUE;
    
    // image matrix dimension, their values are read in from image files
    private static int largeur = 0;
    private static int hauteur = 0;
    
    // computation step length (unit is 1)
    private static final int DELTA_X = 1;
    private static final int DELTA_Y = 1;
    
    // initial points (seeds) for fast marching, there might be multiple points
    // location depends on the input image, and must be specified by user
	// NEED TO MODIFY FOR DIFFERENT IMAGES!!!
    private static final int OI1 = 66;   // row
    private static final int OJ1 = 87;   // colmun
    
    private static final int OI2 = 122;  // row
    private static final int OJ2 = 87;   // column
	// NEED TO MODIFY FOR DIFFERENT IMAGES!!!
    
    // control parameters for fast marching
    private static final int FMSteps = 11000;            // steps of fast marching 
    private static final double ALPHA = 70.0;           // multiplier for exponential term 
    
    // control parameters for level set
    private static final int iter_inner = 5;            // steps of iterations for inner loop 
    private static final int iter_outer = 0;            // steps of iterations for outer loop
    private static final double DELTA_T = 0.0005;       // time step
    private static final double bandWidth = 3;          // narrow band width
    private static final double FA = 1;                 // advection force term
    private static final double EPSILON = 0.025;        // curvature force term
    private static final double BETA = 0.01;            // attraction force term

    // computation matrix and vectors for fast marching
    private double [][] T = null;                       // T values
    private Vector alive = new Vector();                // alive set
    private Vector close = new Vector();                // close set
    private Vector farAway = new Vector();              // farAway set
    private Vector neighbour = new Vector();            // neighbour set

    // computation matrix and vectors for level set
    private double [][] phi = null;                     // level set Phi value
    private double [][] phiOld = null;                  // old level set Phi value for re-initilaizing nonfront points
    private double [][] Io = null;                      // original image pixel value
    private double [][] I = null;                       // stretched and blurred image pixel value
    private double [][] KI = null;                      // KI gradient magnitude values
    private Vector frontLine = new Vector();            // zero level set (image contour)
    private Vector narrowBand = new Vector();           // narrow band
    private Vector nbBoundary = new Vector();           // NB boundary pixels set
    private Vector testNB = new Vector();               // old phi values for pixels in NB boundary 
   
private int lireEntier(StringTokenizer st) 
    {
	return Integer.parseInt(st.nextToken());
    }

    /***
     * read image file 
     ***/
    private void readFile(String filename) throws FileNotFoundException, IOException {
    BufferedReader in = new BufferedReader(new FileReader(filename));
	String ligne;

	// lecture du magic number
	ligne = in.readLine();
	if(!ligne.equals("P2")) 
        {
	    System.err.println("ce n'est pas un fichier pgm ascii");
	    System.exit(1);
	}

	// lecture des commentaires
	do 
        {
	    ligne = in.readLine();
	} 
        while(ligne.charAt(0) == '#');

	// lecture de largeur et hauteur
	StringTokenizer tokenizer = new StringTokenizer(ligne);
	 largeur = lireEntier(tokenizer);
	hauteur = lireEntier(tokenizer);
		// initialize matrix for fast marching
		T = new double[largeur][hauteur];   
	
		// initilaize matrix for level set
		phi = new double[largeur][hauteur];   
		phiOld = new double[largeur][hauteur]; 
		I = new double[largeur][hauteur]; 
		Io = new double[largeur][hauteur];
		KI = new double[largeur][hauteur]; 

	// lecture du ng maximum
	ligne = in.readLine();
	tokenizer = new StringTokenizer(ligne);
	int ngMax = lireEntier(tokenizer);

	// lecture des valeurs 
	int x = 0, y = 0;
	while((ligne = in.readLine()) != null) 
        {
	    tokenizer = new StringTokenizer(ligne);
	    while(tokenizer.hasMoreTokens()) 
            {
	        Io[x][y] = lireEntier(tokenizer);
	     	x ++;
		if(x == largeur) 
                {
		    x = 0;
		    y++;
		}
	    }
	}	

    	
    	
    	
    	
    	
    	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

	
	 
	
	
		// read in original image pixel values			
	/*	String line = null;
		int pixelCount = 0;
		while((line = br.readLine()) != null) {
			try{
				StringTokenizer tokens = new StringTokenizer(line.trim());
				int tokenCount = tokens.countTokens();
		
				//check number of pixles in input file
				if((pixelCount + tokenCount) > largeur * hauteur) {
					System.out.println("Input file: incorrect number of pixels.");
					System.exit(1);
				}
		
				// read pixel value
				for (int index = 0; index < tokenCount; index++) {
					int i = pixelCount / largeur;
					int j = pixelCount % largeur;
					Io[i][j] = Double.parseDouble(tokens.nextToken().trim());
					pixelCount++;
				}
		
			} catch(NumberFormatException e) {
				System.out.println("Input file: incorrect pixel value.");
				System.exit(1);
			} 
		}*/
	
		// check if correct number of pixel values read in
		/*if(pixelCount != largeur * hauteur) {
			System.out.println("Input file: incorrect number of pixels.");
			System.exit(1);
		}*/
	
		// blur input image using gaussian smoothing using (5 * 5) matrix
		// boundaries not blurred
		for (int j = 0; j < hauteur; j++) {
			for (int i = 0; i < largeur; i++) {
				if (j >= 2 && j <= (hauteur - 3) && i >= 2 && i <= (largeur - 3)) {
					I[i][j] = (1.0 / 256.0) * 
						(Io[i - 2][j - 2] + 4 * Io[i - 2][j - 1] + 6 * Io[i - 2][j] + 4 * Io[i - 2][j + 1] + Io[i - 2][j + 2]
						 + 4 * Io[i - 1][j - 2] + 16 * Io[i - 1][ j - 1] + 24 * Io[i - 1][j] + 16 * Io[i - 1][j + 1] + 4 * Io[i - 1][j + 2]
						 + 6 * Io[i][j - 2] + 24 * Io[i][j - 1] + 36 * Io[i][j] + 24 * Io[i][j + 1] + 6 * Io[i][j + 2]
						 + 4 * Io[i + 1][j - 2] + 16 * Io[i + 1][ j - 1] + 24 * Io[i + 1][j] + 16 * Io[i + 1][j + 1] + 4 * Io[i + 1][j + 2]
						 + Io[i + 2][j - 2] + 4 * Io[i + 2][j - 1] + 6 * Io[i + 2][j] + 4 * Io[i + 2][j + 1] + Io[i + 2][j + 2]);	
				} else I[i][j] = Io[i][j];	
			}
		}
	
		// find the minimum and maximum pixel values
		double minPixel = BIG_NUMBER;
		double maxPixel = SMALL_NUMBER;
	
		for (int j = 0; j < hauteur; j++) {
			for (int i = 0; i < largeur; i++) {
				if (minPixel > I[i][j]) minPixel = I[i][j];
				if (maxPixel < I[i][j]) maxPixel = I[i][j];
			}
		}
	
		// normalize image (stretch between 0 and 255)
		for (int j = 0; j < hauteur; j++) {
			for (int i = 0; i < largeur; i++) {
				I[i][j] = (I[i][j] - minPixel) * 255.0 / (maxPixel - minPixel);
			}
		}
	
		// write blurred image to an output file for checking
		// comment out for now
		/*
		  PrintWriter pw = new PrintWriter(new BufferedWriter(new FileWriter(new File(filename + "_blur"))));
	
		  // write first three lines to output file
		  pw.println("P2");
		  pw.println(largeur + " " + hauteur);
		  pw.println("255");
	
		  // write blurred pixel values to output file 
		  for (int i = 0; i < hauteur; i++) {
		  for (int j = 0; j < largeur; j++) {
		  pw.print((int) I[i][j] + " ");
		  }
		  pw.println();	    
		  }
		  pw.close();
		*/	

		// compute gradient magnitude (KI) values
		// KI values at boundaries initialized to 0
		for (int j = 0; j < hauteur; j++) {
			for (int i = 0; i < largeur; i++) {
				if (j >= 2 && j <= hauteur - 3 && i >= 2 && i <= largeur - 3) {
					double Ix = I[i][j] - I[i][j - 1];
					double Iy = I[i][j] - I[i - 1][j];
					KI[i][j] = Math.sqrt(Ix * Ix + Iy * Iy);
				} else KI[i][j] = 0;	
			}
		}
              
		// find the minimum and maximum KI value
		double minGra = BIG_NUMBER;
		double maxGra = SMALL_NUMBER;
	
		for (int j = 0; j < hauteur; j++) {
			for (int i = 0; i< largeur; i++) {
				if (minGra > KI[i][j]) minGra = KI[i][j];
				if (maxGra < KI[i][j]) maxGra = KI[i][j];
			}
		}

		// assign boundary KI values to be max KI
		// to avoid "out of boundary" problem in fast marching
		for (int j = 0; j < hauteur; j++) {
			for (int i = 0; i < largeur; i++) {
				if (j >= 2 && j <= (hauteur - 3) && j >= 2 && j <= (largeur - 3)) ;
				else KI[i][j] = maxGra;	
			}
		}

		// normalize KI values (stretch between 0 and 1)
		// large KI values cause errors when updating T values in fast marching
		for (int j = 0; j < hauteur; j++) {
			for (int i = 0; i < largeur; i++) {
				KI[i][j] = (KI[i][j] - minGra) * 1.0 / (maxGra - minGra);
			}
		}
      
		// write KI values to an output file for checking
		// comment out for now
		/*
		  PrintWriter pwk = new PrintWriter(new BufferedWriter(new FileWriter(new File(filename+ "_ki"))));
	
		  // write first three lines to output file
		  pwk.println("P2");
		  pwk.println(largeur + " " + hauteur);
		  pwk.println("255");
	
		  // write KI values to output file 
		  for (int i = 0; i < hauteur; i++) {
		  for (int j = 0; j < largeur; j++) {
		  pwk.print((int) (KI[i][j] * 255) + " ");  // stretched between 0 ~ 255 for display
		  }
		  pwk.println();	    
		  }
		  pwk.close();
		*/
    }
    
    
    /***
     * after fast marching, KI should be stretched back to 0 ~ 255 for level set computation  
     ***/
    private void restoreKi() {
		for (int j = 0; j < hauteur; j++) {
			for (int i = 0; i < largeur; i++) {
				KI[i][j] *= 255.0;
			}
		}
    }
    

    //*****************************************************************************
    // Fast marching step
    //*****************************************************************************

    /***
     * initialize variables for fast marching   
     ***/
		private void initializeFastMarch() {

			// set all points to farAway
			for (int j = 0; j < hauteur; j++) {
				for (int i = 0; i < largeur; i++) {
					farAway.addElement(new Point(j, i));
					T[i][j] = Integer.MAX_VALUE;
				}
			}

			// initialize seed points
			// NEED TO MODIFY FOR DIFFERENT IMAGES!!!
			initialFMAux(OJ1, OI1);
			initialFMAux(OJ2, OI2);	
			// NEED TO MODIFY FOR DIFFERENT IMAGES!!!
		}


    /***
     * initialize seed point(s)   
     ***/
    private void initialFMAux(int OJ, int OI) {

		// set seed point(s) to alive and its value(s) to 0 
		removeElement(farAway, new Point(OJ, OI));
		alive.addElement(new Point(OJ, OI));
		T[OI][OJ] = 0;

		// set seed's neighbours to close and compute its values
		removeElement(farAway, new Point(OJ - 1, OI));
		removeElement(farAway, new Point(OJ + 1, OI));
		removeElement(farAway, new Point(OJ, OI - 1));
		removeElement(farAway, new Point(OJ, OI + 1));
	
		close.addElement(new Point(OJ - 1, OI));
		close.addElement(new Point(OJ + 1, OI));
		close.addElement(new Point(OJ, OI - 1));
		close.addElement(new Point(OJ, OI + 1));
	
		T[OI][OJ - 1] = DELTA_X / Math.exp(-1 * ALPHA * KI[OI][OJ - 1]);
		T[OI][OJ + 1] = DELTA_X / Math.exp(-1 * ALPHA * KI[OI][OJ + 1]);
		T[OI - 1][OJ] = DELTA_Y / Math.exp(-1 * ALPHA * KI[OI - 1][OJ]);
		T[OI + 1][OJ] = DELTA_Y / Math.exp(-1 * ALPHA * KI[OI + 1][OJ]);
    }
    
    
    /***
     * check if a certain point is alive
     ***/
    private boolean isAlivePoint(Point p) {
		int numAlivePoints = alive.size();

		for (int i = 0; i < numAlivePoints; i++) {
			if (p.equals((Point) alive.elementAt(i))) return true;
		}

		return false;
    }


    /***
     * check if a certain point is close
     ***/
    private boolean isClosePoint(Point p) {
		int numClosePoints = close.size();

		for (int i = 0; i < numClosePoints; i++) {
			if (p.equals((Point) close.elementAt(i))) return true;
		}

		return false;
    }
    
    
    /***
     * fast marching implementation
     ***/
    private void march() {
	
		for (int i = 0; i < FMSteps; i++) {	   
 
			// select smallest T value point in close set
			int numClosePoints = close.size();
			double min = BIG_NUMBER;
			int rowMin = -1;
			int colMin = -1;
	    
			for (int j = 0; j < numClosePoints; j++) {
				Point tmpPoint = (Point) close.elementAt(j);
				if (min > T[tmpPoint.getY()][tmpPoint.getX()]) {
					min = T[tmpPoint.getY()][tmpPoint.getX()];
					rowMin = tmpPoint.getY();
					colMin = tmpPoint.getX();
				}
			}
	    
			// check validity of min T value	    
			if (rowMin == -1 && colMin == -1) {
				System.out.println("finding smallest T value pixel fails!");
				System.exit(1);
			}
	    
			// add smallest value point to alive, and remove it from close
			alive.addElement(new Point(colMin, rowMin));
			removeElement(close, new Point(colMin, rowMin));
	    
			// check its 4 neighbour points
			// first reset neighbour set
			neighbour.removeAllElements();

			checkNeighbour(new Point(colMin - 1, rowMin));
			checkNeighbour(new Point(colMin + 1, rowMin));
			checkNeighbour(new Point(colMin, rowMin - 1));
			checkNeighbour(new Point(colMin, rowMin + 1));
	    
			//recompute T value at neighbour points
			int numNeighbourPoints = neighbour.size();
	    
			for (int k = 0; k < numNeighbourPoints; k++) {
				Point tmpPoint = (Point) neighbour.elementAt(k);
				T[tmpPoint.getY()][tmpPoint.getX()] = updateT(tmpPoint.getY(), tmpPoint.getX()); 
			}
		}
    }
    
    
    /***
     * check neighbours of current alive point
     ***/
    private void checkNeighbour(Point p) {
		if (isAlivePoint(p)) ;
		else if (isClosePoint(p)) {
			neighbour.addElement(p);
		} else {
			removeElement(farAway, p);
			close.addElement(p);
			neighbour.addElement(p);
		}
    }
    
    
    /***
     * remove a point from a set
     ***/
    private void removeElement(Vector set, Point p) {
		int numElements = set.size();