content: package edu.stanford.rsl.conrad.geometry;
  
  
  import ij.process.BinaryProcessor;
  import ij.process.ByteProcessor;
  import ij.process.ImageProcessor;
  
  import java.io.BufferedWriter;
  import java.io.IOException;
  import java.util.ArrayList;
  import java.util.Iterator;
  
  import edu.stanford.rsl.conrad.geometry.shapes.simple.PointND;
  import edu.stanford.rsl.conrad.geometry.shapes.simple.StraightLine;
  import edu.stanford.rsl.conrad.geometry.shapes.simple.Triangle;
  import edu.stanford.rsl.conrad.numerics.SimpleMatrix;
  import edu.stanford.rsl.conrad.numerics.SimpleOperators;
  import edu.stanford.rsl.conrad.numerics.SimpleVector;
  import edu.stanford.rsl.conrad.utils.CONRAD;
  
  
  /**
   * All kinds of geometric routines that are not specific to some geometric object or interact on a set of these.
   * @author Andreas Keil
   */
  public abstract class General {
  
  	public static final SimpleVector E_X = new SimpleVector(1.0, 0.0, 0.0);
  	public static final SimpleVector E_Y = new SimpleVector(0.0, 1.0, 0.0);
  	public static final SimpleVector E_Z = new SimpleVector(0.0, 0.0, 1.0);
  
  	private static boolean normalizeMode = false;
  
  	public static boolean isNormalizeMode(){
  		return normalizeMode;
  	}
  
  	public static void setNormalizeMode(boolean mode){
  		normalizeMode = mode;
  	}
  	
  	/**
  	 * Computes the convex hull of the projection of the eight corners of the volume to the projection.
  	 * Here we assume that the volume is aligned with the world coordinate system and can be described
  	 * by a maximal and a minimal coordinate.
  	 * 
  	 * Note that the returned int array contains the 4 values of the convex hull in projection domain
  	 * as sequence minimum x, maximum x, minimum y, maximum y 
  	 * 
  	 * @param min the minimal 3D coordinate of the volume
  	 * @param max the maximal 3D coordinate of the volume
  	 * @param proj the projection matrix
  	 * @return the convex hull in projection domain.
  	 */
  	public static int[] projectVolumeToProjection(PointND min, PointND max, Projection proj){
  		// bounds in 2D:
  		int [] bounds = new int [4];
  		// eight corners of the volume in homogeneous coordinates:
  		PointND [] corners = new PointND[8];
  		corners[0] = new PointND(min.get(0),min.get(1),min.get(2), 1);
  		corners[1] = new PointND(min.get(0),min.get(1),max.get(2), 1);
  		corners[2] = new PointND(min.get(0),max.get(1),min.get(2), 1);
  		corners[3] = new PointND(min.get(0),max.get(1),max.get(2), 1);
  		corners[4] = new PointND(max.get(0),min.get(1),min.get(2), 1);
  		corners[5] = new PointND(max.get(0),min.get(1),max.get(2), 1);
  		corners[6] = new PointND(max.get(0),max.get(1),min.get(2), 1);
  		corners[7] = new PointND(max.get(0),max.get(1),max.get(2), 1);
  		SimpleVector projected = new SimpleVector(2);
  		@SuppressWarnings("unused")
  		double depth = proj.project(corners[0].getAbstractVector(), projected);
  // TODO: Check depth value to see if point is actually in front of the camera (>0) or behind it (<0)
  //		if (depth <= 0.0) ; // add code for invisible point here
  // old code:
  //		SimpleMatrix p = proj.computeP();
  //		SimpleVector projected = SimpleOperators.multiply(p, corners[0].getAbstractVector());
  //		projected.divideBy(projected.getElement(2));
  		bounds[0] = (int) projected.getElement(0);
  		bounds[1] = (int) projected.getElement(0);
  		bounds[2] = (int) projected.getElement(1);
  		bounds[3] = (int) projected.getElement(1);
  		for (int i = 1; i < 8; i++){
  			depth = proj.project(corners[i].getAbstractVector(), projected);
  // TODO: Check depth value to see if point is actually in front of the camera (>0) or behind it (<0)
  //			if (depth <= 0.0) ; // add code for invisible point here
  // old code:
  //			projected = SimpleOperators.multiply(p, corners[i].getAbstractVector());
  //			projected.divideBy(projected.getElement(2));
  // proposed new code for min/max operations:
  			bounds[0] = java.lang.Math.min(bounds[0], (int)projected.getElement(0)); // TODO: Maybe better use ceil/floor instead of a simple cast here?
  			bounds[1] = java.lang.Math.max(bounds[1], (int)projected.getElement(0)); // TODO: Maybe better use ceil/floor instead of a simple cast here?
  			bounds[2] = java.lang.Math.min(bounds[2], (int)projected.getElement(1)); // TODO: Maybe better use ceil/floor instead of a simple cast here?
  			bounds[3] = java.lang.Math.max(bounds[3], (int)projected.getElement(1)); // TODO: Maybe better use ceil/floor instead of a simple cast here?
  // old code:
  //			bounds[0] = (bounds[0] > projected.getElement(0)) ? (int) projected.getElement(0) : bounds[0];
  //			bounds[1] = (bounds[1] < projected.getElement(0)) ? (int) projected.getElement(0) : bounds[1];
  //			bounds[2] = (bounds[2] > projected.getElement(1)) ? (int) projected.getElement(1) : bounds[2];
  //			bounds[3] = (bounds[3] < projected.getElement(1)) ? (int) projected.getElement(1) : bounds[3];
  		}
  		return bounds;
  	}
  	
  
  	public static SimpleVector crossProduct(final SimpleVector v1, final SimpleVector v2) {
  		// input check
  		assert (v1.getLen() == 3) : new IllegalArgumentException("v1 has to be a 3-vector!");
  		assert (v2.getLen() == 3) : new IllegalArgumentException("v2 has to be a 3-vector!");
  
  		final SimpleVector result = new SimpleVector(3);
  		result.setElementValue(0, v1.getElement(1) * v2.getElement(2) - v1.getElement(2) * v2.getElement(1));
  		result.setElementValue(1, v1.getElement(2) * v2.getElement(0) - v1.getElement(0) * v2.getElement(2));
  		result.setElementValue(2, v1.getElement(0) * v2.getElement(1) - v1.getElement(1) * v2.getElement(0));
  		return result;
  	}
  
  	public static boolean areColinear(final SimpleVector v1, final SimpleVector v2, final double delta) {
  		return (crossProduct(v1, v2).normL2() < delta);
  	}
  
  	/**
  	 * Computes the angle between two vectors;
  	 * @param a
  	 * @param b
  	 * @return the smaller (and positive) angle between a and b in radians
  	 */
  	public static double angle(SimpleVector a, SimpleVector b) {
  		final double norma = a.normL2();
  		final double normb = b.normL2();
  		if (norma == 0 || normb == 0)
  			throw new RuntimeException("At least one of the vectors has zero length!");
  		else
  			return (Math.acos(SimpleOperators.multiplyInnerProd(a, b)/(norma * normb)));
  	}
  	
  	public static double euclideanDistance(final SimpleVector v1, final SimpleVector v2) {
  		return SimpleOperators.subtract(v1, v2).normL2();
  	}
  
  	public static SimpleVector augmentToHomgeneous(final SimpleVector v) {
  		SimpleVector v_hom = new SimpleVector(v.getLen() + 1);
  		v_hom.setSubVecValue(0, v);
  		v_hom.setElementValue(v.getLen(), 1.0);
  		return v_hom;
  	}
  
  	public static SimpleVector normalizeFromHomogeneous(final SimpleVector v) {
  		SimpleVector v_normalized = v.getSubVec(0, v.getLen()-1);
  		assert (v.getElement(v.getLen()-1) > CONRAD.DOUBLE_EPSILON) : new RuntimeException("Cannot de-homogenize a point at infinity!");
  		v_normalized.divideBy(v.getElement(v.getLen()-1));
  		return v_normalized;
  	}
  
  	public static SimpleMatrix createHomAffineMotionMatrix(final SimpleMatrix A, final SimpleVector t) {
  		// input checks
  		final int d = t.getLen(); 
  		assert (d == 2 || d == 3);
  		assert (A.getRows() == d && A.getCols() == d);
  
  		final SimpleMatrix result = new SimpleMatrix(d+1, d+1);
  		result.setSubMatrixValue(0, 0, A);
  		result.setSubColValue(0, d, t);
  		result.setElementValue(d, d, 1.0);
  		return result;
  	}
  
  	public static SimpleMatrix createHomAffineMotionMatrix(final SimpleMatrix A) {
  		// input checks
  		final int d = A.getRows(); 
  		assert (d == 2 || d == 3);
  		assert A.getCols() == d;
  
  		final SimpleMatrix result = new SimpleMatrix(d+1, d+1);
  		result.setSubMatrixValue(0, 0, A);
  		result.setElementValue(d, d, 1.0);
  		return result;
  	}
  
  	public static SimpleMatrix createHomAffineMotionMatrix(final SimpleVector t) {
  		// input checks
  		final int d = t.getLen(); 
  		assert (d == 2 || d == 3);
  
  		final SimpleMatrix result = new SimpleMatrix(d+1, d+1);
  		result.identity();
  		result.setSubColValue(0, d, t);
  		return result;
  	}
  
  	public static double toRadians(double degrees){
  		return (degrees / 180.0) * Math.PI;
  	}
  	
  	public static double toDegrees(double radians){
  		return (radians * 180.0) / Math.PI;
  	}
  		
  	/**
  	 * Convert voxel indexes to world coordinates (in mm), given the spacing and
  	 * origin of a volume.
  	 * 
  	 * The origin in world coordinates is defined as the world coordinate of the
  	 * center of the first voxel (the 0/0/0 voxel) in world coordinates (in mm).
  	 * Example: To center a volume with 512 voxels and 1.5mm spacing in every
  	 * direction around the 0mm/0mm/0mm world coordinate, one has to set all the
  	 * origin parameters to -383.25mm (= -512/2*1.5mm + 1/2*1.5mm).
  	 * 
  	 * @see #worldToVoxel(double[] world, double[] spacing, double[] origin)
  	 */
  	public static double[] voxelToWorld(int[] voxel, double[] spacing, double[] origin) {
  		double[] world = new double[3];
  		for (int i = 0; i < 3; ++i)
  			world[i] = voxelToWorld(voxel[i],spacing[i], origin[i]);
  		return world;
  	}
  
  
  	/**
  	 * Computes the origin in world coordinates from pixel coordinates. Example<BR>
  	 * General.voxelToWorld(-originInPixelsX, getVoxelSpacingX(), 0)
  	 * @param voxel the negative voxel coordinate
  	 * @param spacing the spacing
  	 * @param origin 0
  	 * @return the origin in world coordinates
  	 */
  	public static double voxelToWorld(double voxel, double spacing, double origin) {
  		return voxel*spacing + origin;
  	}
  	
  	
  	/**
  	 * Convert world coordinates (in mm) to voxel indexes, given the spacing and
  	 * origin of a volume.
  	 * 
  	 * @see #voxelToWorld(int[] voxel, double[] spacing, double[] origin) 
  	 */
  	public static double[] worldToVoxel(double[] world, double[] spacing, double[] origin) {
  		double[] voxel = new double[3];
  		for (int i = 0; i < 3; ++i)
  			voxel[i] = worldToVoxel(world[i], spacing[i], origin[i]);
  		return voxel;
  	}
  	
  	/**
  	 * Helper function to convert world coordinates (in mm) to voxel indexes.
  	 * Example: <BR>
  	 * General.worldToVoxel(-originInWorld.get(0), getVoxelSpacingX(), 0)
  	 * 
  	 * @see #worldToVoxel(double[] world, double[] spacing, double[] origin)
  	 */
  	public static double worldToVoxel(double world, double spacing, double origin) {
  		return (world - origin) / spacing;
  	}
  	
  	public static void splitHomAffineMotionMatrix(final SimpleMatrix At, final SimpleMatrix A, final SimpleVector t) {
  		// check input
  		assert (At.getRows() == 4 && At.getCols() == 4 && At.getElement(3, 0) == 0 && At.getElement(3,1) == 0 && At.getElement(3,2) == 0 && At.getElement(3,3) > 10.0*CONRAD.DOUBLE_EPSILON) :
  			new IllegalArgumentException("At has to be a homogeneous rigid motion matrix!");
  
  		// check output
  		assert (A.getRows() == 3 && A.getCols() == 3) :
  			new IllegalArgumentException("A has to be a 3x3 matrix!");
  		assert (t.getLen() == 3) :
  			new IllegalArgumentException("t has to be a 3-vector!");
  
  		// extract normalized A and t
  		final double scale = 1.0/At.getElement(3, 3);
  		A.init(At.getSubMatrix(0, 0, 3, 3).multipliedBy(scale));
  		t.init(At.getSubCol(0, 3, 3).multipliedBy(scale));
  	}
  
  	public static ArrayList<PointND> intersectRayWithCuboid(StraightLine line, PointND min, PointND max){
  		double [] vals = new double [2];
  		ArrayList<PointND> revan = new ArrayList<PointND>();
  		if (General.intersectRayWithCuboid(line.getPoint().getAbstractVector(), line.getDirection(), min.getAbstractVector(), max.getAbstractVector(), vals)){
  			revan.add(line.evaluate(vals[0]));
  			revan.add(line.evaluate(vals[1]));
  		}
  		return revan;
  	}
  	
  	/**
  	 * Method to check whether a point is within a given cubiod defined by min and max. 
  	 * @param point the point
  	 * @param min the minimum coordinate of the cuboid
  	 * @param max the maximum coordinate of the cuboid.
  	 * @return true, if the point is inside
  	 */
  	public static boolean isWithinCuboid(PointND point, PointND min, PointND max){
  		boolean fulfilled = true;
  		for (int i =0; i < point.getDimension(); i++){
  			fulfilled = fulfilled && (point.get(i)<max.get(i)) && (point.get(i)>min.get(i));
  		}
  		return fulfilled;
  	}
  	
  	/**
  	 * Computes the two intersections of a ray with a cuboid, called entry and
  	 * exit point where the ray is specified by the given line origin and ray direction.
  	 *
  	 * The entry and exit points are returned as distances from the line origin along the
  	 * ray. The world coordinates of the entry and exit points may be computed as follows:
  	 * {@latex.ilb %preamble{\\usepackage{amsmath}} \\begin{align*}
  	 *   \\mathbf{v}_n & = \\mathbf{C} + t_n \\cdot \\mathbf{d} \\
  	 * 	 \\mathbf{v}_f & = \\mathbf{C} + t_f \\cdot \\mathbf{d}
  	 * \\end{align*} }
  	 * 
  	 * @param origin  The ray origin in world coordinates.
  	 * @param dir     The normalized(!) ray direction (corresponding to a specific pixel) in world coordinates [rd_x rd_y rd_z].
  	 * @param cubmin  The cuboid's minimal planes given as [min_x, min_y, min_z] in world
  	 *                coordinates.
  	 * @param cubmax  The cuboid's maximal planes given as [max_x, max_y, max_z] in world
  	 *                coordinates.
  	 * @param distanceNearAndFar    Return values. In case of a hit: Positive distances (in world
  	 *                coordinate units) of nearest and farthest plane intersection.
  	 * @return        Boolean value which is true if the ray intersects with the bounding
  	 *                box and false otherwise.
  	 *
  	 * @see <a href="http://www.siggraph.org/education/materials/HyperGraph/raytrace/rtinter0.htm">http://www.siggraph.org/education/materials/HyperGraph/raytrace/rtinter0.htm</a>
  	 * @see <a href="http://dx.doi.org/10.1145/15922.15916">http://dx.doi.org/10.1145/15922.15916</a>
  	 */
  	public static boolean intersectRayWithCuboid(final SimpleVector origin, final SimpleVector dir, final SimpleVector cubmin, final SimpleVector cubmax, final double[] distanceNearAndFar) {
  		// input checks
  		assert (origin.getLen() == 3) : new IllegalArgumentException("origin has to be a 3-vector!");
  		assert (dir.getLen() == 3) : new IllegalArgumentException("dir has to be a 3-vector!");
  		assert (Math.abs(dir.normL2() - 1.0) < 10.0*CONRAD.DOUBLE_EPSILON) : new IllegalArgumentException("Direction vector must be of length 1:" + dir);
  		assert (cubmin.getLen() == 3) : new IllegalArgumentException("cubmin has to be a 3-vector!");
  		assert (cubmax.getLen() == 3) : new IllegalArgumentException("cubmax has to be a 3-vector!");
  
  		// output checks
  		assert (distanceNearAndFar.length == 2) : new IllegalArgumentException("tntf output must be a 2-array!");
  
  		// init near and far intersection distance
  		distanceNearAndFar[0] = Double.NEGATIVE_INFINITY;
  		distanceNearAndFar[1] = Double.POSITIVE_INFINITY;
  
  		// intersect with bounding box's 6 planes
  		for (int i = 0; i < 3; ++i) { // loop over bb's x, y, and z plane pairs
  
  			// check whether ray is parallel to the two i planes
  			if (Math.abs(dir.getElement(i)) < CONRAD.DOUBLE_EPSILON) {
  
  				// check whether ray is not between the two i planes, i.e., if it misses the box
  				if (origin.getElement(i) < cubmin.getElement(i) || origin.getElement(i) > cubmax.getElement(i)) return false;
  
  			} else {
  
  				// compute the two intersections
  				double t1 = (cubmin.getElement(i) - origin.getElement(i))/dir.getElement(i);
  				double t2 = (cubmax.getElement(i) - origin.getElement(i))/dir.getElement(i);
  
  				// sort t1 and t2 so that t1 is intersection with near plane, t2 with far plane (t1 <= t2)
  				if (t1 > t2) {
  					// swap the two values if necessary
  					final double ttmp = t1;
  					t1 = t2;
  					t2 = ttmp;
  				}
  
  				// keep largest near intersection
  				if (t1 > distanceNearAndFar[0]) distanceNearAndFar[0] = t1;
  
  				// keep smallest far intersection
  				if (t2 < distanceNearAndFar[1]) distanceNearAndFar[1] = t2;
  
  				// check if box was missed
  				if (distanceNearAndFar[0] > distanceNearAndFar[1] + CONRAD.FLOAT_EPSILON) return false;
  
  				// check if box is behind the ray
  				if (distanceNearAndFar[1] < 0) return false;
  			}
  		}
  
  		// survived all tests, bounding box was hit by ray
  		return true;
  	}
  
  	/**
  	 * 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 geometric center of an iterator of points
  	 * @param list the set of points
  	 * @return the geometric center
  	 */
  	public static PointND getGeometricCenter(Iterator<PointND> list){
  		PointND p = list.next();
  		int count = 1;
  		int dim = p.getDimension();
  		double [] temp = new double [dim];
  		for (int j = 0; j < dim; j++){
  			temp[j] += p.get(j);
  		}
  		while (list.hasNext()){
  			p = list.next();
  			if (p!= null){
  				count ++;
  				for (int j = 0; j < dim; j++){
  					temp[j] += p.get(j);
  				}
  			}
  		}
  		for (int j = 0; j < dim; j++){
  			temp[j] /= count;
  		}
  		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;
  	}
  	
  	
  	/**
  	 * Threshold image and create binary mask
  	 * 
  	 * @param img the ImageProcessor
  	 * @param offset the threshold
  	 * @return the binary image
  	 */
  	public static ImageProcessor thresholdImage(ImageProcessor img, double offset){
  		byte[] pixels = new byte[img.getWidth()*img.getHeight()];
  		ImageProcessor imp = new BinaryProcessor(new ByteProcessor(img.getWidth(), img.getHeight(), pixels));
  		for (int j = 0; j< img.getHeight(); j++){
  			for (int i = 0; i< img.getWidth(); i++){
  				if (img.getPixelValue(i, j) > offset) {
  					imp.set(i,j,255);
  				}
  			}
  		}
  		return imp;
  	}
  	
  	
  
  	/**
  	 * 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){
  		return extractClusterCenter(pointList, distance, true);
  	}
  	
  	/**
  	 * Extracts cluster centers from an ordered List of points. 
  	 * A cluster center is then computed as the geometric center of the points in the same cluster.
  	 * @param pointList the list of candidate points
  	 * @param distance the minimal distance between clusters
  	 * @param breakOption option to break after distance exceeded once. Only applicable if points are ordered in the correct manner.
  	 * @return the list of cluster centers
  	 */
  	public static ArrayList<PointND> extractClusterCenter(ArrayList<PointND> pointList, double distance, boolean breakOption){
  		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 = pointList.size()-1; i >= 0; i--){
  				PointND current = pointList.get(i);
  				if (current.euclideanDistance(reference) < distance){
  					currentSubset.add(current);
  					pointList.remove(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 (breakOption) if (Math.abs(reference.get(1) - current.get(1)) > distance) break;
  				}
  			}
  			centerPoint.add(getGeometricCenter(currentSubset));
  		}
  		return centerPoint;
  	}
  
  	public static PointND getGeometricCenter(PointND[] pts) {
  		int dim = pts[0].getDimension();
  		double [] temp = new double [pts[0].getDimension()];
  		for (int i = 0; i < pts.length; i++){
  			for (int j = 0; j < dim; j++){
  				temp[j] += pts[i].get(j);
  			}
  		}
  		for (int j = 0; j < dim; j++){
  			temp[j] /= pts.length;
  		}
  		return new PointND(temp);
  	}
  	
  	/**
  	 * Creates a triangle mesh for a planar set of points. The geometric center is estimated and each subsequent set of points is connected with the center to form the mesh.
  	 * Note that the points must be neighboring points in the ArrayList. Such a set of points can be obtained using for example a convex hull algorithm. This method can also be used to close a BSplineSurface.
  	 * @param points the points
  	 * @return the list of triangles
  	 * @see ConvexHull
  	 * @see edu.stanford.rsl.conrad.geometry.splines.SurfaceBSpline
  	 */
  	public static ArrayList<Triangle> createTrianglesFromPlanarPointSet(ArrayList<PointND> points){
  		return createTrianglesFromPlanarPointSet(points, "", null);
  	}
  	
  	/**
  	 * Creates a triangle mesh for a planar set of points. The geometric center is estimated and each subsequent set of points is connected with the center to form the mesh.
  	 * Note that the points must be neighboring points in the ArrayList. Such a set of points can be obtained using for example a convex hull algorithm. This method can also be used to close a BSplineSurface.
  	 * @param points the points
  	 * @return the list of triangles
  	 * @see ConvexHull
  	 * @see edu.stanford.rsl.conrad.geometry.splines.SurfaceBSpline
  	 */
  	public static ArrayList<Triangle> createTrianglesFromPlanarPointSet(ArrayList<PointND> points, String nameTag, BufferedWriter bwpoint){
  		ArrayList<Triangle> mesh = new ArrayList<Triangle>();
  		if (points.size() > 0) {
  			PointND center = General.getGeometricCenter(points);
  			for (int i = 1; i < points.size(); i++){
  				try{
  					Triangle t = new Triangle(center, points.get(i), points.get(i-1));
  					if (bwpoint != null){
  						try {
  							PointND pt = center;
  							bwpoint.write("Center"+nameTag+"\t" + i+"id"+nameTag + "\t"+ +pt.get(0)+"\t"+pt.get(1)+"\t"+pt.get(2)+"\n");
  							pt = points.get(i);
  							bwpoint.write(nameTag+i+"\t"+i+"id"+nameTag + "\t" +pt.get(0)+"\t"+pt.get(1)+"\t"+pt.get(2)+"\n");
  							pt = points.get(i-1);
  							bwpoint.write(nameTag+(i-1)+"\t"+i+"id"+nameTag + "\t" +pt.get(0)+"\t"+pt.get(1)+"\t"+pt.get(2)+"\n");
  						} catch (IOException e) {
  							// TODO Auto-generated catch block
  							e.printStackTrace();
  						}
  					}
  					t.setName(i+"id"+nameTag);
  					mesh.add(t);
  				} catch (Exception e){
  					if (!e.getLocalizedMessage().contains("direction vector")){
  						System.out.println(e.getLocalizedMessage());
  					}
  				}
  			}
  			try{
  				Triangle t = new Triangle(center, points.get(0), points.get(points.size()-1));
  				if (bwpoint != null){
  					try {
  						PointND pt = center;
  						bwpoint.write("Center"+nameTag+"\t" + 0+"id"+nameTag + "\t"+ +pt.get(0)+"\t"+pt.get(1)+"\t"+pt.get(2)+"\n");
  						pt = points.get(0);
  						bwpoint.write(nameTag+0+"\t"+0+"id"+nameTag + "\t" +pt.get(0)+"\t"+pt.get(1)+"\t"+pt.get(2)+"\n");
  						pt = points.get(points.size()-1);
  						bwpoint.write(nameTag+(points.size()-1)+"\t"+0+"id"+nameTag + "\t" +pt.get(0)+"\t"+pt.get(1)+"\t"+pt.get(2)+"\n");
  					} catch (IOException e) {
  						// TODO Auto-generated catch block
  						e.printStackTrace();
  					}
  				}
  				t.setName(0+"id"+nameTag);
  				mesh.add(t);
  			} catch (Exception e){
  				if (!e.getLocalizedMessage().contains("direction vector")){
  					System.out.println(e.getLocalizedMessage());
  				}
  			}
  		}
  		return mesh;
  	}
  
  }
  /*
   * Copyright (C) 2010-2014 Andreas Keil
   * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
  */
  