content: package edu.stanford.rsl.conrad.utils;
  
  import ij.process.ImageProcessor;
  
  import java.util.ArrayList;
  
  import edu.stanford.rsl.conrad.geometry.General;
  import edu.stanford.rsl.conrad.geometry.Rotations;
  import edu.stanford.rsl.conrad.geometry.shapes.simple.PointND;
  import edu.stanford.rsl.conrad.geometry.transforms.ScaleRotate;
  import edu.stanford.rsl.conrad.geometry.transforms.Transform;
  
  public abstract class DoublePrecisionPointUtil {
  
  
  	/*
  	 * PCA:
  	 *  cov = a * a'
  	 *  eiv = eigen(cov)
  	 *  sort eiv system by eigenvalues, largest first.
  	 *  diag = diagonal matrix of eigenvalues
  	 *  psi = a' * eigenvectors * diag
  	 *  eigenimages = psi'
  	 */
  
  
  
  	/**
  	 * Compute the geometric center of a set of points
  	 * @param list the set of points
  	 * @return the geometric center
  	 */
  	public static PointND getGeometricCenter(ArrayList<PointND> list){
  		int dim = list.get(0).getDimension();
  		double [] temp = new double [list.get(0).getDimension()];
  		for (int i = 0; i < list.size(); i++){
  			for (int j = 0; j < dim; j++){
  				temp[j] += list.get(i).get(j);
  			}
  		}
  		for (int j = 0; j < dim; j++){
  			temp[j] /= list.size();
  		}
  		return new PointND(temp);
  	}
  
  	/**
  	 * Compute the standard deviation of a set of points
  	 * @param list the set of points
  	 * @return the geometric center
  	 */
  	public static PointND getStandardDeviation(ArrayList<PointND> list){
  		int dim = list.get(0).getDimension();
  		PointND center = getGeometricCenter(list);
  		double [] temp = new double [list.get(0).getDimension()];
  		for (int i = 0; i < list.size(); i++){
  			for (int j = 0; j < dim; j++){
  				temp[j] += Math.pow(list.get(i).get(j)-center.get(j),2);
  			}
  		}
  		for (int j = 0; j < dim; j++){
  			temp[j] = Math.sqrt(temp[j]) / list.size();
  		}
  		return new PointND(temp);
  	}
  
  
  	/**
  	 * Extract points from an ImageProcessor which exceed a certain value
  	 * 
  	 * @param houghSpace the ImageProcessor
  	 * @param offset the threshold for extraction
  	 * @return the list of candidate points
  	 */
  	public static ArrayList<PointND> extractCandidatePoints(ImageProcessor houghSpace, double offset){
  		ArrayList<PointND> candidate = new ArrayList<PointND>();
  		for (int j = 0; j< houghSpace.getHeight(); j++){
  			for (int i = 0; i< houghSpace.getWidth(); i++){
  				if (houghSpace.getPixelValue(i, j) > offset) {
  					PointND point = new PointND(i, j);
  					candidate.add(point);					
  				}
  			}
  		}
  		return candidate;
  	}
  
  	/**
  	 * Extracts cluster centers from an ordered List of points. Points must be ordered first with respect to x, then to y coordinate. Algorithm assumes that only one point may appear in the same row, i.e.,  all clusters must be separable via the y direction.
  	 * A cluster center is then computed as the geometric center of the points in the same cluster. Algorithm is fast, but very restricted.
  	 * @param pointList the list of candidate points
  	 * @param distance the minimal distance between clusters
  	 * @return the list of cluster centers
  	 */
  	public static ArrayList<PointND> extractClusterCenter(ArrayList<PointND> pointList, double distance){
  		ArrayList<PointND> centerPoint = new ArrayList<PointND>();
  		while (pointList.size() > 0){
  			PointND reference = pointList.get(0);
  			ArrayList<PointND> currentSubset = new ArrayList<PointND>();
  			//currentSubset.add(reference);
  			for (int i = 0; i < pointList.size(); i++){
  				PointND current = pointList.get(i);
  				if (current.euclideanDistance(reference) < distance){
  					currentSubset.add(current);
  					pointList.remove(i);
  					i--;
  				} else {
  					// points are ordered first in x and then in y direction
  					// hence, end of current cluster if more than distance away in y direction
  					if (Math.abs(reference.get(1) - current.get(1)) > distance) break;
  				}
  			}
  			centerPoint.add(getGeometricCenter(currentSubset));
  		}
  		return centerPoint;
  	}
  
  
  	/**
  	 * Computes the total distance between two list of points.
  	 * We assume that the respectively same entry of each list refers to the point to be compared in the other list.
  	 * @param list1 the one list
  	 * @param list2 the other list
  	 * @return the sum of all euclidian distances.
  	 */
  	public static double computePointWiseDifference(ArrayList<PointND> list1, ArrayList<PointND> list2){
  		double revan =0;
  		for (int i=0;i<list1.size();i++){
  			revan += list1.get(i).euclideanDistance(list2.get(i));
  		}
  		return revan;
  	}
  
  	/**
  	 * Transforms a list of given points and returns them as new instances in a new list of points.
  	 * @param list the list
  	 * @param t the transform
  	 * @return the new list of points
  	 */
  	public static ArrayList<PointND> transformPoints(ArrayList<PointND> list, Transform t){
  		ArrayList<PointND> revan = new ArrayList<PointND>();
  		for (int i=0; i<list.size();i++){
  			PointND newP = t.transform(list.get(i));
  			revan.add(newP);
  		}
  		return revan;
  	}
  
  	/**
  	 * Estimates an optimal rotation transform to transform list1 onto list2.
  	 * @param list1 the first list of points
  	 * @param list2 the second list of points
  	 * @param maxAngle the maximal angle that is searched (in radians)
  	 * @param iterations the maximal number of iterations per step.
  	 * @return the scale rotation
  	 */
  	public static ScaleRotate estimateRotation(ArrayList<PointND> list1, ArrayList<PointND> list2, double maxAngle, int iterations){
  		double angleX = 0;
  		double angleY = 0;
  		double angleZ = 0;
  		ScaleRotate transform = new ScaleRotate(Rotations.createRotationMatrix(angleX, angleY, angleZ));
  		for (int i = 0; i < iterations; i++){
  			double angleMinX = -maxAngle;
  			double angleMaxX = maxAngle;
  			angleX = 0;
  			double errorLeft = DoublePrecisionPointUtil.computePointWiseDifference(
  					DoublePrecisionPointUtil.transformPoints(
  							list1, new ScaleRotate(
  									Rotations.createRotationMatrix(angleMinX, angleY, angleZ))), 
  									list2);
  			double errorCenter = DoublePrecisionPointUtil.computePointWiseDifference(
  					DoublePrecisionPointUtil.transformPoints(
  							list1, new ScaleRotate(
  									Rotations.createRotationMatrix(angleX, angleY, angleZ))), 
  									list2);
  			double errorRight = DoublePrecisionPointUtil.computePointWiseDifference(
  					DoublePrecisionPointUtil.transformPoints(
  							list1, new ScaleRotate(
  									Rotations.createRotationMatrix(angleMaxX, angleY, angleZ))), 
  									list2);
  			for (int j = 0; j <iterations; j++){
  				if (errorLeft < errorRight){
  					errorRight = errorCenter;
  					angleMaxX = angleX;
  					angleX = (angleX+angleMinX) / 2.0;
  				} else {
  					errorLeft = errorCenter;
  					angleMinX = angleX;
  					angleX = (angleX+angleMaxX) / 2.0;
  
  				}	
  				errorCenter = DoublePrecisionPointUtil.computePointWiseDifference(
  						DoublePrecisionPointUtil.transformPoints(
  								list1, new ScaleRotate(
  										Rotations.createRotationMatrix(angleX, angleY, angleZ))), 
  										list2);
  			}
  			double angleMinY = -maxAngle;
  			double angleMaxY = maxAngle;
  			angleY = 0;
  			errorLeft = DoublePrecisionPointUtil.computePointWiseDifference(
  					DoublePrecisionPointUtil.transformPoints(
  							list1, new ScaleRotate(
  									Rotations.createRotationMatrix(angleX, angleMinY, angleZ))), 
  									list2);
  			errorCenter = DoublePrecisionPointUtil.computePointWiseDifference(
  					DoublePrecisionPointUtil.transformPoints(
  							list1, new ScaleRotate(
  									Rotations.createRotationMatrix(angleX, angleY, angleZ))), 
  									list2);
  			errorRight = DoublePrecisionPointUtil.computePointWiseDifference(
  					DoublePrecisionPointUtil.transformPoints(
  							list1, new ScaleRotate(
  									Rotations.createRotationMatrix(angleX, angleMaxY, angleZ))), 
  									list2);
  			for (int j = 0; j <iterations; j++){
  				if (errorLeft < errorRight){
  					errorRight = errorCenter;
  					angleMaxY = angleY;
  					angleY = (angleY+angleMinY) / 2.0;
  				} else {
  					errorLeft = errorCenter;
  					angleMinY = angleY;
  					angleY = (angleY+angleMaxY) / 2.0;
  				}
  				errorCenter = DoublePrecisionPointUtil.computePointWiseDifference(
  						DoublePrecisionPointUtil.transformPoints(
  								list1, new ScaleRotate(
  										Rotations.createRotationMatrix(angleX, angleY, angleZ))), 
  										list2);
  			}
  			double angleMinZ = -maxAngle;
  			double angleMaxZ = maxAngle;
  			angleZ = 0;
  			errorLeft = DoublePrecisionPointUtil.computePointWiseDifference(
  					DoublePrecisionPointUtil.transformPoints(
  							list1, new ScaleRotate(
  									Rotations.createRotationMatrix(angleX, angleY, angleMinZ))), 
  									list2);
  			errorCenter = DoublePrecisionPointUtil.computePointWiseDifference(
  					DoublePrecisionPointUtil.transformPoints(
  							list1, new ScaleRotate(
  									Rotations.createRotationMatrix(angleX, angleY, angleZ))), 
  									list2);
  			errorRight = DoublePrecisionPointUtil.computePointWiseDifference(
  					DoublePrecisionPointUtil.transformPoints(
  							list1, new ScaleRotate(
  									Rotations.createRotationMatrix(angleX, angleY, angleMaxZ))), 
  									list2);
  			for (int j = 0; j <iterations; j++){
  				if (errorLeft < errorRight){
  					errorRight = errorCenter;
  					angleMaxZ = angleZ;
  					angleZ = (angleZ+angleMinZ) / 2.0;
  				} else {
  					errorLeft = errorCenter;
  					angleMinZ = angleZ;
  					angleZ = (angleZ+angleMaxZ) / 2.0;
  				}
  				errorCenter = DoublePrecisionPointUtil.computePointWiseDifference(
  						DoublePrecisionPointUtil.transformPoints(
  								list1, new ScaleRotate(
  										Rotations.createRotationMatrix(angleX, angleY, angleZ))), 
  										list2);
  			}	
  		}
  		transform = new ScaleRotate(Rotations.createRotationMatrix(angleX, angleY, angleZ));
  		if (CONRAD.DEBUGLEVEL > 0) {
  			double finalError = DoublePrecisionPointUtil.computePointWiseDifference(
  					DoublePrecisionPointUtil.transformPoints(
  							list1, transform), list2);
  			System.out.println(General.toDegrees(angleX) + " " + General.toDegrees(angleY) + " " + General.toDegrees(angleZ) + " " + finalError / list1.size());
  		}
  		return transform;
  	}
  
  
  }
  /*
   * Copyright (C) 2010-2014  Andreas Maier
   * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
  */