levelset.java

LevelSet algorithme written in java language

		for (int i = 0; i < numElements; i++) {
			Point tmpPoint = (Point) set.elementAt(i);
			if(p.equals(tmpPoint)) { 
				set.removeElement(tmpPoint);

				// check removal successful
				if (set.size() != numElements - 1) {
					System.out.println("Remove point fails!");
					System.exit(1);
				}
		
				return;
			}
		}
	
		System.out.println("Remove point fails!");
		System.exit(1);
    }
    

    /***
     * update T values
     ***/    
    private double updateT (int i, int j) {
		double root [] = new double[16]; 
		boolean valid [] = new boolean[16];

		// check validity of row and column 
		if (i <= 0 || i >= (hauteur - 1) || j <= 0 || j >= (largeur - 1)) {
			System.out.println("update alter: image boundary reached!  i = " + i + ", j = " + j);
			System.exit(1);
		}

		// compute speed term
		double K = 1.0 / Math.exp(-1 * ALPHA * KI[i][j]);
	
		// case 1
		root[0] = f1(2, -2 * (T[i][j - 1] + T[i - 1][j]), T[i][j - 1] * T[i][j - 1] + T[i - 1][j] * T[i - 1][j] - K * K);
		root[1] = f2(2, -2 * (T[i][j - 1] + T[i - 1][j]), T[i][j - 1] * T[i][j - 1] + T[i - 1][j] * T[i - 1][j] - K * K);
	
		// case 2
		root[2] = f1(2, -2 * (T[i][j - 1] + T[i + 1][j]), T[i][j - 1] * T[i][j - 1] + T[i + 1][j] * T[i + 1][j] - K * K);
		root[3] = f2(2, -2 * (T[i][j - 1] + T[i + 1][j]), T[i][j - 1] * T[i][j - 1] + T[i + 1][j] * T[i + 1][j] - K * K);
	
		// case 3
		root[4] = f1(2, -2 * (T[i][j + 1] + T[i - 1][j]), T[i][j + 1] * T[i][j + 1] + T[i - 1][j] * T[i - 1][j] - K * K);
		root[5] = f2(2, -2 * (T[i][j + 1] + T[i - 1][j]), T[i][j + 1] * T[i][j + 1] + T[i - 1][j] * T[i - 1][j] - K * K);

		// case 4
		root[6] = f1(2, -2 * (T[i][j + 1] + T[i + 1][j]), T[i][j + 1] * T[i][j + 1] + T[i + 1][j] * T[i + 1][j] - K * K);
		root[7] = f2(2, -2 * (T[i][j + 1] + T[i + 1][j]), T[i][j + 1] * T[i][j + 1] + T[i + 1][j] * T[i + 1][j] - K * K);

		// case 5
		root[8] = T[i][j - 1] + K;
		root[9] = T[i][j - 1] - K;

		// case 6
		root[10] = T[i][j + 1] + K;
		root[11] = T[i][j + 1] - K;

		// case 7
		root[12] = T[i - 1][j] + K;
		root[13] = T[i - 1][j] - K;

		// casr 8
		root[14] = T[i + 1][j] + K;
		root[15] = T[i + 1][j] - K;

		// check validity 
		for (int m = 0; m < 16; m++) valid[m] = false;

		// case 1	
		if (T[i][j - 1] <= T[i][j + 1] && root[0] >= T[i][j - 1] && T[i - 1][j] <= T[i + 1][j] && root[0] >= T[i - 1][j]) valid[0] = true;
		if (T[i][j - 1] <= T[i][j + 1] && root[1] >= T[i][j - 1] && T[i - 1][j] <= T[i + 1][j] && root[1] >= T[i - 1][j]) valid[1] = true;
	
		// case 2
		if (T[i][j - 1] <= T[i][j + 1] && root[2] >= T[i][j - 1] && T[i - 1][j] >= T[i + 1][j] && root[2] >= T[i + 1][j]) valid[2] = true;
		if (T[i][j - 1] <= T[i][j + 1] && root[3] >= T[i][j - 1] && T[i - 1][j] >= T[i + 1][j] && root[3] >= T[i + 1][j]) valid[3] = true;

		// case 3
		if (T[i][j - 1] >= T[i][j + 1] && root[4] >= T[i][j + 1] && T[i - 1][j] <= T[i + 1][j] && root[4] >= T[i - 1][j]) valid[4] = true;
		if (T[i][j - 1] >= T[i][j + 1] && root[5] >= T[i][j + 1] && T[i - 1][j] <= T[i + 1][j] && root[5] >= T[i - 1][j]) valid[5] = true;

		// case 4
		if (T[i][j - 1] >= T[i][j + 1] && root[6] >= T[i][j + 1] && T[i - 1][j] >= T[i + 1][j] && root[6] >= T[i + 1][j]) valid[6] = true;
		if (T[i][j - 1] >= T[i][j + 1] && root[7] >= T[i][j + 1] && T[i - 1][j] >= T[i + 1][j] && root[7] >= T[i + 1][j]) valid[7] = true;

		// case 5
		if (T[i][j - 1] <= T[i][j + 1] && root[8] >= T[i][j - 1] && root[8] <= T[i - 1][j] && root[8] <= T[i + 1][j]) valid[8] = true;
		if (T[i][j - 1] <= T[i][j + 1] && root[9] >= T[i][j - 1] && root[9] <= T[i - 1][j] && root[9] <= T[i + 1][j]) valid[9] = true;

		// case 6	
		if (T[i][j - 1] >= T[i][j + 1] && root[10] >= T[i][j + 1] && root[10] <= T[i - 1][j] && root[10] <= T[i + 1][j]) valid[10] = true;
		if (T[i][j - 1] >= T[i][j + 1] && root[11] >= T[i][j + 1] && root[11] <= T[i - 1][j] && root[11] <= T[i + 1][j]) valid[11] = true;
    
		// case 7	
		if (root[12] <= T[i][j - 1] && root[12] <= T[i][j + 1] && T[i - 1][j] <= T[i + 1][j] && root[12] >= T[i - 1][j]) valid[12] = true;
		if (root[13] <= T[i][j - 1] && root[13] <= T[i][j + 1] && T[i - 1][j] <= T[i + 1][j] && root[13] >= T[i - 1][j]) valid[13] = true;

		// case 8
		if (root[14] <= T[i][j - 1] && root[14] <= T[i][j + 1] && T[i - 1][j] >= T[i + 1][j] && root[14] >= T[i + 1][j]) valid[14] = true;
		if (root[15] <= T[i][j - 1] && root[15] <= T[i][j + 1] && T[i - 1][j] >= T[i + 1][j] && root[15] >= T[i + 1][j]) valid[15] = true;
	
		// find the largest T value
		double largestRoot = SMALL_NUMBER;
		for (int n = 0; n < 16; n++) {
			if (valid[n]) {
				if (largestRoot < root[n]) largestRoot = root[n];
			}	
		}
	
		// check validity of new T value 
		if (largestRoot == SMALL_NUMBER) {
			System.out.println(" Fail to find a valid root");
			System.exit(1);
		}
	
		return largestRoot;
    }


    /***
     * compute first root of quadratic equation 
     ***/ 
    private double f1 (double a, double b, double c) {
		return (-1 * b + Math.sqrt(b * b - 4 * a * c)) / (2 * a);
    }


    /***
     * compute second root of quadratic equation 
     ***/ 
    private double f2 (double a, double b, double c) {
		return (-1 * b - Math.sqrt(b * b - 4 * a * c)) / (2 * a);
    }
    
    

    //*****************************************************************************
    // Level set step
    //*****************************************************************************
	
    /***
     * initialize Phi values 
     ***/
		private void initializePhi() {
			for (int j = 0; j < hauteur; j++) {
				for (int i = 0; i < largeur; i++) phi[i][j] = 0.0;
			}
		}
    
    
    /***
     * draw the initial contour, which is actually the close points set from fast marching step
     ***/
    private void initializeFront() {
		int numClosePoints = close.size();
		for (int i = 0; i < numClosePoints; i++) frontLine.addElement((Point) close.elementAt(i));
    }
    
    
    /***
     * initialize phi values for all non front points 
     ***/	
    private void initializeNonFront() {
		double min_distance = 0.0;
		int numFrontPoints = frontLine.size();
		int numAlivePoints = alive.size();
		int numFarAwayPoints = farAway.size();
	
		for (int i = 0; i < numAlivePoints; i++) {
	    
			min_distance = BIG_NUMBER;
			Point tmpPoint = (Point) alive.elementAt(i);
	    
			// traverse the front points to find minimum distance
			for (int k = 0; k < numFrontPoints; k++) { 
				double dist = tmpPoint.distance((Point) frontLine.elementAt(k));
				if (dist < min_distance) min_distance = dist;
			}
	    
			// set image_phi value at this image point
			phi[tmpPoint.getY()][tmpPoint.getX()] = -1 * min_distance;  // inner pixel point
		}
	
		for (int i = 0; i < numFarAwayPoints; i++) {
	    
			min_distance = BIG_NUMBER;
			Point tmpPoint = (Point) farAway.elementAt(i);
	    
			// traverse the front points to find minimum distance
			for (int k = 0; k < numFrontPoints; k++) { 
				double dist = tmpPoint.distance((Point) frontLine.elementAt(k));
				if (dist < min_distance) min_distance = dist;
			}
	    
			// set image_phi value at this image point
			phi[tmpPoint.getY()][tmpPoint.getX()] = min_distance;  // outer pixel point
		}
    }
    

    /***
     * check if a certain point is in given front line 
     ***/
    private boolean isFrontPoint(Point p) {
		int numFrontPoints = frontLine.size();

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

		return false;
    }

    
    /***
     * check if point (i, j) is front point 
     ***/
    private boolean newFrontPoint(int i, int j) {
		// check image boundary
		if (i == hauteur - 1 || j == largeur - 1) return false; 
	
		double max_phi = Math.max(Math.max(phi[i][j], phi[i + 1][j]), Math.max(phi[i][j + 1], phi[i + 1][j + 1]));
		double min_phi = Math.min(Math.min(phi[i][j], phi[i + 1][j]), Math.min(phi[i][j + 1], phi[i + 1][j + 1]));
	
		if (max_phi > 0.0 && min_phi < 0.0) return true;
	
		return false;
    }
    
    
    /***
     * check all points, add front points to front line 
     ***/
    private void setFront() {
	
		// reset new front line
		frontLine.removeAllElements();
	
		for (int i = 0; i < hauteur; i++) {
			for (int j = 0; j < largeur; j++) {
				if (newFrontPoint(i, j)) frontLine.addElement(new Point(j, i));  // j = x (col), i = y (row) 
			}
		}
	
		// reset front point values 
		int numFrontPoints = frontLine.size();
		for (int i = 0; i < numFrontPoints; i++) {
			Point tmpPoint = (Point) frontLine.elementAt(i);
			phi[tmpPoint.getY()][tmpPoint.getX()] = 0;
		}
    }
    
    
    /***
     * compute phi values for all non front points 
     ***/	
    private void setNonFront() {
		double min_distance = 0.0;
		int numFrontPoints = frontLine.size();
	
		for (int i = 0; i < hauteur; i++) { // traverse the whole image
			for (int j = 0; j < largeur; j++) {
		
				Point tmpPoint = new Point(j, i);  // j = x (col), i = y (row)
		
				if (!isFrontPoint(tmpPoint)) {
					min_distance = BIG_NUMBER;
		    
					// traverse the front points to find minimum distance
					for (int k = 0; k < numFrontPoints; k++) { 
						double dist = tmpPoint.distance((Point) frontLine.elementAt(k));
						if (dist < min_distance) min_distance = dist;
					}
		    
					// set image_phi value at this image point
					if (phiOld[i][j] < 0) phi[i][j] = -1 * min_distance;  // inner pixel point
					else phi[i][j] = min_distance;  // outer pixel point
				}
			}  	
		}
    }
    
    
    /***
     * check phi values and reset narrow band 
     ***/
    private void setNarrowBand() {
		// reset narrow band
		narrowBand.removeAllElements();
		nbBoundary.removeAllElements();
		testNB.removeAllElements();	
	
		for (int j = 0; j < hauteur; j++) {
			for (int i = 0; i < largeur; i++) {
				if (Math.abs(phi[i][j]) <= bandWidth) narrowBand.addElement(new Point(j, i));
				if (Math.abs(phi[i][j]) >= (bandWidth - 2) && Math.abs(phi[i][j]) <= (bandWidth - 1)) {
					nbBoundary.addElement(new Point(j, i));
					testNB.addElement(new Double(phi[i][j]));
				}
			}
		}
    }
    
    
    /***
     * check if narrow band boundary reached 
     ***/
    private boolean reachNBBoundary() {
		int numNBBPoints = nbBoundary.size();
	
		for (int i = 0; i < numNBBPoints; i++) {
			Point tmpPoint = (Point) nbBoundary.elementAt(i);
			double oldValue = ((Double) testNB.elementAt(i)).doubleValue();
	    
			if (phi[tmpPoint.getY()][tmpPoint.getX()] * oldValue < 0) return true;
		}

		return false;
    }
    
    
    /***
     * update phi values at points within narrow band
     ***/    
    private void updatePhi() {
		int numNBPoints = narrowBand.size(); 
	
		for (int n = 0; n < numNBPoints; n++) {

			int i = ((Point) narrowBand.elementAt(n)).getY();  // row number 
			int j = ((Point) narrowBand.elementAt(n)).getX();  // column number
	    
			if (i > 2  && i < (hauteur - 3) && j > 2 && j < (largeur - 3)) {
		
				// first term (advection)
				double Dx_minus = (double) (phi[i][j] - phi[i][j - 1]) / DELTA_X;
				double Dx_plus  = (double) (phi[i][j + 1] - phi[i][j]) / DELTA_X;
				double Dy_minus = (double) (phi[i][j] - phi[i - 1][j]) / DELTA_Y;
				double Dy_plus  = (double) (phi[i + 1][j] - phi[i][j]) / DELTA_Y;
		
				double Dx_max = Math.max(Dx_minus, 0);
				double Dx_min = Math.min(Dx_plus, 0);
				double Dy_max = Math.max(Dy_minus, 0);
				double Dy_min = Math.min(Dy_plus, 0);
		
				double sum_adv = Math.sqrt(Dx_max * Dx_max + Dx_min * Dx_min + Dy_max * Dy_max + Dy_min * Dy_min);	    
		
				// second term (curvature)
				double delta_phi = Math.sqrt((phi[i + 1][j] - phi[i - 1][j]) * (phi[i + 1][j] - phi[i - 1][j]) / (4.0 * DELTA_X * DELTA_X)
											 + (phi[i][j + 1] - phi[i][j - 1]) * (phi[i][j + 1] - phi[i][j - 1]) / (4.0 * DELTA_Y * DELTA_Y));
		
				double phi_x = (phi[i][j + 1] - phi[i][j - 1]) / (2.0 * DELTA_X);
				double phi_y = (phi[i + 1][j] - phi[i - 1][j]) / (2.0 * DELTA_Y);
				double phi_xx = (phi[i][j + 1] + phi[i][j - 1] - 2 * phi[i][j]) / (1.0 * DELTA_X);
				double phi_yy = (phi[i + 1][j] + phi[i - 1][j] - 2 * phi[i][j]) / (1.0 * DELTA_Y);
				double phi_xy = (phi[i + 1][j + 1] + phi[i - 1][j - 1] - phi[i - 1][j + 1] - phi[i + 1][j - 1]) / (4.0 * DELTA_X * DELTA_Y);
		
				double k = 0.0;
				if (phi_x == 0 && phi_y == 0) k = 0;
				else k = (phi_xx * phi_y * phi_y - 2 * phi_x * phi_y * phi_xy + phi_yy * phi_x * phi_x) 
						 / Math.pow((phi_x * phi_x + phi_y * phi_y), 3.0/2.0);
		
				double K_I = 1.0 / (1.0 + Math.abs(KI[i][j]));
		
				// third term (attraction)
				double Ix_KI = -1.0 * (KI[i][j] - KI[i][j - 1]);
				double Iy_KI = -1.0 * (KI[i][j] - KI[i - 1][j]);
				double Ix_phi = phi[i][j] - phi[i][j - 1];
				double Iy_phi = phi[i][j] - phi[i - 1][j];	    
		
				double sum_att = Ix_KI * Ix_phi + Iy_KI * Iy_phi;	    
		
				// update phi
				phi[i][j] = phi[i][j] - DELTA_T * K_I * (sum_adv * FA - EPSILON * k * delta_phi) + DELTA_T * BETA * sum_att;
			}
		}
    }
    
    
    /***
     * store Phi values, i.e. save phi values for reinitialization 
     * after each outer level set iteration 
     ***/
    private void storePhi() {
		for (int i = 0; i < hauteur; i++) {
			for (int j = 0; j < largeur; j++) {
				phiOld[i][j] = phi[i][j];
			}
		}
    }


    /***
     * draw front lines on the image file for visualization
     ***/
    private void redraw(String filename1, String filename2) throws FileNotFoundException, IOException {
	
		// for reading from input image file	
		BufferedReader br = new BufferedReader(new FileReader(new File(filename1))); 

		// for writting to output image file
		PrintWriter pw = new PrintWriter(new BufferedWriter(new FileWriter(new File(filename2))));
	
		// read first three lines from input file and write to output file
		for (int i = 0; i < 3; i++) pw.println(br.readLine());
	
		// read input image and write proper values to output image 
		for (int j = 0 ; j < hauteur; j++) {
			String line = br.readLine();
			StringTokenizer tokens = new StringTokenizer(line);
	    
			for (int i = 0; i < largeur; i++) {
				int value = Integer.parseInt(tokens.nextToken());
	
				if (isFrontPoint(new Point(i, j))) pw.print(0 + " ");
				else pw.print(value + " ");
			}
	    
			pw.println();	    
		}
	
		pw.close();
    }
    

    public static void main(String [] args) throws FileNotFoundException, IOException {
	/*	if (args.length != 2) {
			System.err.println("Usage: LevelSet <input image file> <output image file>");
			System.exit(1);
		}*/
	
		LevelSet ls = new LevelSet();
	
		// read KI values
		ls.readFile("nomFichier");
	
		// first use fast marching for rough contour
		System.out.println("fast marching: begins!!!");
		ls.initializeFastMarch();
		ls.march();
		System.out.println("fast marching: finishes!!!");

		//restore KI value
		ls.restoreKi();

		// seconfly, use level set for fine tuning
		// initialize 
		System.out.println("level set: begins!!!");
		ls.initializePhi();	  	
		ls.initializeFront();
		ls.initializeNonFront();
		ls.setNarrowBand();
		System.out.println("level set: initialize finishes!");
	
		// iterate 	
		for (int iter = 0; iter < iter_outer; iter++) {
			int m = 0;
	    
			while (m < iter_inner) { 
				ls.updatePhi();
				m++;
			} // while

			// save phi value for re-initialization
			ls.storePhi();
	    
			ls.setFront();
			ls.setNonFront();
			ls.setNarrowBand();
			System.out.println("Iteration #" + iter + ": finished! m = " + m);
		} // for
	
		// if original image is not gray scaled, its blurred version should be used for display
		ls.redraw("nomFichier", "h2");
		System.out.println("level set: finishes!");
    }
    
}