/*--- formatted by Jindent 2.1, (www.c-lab.de/~jindent) ---*/

/*
 * The package jm.segmentation is part of the Renal Function Project
 * for analysis of dynamic contrast medium evalutions in MRT-Images
 * of the Kidneys.
 * 
 * Copyright (C) 1999 / 2000 Jens Martin
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License 2
 * as published by the Free Software Foundation.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 */

package jm.segmentation;

import ij.*;
import ij.gui.*;
import ij.plugin.filter.PlugInFilter;
import ij.process.*;
import java.awt.*;
import java.awt.image.*;
import java.util.*;


/**
 * Diese Klasse versucht, die Kontur der Niere durch ein Strahlverfolgungs-Verfahren
 * zu realisieren. Dabei wird das Intensitaetsprofil entlang des STrahls untersucht und
 * anhand definierte Intensitaetswechsel die Konturpunkte festgelegt. In einem zweiten
 * Durchgang werden die Konturpunkte auf ihre Plausiblitaet ueberprueft und anhand
 * von Kriterien korrigiert, die von einem Snake-Verfahren entliehen sind.
 * <BR> Der Aufrug erfolgt aus der Klasse AnalysisDialog
 * @see jm.kidney.AnalysisDialog
 * @version  6.53, 11/01/2000
 * @author   Jens Martin
 */
public class KidneyExtraction {
	private ImagePlus				curImage;
	private final double		fDegreeStep = 4.0;
	private Point						ptStart;
	private PolygonRoi			roi;
	private boolean					bSetRoi = true;
	private boolean					bMaskArea = false;
	private ShortProcessor	vertikalKernel = null;
	private ShortProcessor	diagonalKernel = null;
	private ShortProcessor	sobelKernel = null;
	private ShortProcessor	aktSobelIp = null;
	private ShortProcessor	aktIp = null;
	private short[]					aktPixel = null;
	private int							lastDeg;
	private double					lastLength;
	private int[]						xpoints = null;
	private int[]						ypoints = null;
	private int							numPoints = 0;
	private double[]				ptDistances = null;
	private double[]				ptAngles = null;
	private double[]				ptLengths = null;
	private double					midDistance = 0.0;
	private double					midAngle = 0.0;
	private double					midLength = 0.0;


	/**
	 * Der Konstruktor.
	 * @param    imp    ImagePlus-Objekt von ImageJ
	 */
	public KidneyExtraction(ImagePlus imp) {
		curImage = imp;
	}


	/**
	 * Berechnet die Winkel fuer ein Array von festgelegten Kontupunkten.
	 * Der Winkel in einem Konturpunkt wird jeweils durch den Punkt selbst sowie
	 * seinen Vorgaenger ung Nachfolger festgelegt.
	 * Zusaetzlich wird der mittlere Winkel aller Konturpunkte ermittelt.
	 */
	private void calcAllAngles() {
		midAngle = 0.0;
		ptAngles = new double[numPoints];
		for (int iPoint = 0; iPoint < numPoints; iPoint++) {
			int prevIndex = iPoint - 1;

			if (prevIndex < 0) {
				prevIndex = numPoints - 1;
			} 
			int			nextIndex = (iPoint + 1) % numPoints;
			double	g = calcAngle(xpoints[prevIndex], ypoints[prevIndex], xpoints[iPoint], ypoints[iPoint], xpoints[nextIndex], ypoints[nextIndex]);

			if ((ptLengths[prevIndex] >= ptLengths[iPoint]) && (ptLengths[nextIndex] >= ptLengths[iPoint])) {
				g = 360.0 - g;
			} 
			midAngle = midAngle + g;
			ptAngles[iPoint] = g;
		} 
		midAngle = midAngle / (double) numPoints;

		// IJ.write("mittlerer Winkel: " + midAngle);
	} 


	/**
	 * Berechnet fuer ein Array von festgelegten Kontupunkten deren Abstand zum
	 * Schwerpunkt der ROI. Der berechnete Abstand wird fuer jeden Punkt in ein Array
	 * eingetragen, das durch die Klassenvariable ptLengths referenziert ist.
	 * Zusaetzlich wird der mittlere Abstand aller Konturpunkte ermittelt.
	 */
	private void calcAllLengths() {
		midLength = 0.0;
		ptLengths = new double[numPoints];
		for (int iPoint = 0; iPoint < numPoints; iPoint++) {
			double	d = calcDistance(xpoints[iPoint], ypoints[iPoint], ptStart.x, ptStart.y);

			midLength = midLength + d;
			ptLengths[iPoint] = d;
		} 
		midLength = midLength / (double) numPoints;
		IJ.write("mittlere Laenge: " + midLength);
	} 


	/**
	 * Berechnet fuer ein Array von festgelegten Kontupunkten deren Abstand zum
	 * jeweiligen Vorgaengerpunkt (im Uhrzeigersinn). Der berechnete Abstand wird
	 * fuer jeden Punkt in ein Array
	 * eingetragen, das durch die Klassenvariable ptDistances referenziert ist.
	 * Zusaetzlich wird die mittlere Distanz ermittelt.
	 */
	private void calcAllDistances() {
		midDistance = 0.0;
		ptDistances = new double[numPoints];
		for (int iPoint = 0; iPoint < numPoints; iPoint++) {
			int prevIndex = iPoint - 1;

			if (prevIndex < 0) {
				prevIndex = numPoints + prevIndex;
			} 
			double	d = calcDistance(xpoints[iPoint], ypoints[iPoint], xpoints[prevIndex], ypoints[prevIndex]);

			midDistance = midDistance + d;
			ptDistances[iPoint] = d;
		} 
		midDistance = midDistance / (double) numPoints;
		IJ.write("mittlere Distanz: " + midDistance);
	} 


	/**
	 * Hilfsmethode. Rechnet einen uebergebenen Winkel von Winkelmass in
	 * das Radialmass um.
	 * @param f Winkel (0-359°)
	 * @return double Winkel als radailes Mass
	 */
	private double degToRad(float f) {
		return f * Math.PI / 180.0;
	} 


	/**
	 * Hilfsmethode. Rechnet einen uebergebenen Winkel von Radialmass in
	 * das Winkelmass um.
	 * @param f Winkel im radialen Mass
	 * @return double Winkel im Winkelmass (0-359°)
	 */
	private double radToDeg(double f) {
		return f * 180.0 / Math.PI;
	} 


	/**
	 * Hilfsmethode. Brechnet die Distanz zwischen zwei Rasterpunkten.
	 * @param pointX1 die x-Koordinate der ersten Punktes (ganzzahlig)
	 * @param pointY1 die y-Koordinate der ersten Punktes (ganzzahlig)
	 * @param pointX2 die x-Koordinate der zweiten Punktes (ganzzahlig)
	 * @param pointY2 die y-Koordinate der zweiten Punktes (ganzzahlig)
	 * @return double die Distanz zwischen den Punkten
	 */
	private double calcDistance(int pointX1, int pointY1, int pointX2, int pointY2) {
		return Math.sqrt((double) ((pointX1 - pointX2) * (pointX1 - pointX2)) + (double) ((pointY1 - pointY2) * (pointY1 - pointY2)));
	} 


	/**
	 * Hilfsmethode. Brechnet den Winkel der Nierenkontur in einem Punkt.
	 * Es muessen drei Punkte uebergeben werden, die den Winkel definieren.
	 * @param pointX1 die x-Koordinate der ersten Punktes (ganzzahlig)
	 * @param pointY1 die y-Koordinate der ersten Punktes (ganzzahlig)
	 * @param pointX2 die x-Koordinate der zweiten Punktes (ganzzahlig)
	 * @param pointY2 die y-Koordinate der zweiten Punktes (ganzzahlig)
	 * @param pointX3 die x-Koordinate der dritten Punktes (ganzzahlig)
	 * @param pointY3 die y-Koordinate der dritten Punktes (ganzzahlig)
	 * @return double die Distanz zwischen den Punkten
	 */
	private double calcAngle(int pointX1, int pointY1, int pointX2, int pointY2, int pointX3, int pointY3) {

		// double c = Math.sqrt(  (double)((xpoints[ nextIndex] - xpoints[ prevIndex]) * (xpoints[ nextIndex] - xpoints[ prevIndex])) +
		// (double)((ypoints[ nextIndex] - ypoints[ prevIndex]) * (ypoints[ nextIndex] - ypoints[ prevIndex]))  );
		double	c = calcDistance(pointX1, pointY1, pointX3, pointY3);

		// double a = ptDistances[iPoint];
		double	a = calcDistance(pointX1, pointY1, pointX2, pointY2);

		// double b = ptDistances[nextIndex];
		double	b = calcDistance(pointX2, pointY2, pointX3, pointY3);
		double	g = 180.0;

		if (!(a == 0.0 || b == 0.0)) {
			double	cosg = ((a * a) + (b * b) - (c * c)) / (2.0 * a * b);

			g = radToDeg(Math.acos(cosg));
		} 
		return g;
	} 


	/**
	 * Wendet einen Medialfilter zur Glaettung auf das uebergebene Bild an.
	 * Es wird ein ImageJ-ImageProcessor uebergeben, um einfacheren
	 * Zugriff auf die Pixelwerte zu haben.
	 * @param ip ImageProcessor des Eingabebildes (ImageJ-Datentyp)
	 * @return ImageProcessor das sobel-gefilterte Bild
	 */
	private ImageProcessor medianFilter(ImageProcessor ip) {
		int							width = ip.getWidth();
		int							height = ip.getHeight();

		ImageProcessor	result = null;
		short						median[] = new short[width * height];
		short						pix[] = (short[]) ip.getPixels();

		int[]						p = new int[10];
		int							tmp = 0;
		int							i, x, y;

		for (y = 1; y < height - 1; y++) {
			for (x = 1; x < width - 1; x++) {
				if (x == 1) {
					p[1] = pix[((y - 1) * width) + (x - 1)];	// img.get(xx-1,yy-1);
					p[2] = pix[((y - 1) * width) + (x)];			// img.get(xx  ,yy-1);
					p[3] = pix[((y - 1) * width) + (x + 1)];	// img.get(xx+1,yy-1);
					p[4] = pix[((y) * width) + (x - 1)];			// img.get(xx-1,yy  );
					p[5] = pix[((y) * width) + (x)];					// img.get(xx  ,yy  );
					p[6] = pix[((y) * width) + (x + 1)];			// img.get(xx+1,yy  );
					p[7] = pix[((y) * width) + (x - 1)];			// img.get(xx-1,yy+1);
					p[8] = pix[((y + 1) * width) + (x)];			// img.get(xx  ,yy+1);
					p[9] = pix[((y + 1) * width) + (x + 1)];	// img.get(xx+1,yy+1);
				} else {
					p[1] = p[2];
					p[2] = p[3];
					p[3] = pix[((y - 1) * width) + (x + 1)];	// img.get(xx+1,yy-1);
					p[4] = p[5];
					p[5] = p[6];
					p[6] = pix[((y) * width) + (x + 1)];			// img.get(xx+1,yy  );
					p[7] = p[8];
					p[8] = p[9];
					p[9] = pix[((y + 1) * width) + (x + 1)];	// img.get(xx+1,yy+1);
				} 
				for (i = 1; i < 9; i++) {
					if (p[i] > p[i + 1]) {
						tmp = p[i];
						p[i] = p[i + 1];
						p[i + 1] = tmp;
					} 
				} 
				median[((y) * width) + (x)] = (short) p[4];
			} 
		} 
		result = new ShortProcessor(width, height, median, ip.getColorModel(), false);
		return result;
	} 


	/**
	 * Hauptberechnungsklasse dieser Methode. Hier werden im ersten Durchlauf
	 * die Konturpunkte berechnet und im zweiten Durchlauf auf ihre Gueltigkeit geprueft.
	 * @param ip ImageProcessor des Eingabebildes (ImageJ-Datentyp)
	 */
	public void doSegmentation(ImageProcessor ip) {
		try {

			// dimensions
			int width = ip.getWidth();
			int height = ip.getHeight();
			int i;

			// ip = medianFilter(ip);
			aktIp = (ShortProcessor) medianFilter(ip);
			aktPixel = (short[]) aktIp.getPixels();

			initKernels();

			Roi				roi = curImage.getRoi();
			Rectangle rr = roi.getBoundingRect();

			// startpunkt
			ptStart = new Point();
			ptStart.x = rr.x + rr.width / 2;
			ptStart.y = rr.y + rr.height / 2;

			double	startDegree = 0.0;

			if (ptStart.x < width / 2) {
				startDegree = 180.0;

			} 
			double	maxSteig = (double) (rr.y - ptStart.y) / (double) ((rr.x + rr.width) - ptStart.x);

			if (maxSteig < 0) {
				maxSteig = maxSteig * (-1.0);

				// 1x rumdrehen
			} 
			numPoints = (int) ((double) 360 / fDegreeStep);

			xpoints = new int[numPoints];
			ypoints = new int[numPoints];
			double	degreeValues[] = new double[numPoints];

			int			iPoint = 0;

			float		fDegree = (float) startDegree;

			// while( fDegree < 360)
			// erster Durchlauf
			while (iPoint < numPoints) {
				double					rad = degToRad(fDegree);
				double					si = Math.sin((double) rad);
				double					cs = Math.cos((double) rad);

				double					len = Math.max(rr.height, rr.width) / 2;

				int							xd = (int) (len * cs);
				int							yd = (int) (len * si);
				Point						ptEnd = new Point(ptStart.x + xd, ptStart.y + yd);

				double[]				line1 = aktIp.getLine(ptStart.x, ptStart.y, ptEnd.x, ptEnd.y);

				// ScanLine sl = new ScanLine( line);

				ShortProcessor	sIp = getRotateSobel((double) fDegree);
				double[]				line2 = sIp.getLine(ptStart.x, ptStart.y, ptEnd.x, ptEnd.y);

				// ScanLine sobelsl = new ScanLine( line2);

				DoubleLine			sl = new DoubleLine(line1, line2);

				// smoothen
				sl.smoothLines(3);

				int			endll = rr.height / 2;

				// int le = rr.width / 2;
				double	steigung;
				int			d = (ptEnd.x - ptStart.x);

				if (d != 0) {
					steigung = (double) (ptEnd.y - ptStart.y) / (double) d;
				} else {
					steigung = Double.MAX_VALUE;
				} 
				if ((steigung >= (-1.0)) && (steigung <= 1.0)) {
					endll = rr.width / 2;

					// le = rr.height / 2;
				} 

				// hier die line-Arrays abgrasen.
				int startll = 3 * (endll / 5);
				int ll = sl.findBestPoint(startll, endll);

				// ll ist "ziel"
				int xd2 = (int) (((double) ll) * cs);
				int yd2 = (int) (((double) ll) * si);

				xpoints[iPoint] = ptStart.x + xd2;
				ypoints[iPoint] = ptStart.y + yd2;
				degreeValues[iPoint] = fDegree;

				fDegree += fDegreeStep;
				if (fDegree >= 360.0) {
					fDegree = (float) (fDegree % 360);
				} 
				iPoint++;
			} 

			// zweiter Durchlauf
			for (int grosseSchleife = 1; grosseSchleife <= 1; grosseSchleife++) {
				java.awt.Polygon	poly = new java.awt.Polygon(xpoints, ypoints, numPoints);

				// 2000.8.3 tha: In ImageJ 1.17 wird nicht mehr per boolean entschieden,
				// ob die Koordinaten per WandTool erzeugt wurden. Wir zeichnen
				// Polygone.
				// roi = new PolygonRoi( poly.xpoints, poly.ypoints, poly.npoints, curImage,false);
				roi = new PolygonRoi(poly.xpoints, poly.ypoints, poly.npoints, curImage, Roi.POLYGON);
				curImage.setRoi(roi);
				curImage.updateAndRepaintWindow();

				calcAllDistances();

				calcAllLengths();

				calcAllAngles();

				for (int index = 0; index < numPoints; index++) {
					int prevIndex = index - 1;

					if (prevIndex < 0) {
						prevIndex = numPoints + prevIndex;
					} 
					int			nextIndex = (index + 1) % numPoints;
					double	maxPercentDev = 20.0;

					// double d = calcDistance(xpoints[ iPoint],ypoints[ iPoint],ptStart.x,ptStart.y);
					double	diffPercentPrev = Math.abs(ptLengths[prevIndex] - ptLengths[index]) * 100.0 / ptLengths[prevIndex];
					double	diffPercentNext = Math.abs(ptLengths[nextIndex] - ptLengths[index]) * 100.0 / ptLengths[nextIndex];

					if ((Math.abs(ptAngles[index] - 180.0) > 130.0) && ((diffPercentPrev > maxPercentDev) || (diffPercentNext > maxPercentDev))) {
						int p1;

						p1 = ip.getPixel(xpoints[index], ypoints[index]);
						ip.putPixel(xpoints[index], ypoints[index], 255);
						curImage.updateAndRepaintWindow();

						// IJ.error ("siehe Pixel:");
						ip.putPixel(xpoints[index], ypoints[index], p1);

						double	d = calcDistance(xpoints[prevIndex], ypoints[prevIndex], xpoints[nextIndex], ypoints[nextIndex]);

						i = 0;

						// pos = (ptBadStart + i) % numPoints;
						double	xvecC = (double) (xpoints[nextIndex] - xpoints[prevIndex]);
						double	yvecC = (double) (ypoints[nextIndex] - ypoints[prevIndex]);
						double	partofd = d / 2.0;
						int			newposx = xpoints[prevIndex] + (int) (partofd / d * xvecC);
						int			newposy = ypoints[prevIndex] + (int) (partofd / d * yvecC);

						// IJ.error("-----\nstartPos: " + xpoints[startPos] + "/" + ypoints[startPos] + "\niPoint:" + xpoints[iPoint] + "/" + ypoints[iPoint] + "\nVectorC: " + xvecC + "/" + yvecC + "\nd: " + d + "\npartofd: " + partofd + "\ni: " + i + "\nptBadAnz: " + ptBadAnz + "\nnewpos: " + newposx + "/" + newposy);

						xpoints[index] = newposx;		// points[startPos] + (int)( (double)(i+1) * (d / ((double)(ptBadAnz + 1)) ) * ( (double)(xpoints[iPoint] - xpoints[startPos])));
						ypoints[index] = newposy;		// ypoints[startPos] + (int)( (double)(i+1) * (d / ((double)(ptBadAnz + 1)) ) * ( (double)(ypoints[iPoint] - ypoints[startPos])));

						// ptLengths[pos] = calcDistance(xpoints[ pos],ypoints[ pos],ptStart.x,ptStart.y);

					} 
				} 

				IJ.write("Suche laengste Strecke...");
				int		pos = 0;
				int[] bestStart = new int[numPoints];
				int[] bestEnd = new int[numPoints];

				for (iPoint = 1; iPoint < numPoints; iPoint++) {
					bestStart[pos] = iPoint - 1;
					bestEnd[pos] = iPoint - 1;

					double	d = ptDistances[iPoint];
					int			counter = 0;

					while ((d <= midDistance) && (counter < numPoints)) {
						bestEnd[pos] = (iPoint + counter) % numPoints;

						// iPoint++;
						counter++;
						int arrayPos = (iPoint + counter) % numPoints;

						d = ptDistances[arrayPos];

					} 

					// IJ.write("Einstieg in While mit iPoint " + iPoint + ",Laenge Segment: " + counter + ", best Start/End: " + bestStart[pos] + "/" + bestEnd[pos]);
					iPoint = iPoint + counter;
					pos++;
				} 

				IJ.write("Suche besten Positionen...");
				int alen;
				int maxLength = 0;
				int bestPosStart = 0;
				int bestPosEnd = 0;

				for (i = 0; i < pos; i++) {

					// IJ.write("best Start/End: " + bestStart[i] + "/" + bestEnd[i]);

					if (bestStart[i] <= bestEnd[i]) {
						alen = bestEnd[i] - bestStart[i];
					} else {
						alen = bestEnd[i] + numPoints - bestStart[i];
					} 

					if (alen > maxLength) {
						maxLength = alen;
						bestPosStart = bestStart[i];
						bestPosEnd = bestEnd[i];
					} 
				} 
				if (maxLength < 5) {
					IJ.error("Segmentation of the kidneys was not succesful.\nPlease select the ROI manually !!!");
				} 
				IJ.write("Start/End/numPoints: " + bestPosStart + "/" + bestPosEnd + "/" + numPoints);

				// IJ.write("und noch kurz einzeichnen...");
				// ImageProcessor ip = curImage.getProcessor();
				iPoint = bestPosStart;
				pos = 0;
				while (pos <= maxLength) {

					// ip.putPixel(xpoints[ iPoint],ypoints[ iPoint],255);
					iPoint = (iPoint + 1) % numPoints;
					pos++;

					// newIp2.putPixel(xpoints[ arrayPos],ypoints[ arrayPos],255);
				} 

				// ab hier der zweite Durchgang...
				boolean[] ptBad = new boolean[numPoints];

				for (i = 0; i < numPoints; i++) {
					ptBad[i] = false;
				} 

				iPoint = bestPosEnd;

				// iPoint = 200000;
				startDegree = degreeValues[bestPosEnd];
				fDegree = (float) startDegree;

				// IJ.write("...fertig");
				int endPos = bestPosStart - 1;

				if (endPos < 0) {
					endPos = numPoints - endPos;

					// iPoint = endPos;
				} 
				while (iPoint != endPos) {
					poly = new java.awt.Polygon(xpoints, ypoints, numPoints);

					// 2000.8.3 tha: In ImageJ 1.17 wird nicht mehr per boolean entschieden,
					// ob die Koordinaten per WandTool erzeugt wurden. Wir zeichnen
					// Polygone.
					// roi = new PolygonRoi( poly.xpoints, poly.ypoints, poly.npoints, curImage,false);
					roi = new PolygonRoi(poly.xpoints, poly.ypoints, poly.npoints, curImage, Roi.POLYGON);

					curImage.setRoi(roi);
					curImage.updateAndRepaintWindow();

					double	midLenLast5 = 0.0;
					double	midWinkelLast5 = 0.0;

					for (pos = 0; pos < 5; pos++) {
						int arrayPos = iPoint - pos;

						if (arrayPos < 0) {
							arrayPos = numPoints + arrayPos;
						} 
						midLenLast5 = midLenLast5 + ptLengths[arrayPos];
						midWinkelLast5 = midWinkelLast5 + ptAngles[arrayPos];

					} 
					midLenLast5 = midLenLast5 / 5.0;
					midWinkelLast5 = midWinkelLast5 / 5.0;

					iPoint = (iPoint + 1) % numPoints;

					fDegree = (float) degreeValues[iPoint];

					// IJ.write("Zweite Druchgang: Winkel = " + fDegree + ", MidLenLast5: " + midLenLast5 + ", MidWineklLast5: " + midWinkelLast5);

					double	rad = degToRad(fDegree);
					double	si = Math.sin((double) rad);
					double	cs = Math.cos((double) rad);
					double	len = Math.max(rr.height, rr.width) / 2;
					int			xd = (int) (len * cs);
					int			yd = (int) (len * si);
					Point		ptEnd = new Point(ptStart.x + xd, ptStart.y + yd);
					int			xgood = 0;
					int			ygood = 0;
					double	factor = 1.4;
					double	maxAngleDev = 45.0;
					double	maxLengthPercentDev = 25.0;
					int			ptBadStart = 0;
					int			ptBadAnz = 0;
					double	d = calcDistance(xpoints[iPoint], ypoints[iPoint], ptStart.x, ptStart.y);
					double	percentDev = Math.abs(midLenLast5 - d) * 100.0 / midLenLast5;
					int			prev = iPoint - 1;

					if (prev < 0) {
						prev = numPoints + prev;
					} 
					int prev2 = iPoint - 2;

					if (prev2 < 0) {
						prev = numPoints + prev2;
					} 
					int prev3 = iPoint - 3;

					if (prev3 < 0) {
						prev = numPoints + prev3;
					} 
					double	prevw = calcAngle(xpoints[prev3], ypoints[prev3], xpoints[prev2], ypoints[prev2], xpoints[prev], ypoints[prev]);

					if ((ptLengths[prev3] >= ptLengths[prev2]) && (ptLengths[prev] >= ptLengths[prev2])) {
						prevw = 360.0 - prevw;

						// IJ.write("nu aber");
					} 
					double	w = calcAngle(xpoints[iPoint], ypoints[iPoint], xpoints[prev], ypoints[prev], xpoints[prev2], ypoints[prev2]);

					if ((ptLengths[iPoint] >= ptLengths[prev]) && (ptLengths[prev2] >= ptLengths[prev])) {
						int p1, p2, p3, p4;

						p1 = ip.getPixel(xpoints[iPoint], ypoints[iPoint]);
						ip.putPixel(xpoints[iPoint], ypoints[iPoint], 255);
						p2 = ip.getPixel(xpoints[prev], ypoints[prev]);
						ip.putPixel(xpoints[prev], ypoints[prev], 255);
						p3 = ip.getPixel(xpoints[prev2], ypoints[prev2]);
						ip.putPixel(xpoints[prev2], ypoints[prev2], 255);
						p4 = ip.getPixel(ptStart.x, ptStart.y);
						ip.putPixel(ptStart.x, ptStart.y, 255);
						curImage.updateAndRepaintWindow();

						// IJ.error ("siehe Pixel:\n" + ptLengths[iPoint] + "\n" + ptLengths[prev] + "\n" + ptLengths[prev2]);
						ip.putPixel(xpoints[iPoint], ypoints[iPoint], p1);
						ip.putPixel(xpoints[prev], ypoints[prev], p2);
						ip.putPixel(xpoints[prev2], ypoints[prev2], p3);
						ip.putPixel(ptStart.x, ptStart.y, p4);
						curImage.updateAndRepaintWindow();
						w = 360.0 - w;

						// IJ.write("nu aber");
					} 
					w = w + ((midWinkelLast5 - prevw) / 2.0);
					double	wDev = Math.abs(midWinkelLast5 - w);

					if (((ptDistances[iPoint] > (midDistance * factor)) && (wDev > maxAngleDev)) || (percentDev > maxLengthPercentDev)) {
						IJ.write("__________________\n  Nummer der ersten ptBad: " + iPoint);
						IJ.write("Winkel neu: " + w + "\n  Winkel Deviation: " + wDev + "\n  Mittlerer Winkel Last5: " + midWinkelLast5);
						IJ.write("Winkel alt: " + (w - ((midWinkelLast5 - prevw) / 2.0)) + "\n  Laenge Last5:" + midLenLast5 + "\n  Prozent Laegenabweichung: " + percentDev + "\n  Voriger Winkel: " + prevw);
						IJ.write("Mittlere Distanz: " + midDistance + "\n  midDistance * factor: " + (midDistance * factor));
						ptBadStart = iPoint;
						ptBadAnz++;
						ptBad[iPoint] = true;

						// ip.putPixel(xpoints[ iPoint],ypoints[ iPoint],255);
						curImage.updateAndRepaintWindow();
						iPoint = (iPoint + 1) % numPoints;
						d = calcDistance(xpoints[iPoint], ypoints[iPoint], ptStart.x, ptStart.y);
						percentDev = Math.abs(midLenLast5 - d) * 100.0 / midLenLast5;

						// prev = iPoint - 1;
						// if (prev < 0) prev = numPoints + prev;
						// prev2 = iPoint - 1;
						// if (prev2 < 0) prev = numPoints + prev2;
						w = calcAngle(xpoints[iPoint], ypoints[iPoint], xpoints[prev], ypoints[prev], xpoints[prev2], ypoints[prev2]);
						if ((ptLengths[iPoint] >= ptLengths[prev]) && (ptLengths[prev2] >= ptLengths[prev])) {
							w = 360.0 - w;

							// IJ.write("nu aber");
						} 

						// IJ.write("w/wDev: " + w + "/" + wDev + "/" + midWinkelLast5);
						w = w + ((midWinkelLast5 - prevw) / 2.0);
						wDev = Math.abs(midWinkelLast5 - w);

						// IJ.write("w/wDev: " + w + "/" + wDev + "/" + midWinkelLast5);

						while ( /* (ptBadAnz < 4) && */(ptDistances[iPoint] <= (midDistance * factor)) && (wDev > maxAngleDev) && (iPoint != endPos) && (percentDev > maxLengthPercentDev)) {
							IJ.write("(o)w: " + w + " wDev: " + wDev + " mw5: " + midWinkelLast5);
							IJ.write("(n)w:" + (w - ((midWinkelLast5 - prevw) / 2.0)) + " %d: " + percentDev + " pw: " + prevw);
							ptBadAnz++;

							// IJ.error("Next Point...");
							ptBad[iPoint] = true;

							// ip.putPixel(xpoints[ iPoint],ypoints[ iPoint],255);
							// curImage.updateAndRepaintWindow();
							iPoint = (iPoint + 1) % numPoints;

							d = calcDistance(xpoints[iPoint], ypoints[iPoint], ptStart.x, ptStart.y);
							percentDev = Math.abs(midLenLast5 - d) * 100.0 / midLenLast5;

							// prev = iPoint - 1;
							// if (prev < 0) prev = numPoints + prev;
							// prev2 = iPoint - 1;
							// if (prev2 < 0) prev = numPoints + prev2;
							w = calcAngle(xpoints[iPoint], ypoints[iPoint], xpoints[prev], ypoints[prev], xpoints[prev2], ypoints[prev2]);
							if ((ptLengths[iPoint] >= ptLengths[prev]) || (ptLengths[prev] >= ptLengths[prev2])) {
								w = 360.0 - w;

								// IJ.write("nu aber");
							} 
							w = w + ((midWinkelLast5 - prevw) / 2.0);
							wDev = Math.abs(midWinkelLast5 - w);

							// IJ.write("w/wDev: " + w + "/" + wDev + "/" + midWinkelLast5);
						} 
					} 

					if (ptBadAnz > 0) {
						int startPos = ptBadStart - 1;

						if (startPos < 0) {
							startPos = numPoints - startPos;
						} 
						d = calcDistance(xpoints[startPos], ypoints[startPos], xpoints[iPoint], ypoints[iPoint]);
						i = 0;
						while (i < ptBadAnz) {
							pos = (ptBadStart + i) % numPoints;
							double	xvecC = (double) (xpoints[iPoint] - xpoints[startPos]);
							double	yvecC = (double) (ypoints[iPoint] - ypoints[startPos]);
							double	partofd = (double) (i + 1) * d / ((double) (ptBadAnz + 1));
							int			newposx = xpoints[startPos] + (int) (partofd / d * xvecC);
							int			newposy = ypoints[startPos] + (int) (partofd / d * yvecC);

							// IJ.error("-----\nstartPos: " + xpoints[startPos] + "/" + ypoints[startPos] + "\niPoint:" + xpoints[iPoint] + "/" + ypoints[iPoint] + "\nVectorC: " + xvecC + "/" + yvecC + "\nd: " + d + "\npartofd: " + partofd + "\ni: " + i + "\nptBadAnz: " + ptBadAnz + "\nnewpos: " + newposx + "/" + newposy);

							xpoints[pos] = newposx;		// points[startPos] + (int)( (double)(i+1) * (d / ((double)(ptBadAnz + 1)) ) * ( (double)(xpoints[iPoint] - xpoints[startPos])));
							ypoints[pos] = newposy;		// ypoints[startPos] + (int)( (double)(i+1) * (d / ((double)(ptBadAnz + 1)) ) * ( (double)(ypoints[iPoint] - ypoints[startPos])));
							ptLengths[pos] = calcDistance(xpoints[pos], ypoints[pos], ptStart.x, ptStart.y);
							prev = pos - 1;
							if (prev < 0) {
								prev = numPoints + prev;
							} 
							ptDistances[pos] = calcDistance(xpoints[prev], ypoints[prev], xpoints[pos], ypoints[pos]);
							int prevIndex = pos - 1;

							if (prevIndex < 0) {
								prevIndex = numPoints + prevIndex;
							} 
							int prevprevIndex = pos - 2;

							if (prevprevIndex < 0) {
								prevprevIndex = numPoints + prevprevIndex;
							} 
							int			nextIndex = (pos + 1) % numPoints;
							int			nextnextIndex = (pos + 2) % numPoints;
							double	g;

							g = calcAngle(xpoints[prevprevIndex], ypoints[prevprevIndex], xpoints[prevIndex], ypoints[prevIndex], xpoints[pos], ypoints[pos]);
							if ((ptLengths[prevprevIndex] >= ptLengths[prevIndex]) || (ptLengths[prevIndex] >= ptLengths[pos])) {
								g = 360.0 - g;
							} 
							ptAngles[prevIndex] = g;

							g = calcAngle(xpoints[prevIndex], ypoints[prevIndex], xpoints[pos], ypoints[pos], xpoints[nextIndex], ypoints[nextIndex]);
							if ((ptLengths[prevIndex] >= ptLengths[pos]) || (ptLengths[pos] >= ptLengths[nextIndex])) {
								g = 360.0 - g;
							} 
							ptAngles[pos] = g;

							g = calcAngle(xpoints[pos], ypoints[pos], xpoints[nextIndex], ypoints[nextIndex], xpoints[nextnextIndex], ypoints[nextnextIndex]);
							if ((ptLengths[pos] >= ptLengths[nextIndex]) || (ptLengths[nextIndex] >= ptLengths[nextnextIndex])) {
								g = 360.0 - g;
							} 
							ptAngles[nextIndex] = g;

							// ip.putPixel(xpoints[ pos],ypoints[ pos],255);
							poly = new java.awt.Polygon(xpoints, ypoints, numPoints);

							// 2000.8.3 tha: In ImageJ 1.17 wird nicht mehr per boolean entschieden,
							// ob die Koordinaten per WandTool erzeugt wurden. Wir zeichnen
							// Polygone.
							// roi = new PolygonRoi( poly.xpoints, poly.ypoints, poly.npoints, curImage,false);
							roi = new PolygonRoi(poly.xpoints, poly.ypoints, poly.npoints, curImage, Roi.POLYGON);

							curImage.setRoi(roi);
							curImage.updateAndRepaintWindow();
							i++;
						} 
					} 
				} 
			} 


			// in Polygon-Objekt umsetzen
			java.awt.Polygon	poly_org = new java.awt.Polygon(xpoints, ypoints, numPoints);

			// reduzieren...
			// java.awt.Polygon poly = smoothPolygon( poly_org);
			java.awt.Polygon	poly = poly_org;

			// java.awt.Polygon poly = reducePolygon( poly2);
			// java.awt.Polygon poly = poly_org;

			if (bSetRoi == true) {

				// 2000.8.3 tha: In ImageJ 1.17 wird nicht mehr per boolean entschieden,
				// ob die Koordinaten per WandTool erzeugt wurden. Wir zeichnen
				// Polygone.
				// roi = new PolygonRoi( poly.xpoints, poly.ypoints, poly.npoints, curImage,false);
				roi = new PolygonRoi(poly.xpoints, poly.ypoints, poly.npoints, curImage, Roi.POLYGON);

				curImage.setRoi(roi);

				curImage.draw();
			} 

		} catch (Exception ex0) {
			ex0.printStackTrace(System.out);
		} 
	} 


	/**
	 * Hilfsmethode. Glaettet die berechnete Polygon-ROI, der die Nierenkontur
	 * markiert.
	 * @param p Eingabe der ROI als java.awt.Polygon-Datentyp
	 * @return java.awt.Polygon das geglaettete Polygon
	 */
	private java.awt.Polygon smoothPolygon(java.awt.Polygon p) {
		int i = 1;
		int n = p.npoints;
		int r = 1;
		int xp[] = new int[n];
		int yp[] = new int[n];

		xp[0] = p.xpoints[0];
		yp[0] = p.ypoints[0];

		while (i < n - 1) {

			// betrachte dreieck der Punkte i-1, i, i+1
			int x1 = p.xpoints[i - 1];
			int y1 = p.ypoints[i - 1];
			int x2 = p.xpoints[i];
			int y2 = p.ypoints[i];
			int x3 = p.xpoints[i + 1];
			int y3 = p.ypoints[i + 1];

			// flaecheninhalt
			// smoothing
			int x2a = (x1 + x3) / 2;
			int y2a = (y1 + y3) / 2;
			int x2b = (x2 + x2a) / 2;
			int y2b = (y2 + y2a) / 2;

			xp[r] = x2b;
			yp[r] = y2b;

			r++;
			i++;
		} 
		xp[n - 1] = p.xpoints[n - 1];
		yp[n - 1] = p.ypoints[n - 1];
		r++;

		return new java.awt.Polygon(xp, yp, r);
	} 


	/**
	 * Hilfsmethode. Reduziert die Anzahl der Punkte im Polygon. Redundante Punkte
	 * sind solche, die genau auf der Strecke zwischen ihren zwei Nachbarpunkten liegen.
	 * @param p Eingabe der ROI als java.awt.Polygon-Datentyp
	 * @return java.awt.Polygon das geglaettete Polygon
	 */
	private java.awt.Polygon reducePolygon(java.awt.Polygon p) {
		int i = 1;
		int n = p.npoints;
		int r = 1;
		int xp[] = new int[n];
		int yp[] = new int[n];

		xp[0] = p.xpoints[0];
		yp[0] = p.ypoints[0];

		int x1 = p.xpoints[i - 1];
		int y1 = p.ypoints[i - 1];

		while (i < n - 1) {

			// betrachte dreieck der Punkte i-1, i, i+1
			int			x2 = p.xpoints[i];
			int			y2 = p.ypoints[i];
			int			x3 = p.xpoints[i + 1];
			int			y3 = p.ypoints[i + 1];

			// flaecheninhalt
			double	a = 0.5 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2));

			// wenn flaecheninhalt < x, dann punkte vergessen
			boolean bDiscard = false;

			if (Math.abs(a) < 6) {
				bDiscard = true;

			} 
			if (bDiscard) {
				i++;
			} else {
				xp[r] = x2;
				yp[r] = y2;

				r++;
				i++;
				x1 = x2;
				y1 = y2;
			} 
		} 

		// System.out.println( "Reducing: "+r+" / "+n);
		return new java.awt.Polygon(xp, yp, r);
	} 


	/**
	 * Hilfsmethode. Initalisiert die Kernel, die zur schnellen Berechnung des Kompass-Gradienten
	 * benoetigt werden. Es wird fuer jede Himmelsrichtung ein Kernel in einem eigenen Array angelegt.
	 */
	private void initKernels() {
		vertikalKernel = new ShortProcessor(5, 5, false);
		diagonalKernel = new ShortProcessor(5, 5, false);
		int x, y;

		for (x = 0; x < 5; x++) {
			for (y = 0; y < 5; y++) {
				diagonalKernel.putPixel(x, y, 0);
				vertikalKernel.putPixel(x, y, 0);
			} 
		} 

		// init the vertikal Kernel
		for (x = 0; x < 2; x++) {
			for (y = 0; y < 5; y++) {
				vertikalKernel.putPixel(x, y, 1);
			} 
		} 
		vertikalKernel.putPixel(1, 2, 2);

		for (x = 3; x < 5; x++) {
			for (y = 0; y < 5; y++) {
				vertikalKernel.putPixel(x, y, -1);
			} 
		} 
		vertikalKernel.putPixel(3, 2, -2);

		// init the diagonal Kernel
		y = 0;
		for (x = 0; x < 4; x++) {
			diagonalKernel.putPixel(x, y, 1);
		} 
		y = 1;
		for (x = 0; x < 3; x++) {
			diagonalKernel.putPixel(x, y, 1);
		} 
		diagonalKernel.putPixel(0, 2, 1);
		diagonalKernel.putPixel(1, 2, 1);
		diagonalKernel.putPixel(0, 3, 1);

		y = 4;
		for (x = 1; x < 5; x++) {
			diagonalKernel.putPixel(x, y, -1);
		} 
		y = 3;
		for (x = 2; x < 5; x++) {
			diagonalKernel.putPixel(x, y, -1);
		} 
		diagonalKernel.putPixel(4, 2, -1);
		diagonalKernel.putPixel(3, 2, -1);
		diagonalKernel.putPixel(4, 1, -1);

	} 


	/**
	 * Hilfsmethode. Berechnet fuer den durch die ROI festgelegten Bildausschnitt um die Niere
	 * den Kompassgradienten. Dies entspricht einem rechtungsabhaengigen Sobel-Filter.
	 * die Methode braucht dazu lediglich den Winkel (0-359°), die Bilddaten koennen
	 * aus Klassenvariablen uebernommen werden.
	 * @param degrees Der Winkel der Bilckrichtung.
	 * @return ShortProcessor das berechnete Kompassgradienten-Bild als ImageJ-ShortProcessor-Datentyp.
	 */
	private ShortProcessor getRotateSobel(double degrees) {
		int x, y, nx, ny;
		int deg = (int) (degrees + 22.5) % 360;

		// -----------------------------------------------------------------
		if (deg >= 0 && deg < 45) {
			if (sobelKernel != null && lastDeg >= 0 && lastDeg < 45) {
				lastDeg = deg;
				return aktSobelIp;
			} else {
				sobelKernel = new ShortProcessor(5, 5, false);
				nx = 4;
				for (x = 0; x < 5; x++) {
					for (y = 0; y < 5; y++) {
						sobelKernel.putPixel(nx, y, vertikalKernel.getPixel(x, y));
					} 
					nx--;
				} 
				lastDeg = deg;
				aktSobelIp = getSobel(aktIp, sobelKernel);
				return aktSobelIp;
			} 
		} 

		// -----------------------------------------------------------------
		if (deg >= 45 && deg < 90) {
			if (sobelKernel != null && lastDeg >= 45 && lastDeg < 90) {
				lastDeg = deg;
				return aktSobelIp;
			} else {
				sobelKernel = new ShortProcessor(5, 5, false);
				nx = 4;
				for (x = 0; x < 5; x++) {
					ny = 4;
					for (y = 0; y < 5; y++) {
						sobelKernel.putPixel(nx, ny, diagonalKernel.getPixel(x, y));
						ny--;
					} 
					nx--;
				} 
				lastDeg = deg;
				aktSobelIp = getSobel(aktIp, sobelKernel);
				return aktSobelIp;
			} 
		} 

		// -----------------------------------------------------------------
		if (deg >= 90 && deg < 135) {
			if (sobelKernel != null && lastDeg >= 90 && lastDeg < 135) {
				lastDeg = deg;
				return aktSobelIp;
			} else {
				sobelKernel = new ShortProcessor(5, 5, false);
				nx = 4;
				for (x = 0; x < 5; x++) {
					ny = 0;
					for (y = 0; y < 5; y++) {
						sobelKernel.putPixel(ny, nx, vertikalKernel.getPixel(x, y));
						ny++;
					} 
					nx--;
				} 
				lastDeg = deg;
				aktSobelIp = getSobel(aktIp, sobelKernel);
				return aktSobelIp;
			} 
		} 

		// -----------------------------------------------------------------
		if (deg >= 135 && deg < 180) {
			if (sobelKernel != null && lastDeg >= 135 && lastDeg < 180) {
				lastDeg = deg;
				return aktSobelIp;
			} else {
				sobelKernel = new ShortProcessor(5, 5, false);
				nx = 4;
				for (x = 0; x < 5; x++) {
					ny = 0;
					for (y = 0; y < 5; y++) {
						sobelKernel.putPixel(ny, nx, diagonalKernel.getPixel(x, y));
						ny++;
					} 
					nx--;
				} 
				lastDeg = deg;
				aktSobelIp = getSobel(aktIp, sobelKernel);
				return aktSobelIp;
			} 
		} 

		// -----------------------------------------------------------------
		if (deg >= 180 && deg < 225) {
			if (sobelKernel != null && lastDeg >= 180 && lastDeg < 225) {
				lastDeg = deg;
				return aktSobelIp;
			} else {
				sobelKernel = vertikalKernel;
				lastDeg = deg;
				aktSobelIp = getSobel(aktIp, sobelKernel);
				return aktSobelIp;
			} 
		} 

		// -----------------------------------------------------------------
		if (deg >= 225 && deg < 270) {
			if (sobelKernel != null && lastDeg >= 225 && lastDeg < 270) {
				lastDeg = deg;
				return aktSobelIp;
			} else {
				sobelKernel = diagonalKernel;
				lastDeg = deg;
				aktSobelIp = getSobel(aktIp, sobelKernel);
				return aktSobelIp;
			} 
		} 

		// -----------------------------------------------------------------
		if (deg >= 270 && deg < 315) {
			if (sobelKernel != null && lastDeg >= 270 && lastDeg < 315) {
				lastDeg = deg;
				return aktSobelIp;
			} else {
				sobelKernel = new ShortProcessor(5, 5, false);
				for (x = 0; x < 5; x++) {
					ny = 4;
					for (y = 0; y < 5; y++) {
						sobelKernel.putPixel(ny, x, vertikalKernel.getPixel(x, y));
						ny--;
					} 
				} 
				lastDeg = deg;
				aktSobelIp = getSobel(aktIp, sobelKernel);
				return aktSobelIp;
			} 
		} 

		// -----------------------------------------------------------------
		if (deg >= 315 && deg < 360) {
			if (sobelKernel != null && lastDeg >= 315 && lastDeg < 360) {
				lastDeg = deg;
				return aktSobelIp;
			} else {
				sobelKernel = new ShortProcessor(5, 5, false);
				nx = 4;
				for (x = 0; x < 5; x++) {
					ny = 4;
					for (y = 0; y < 5; y++) {
						sobelKernel.putPixel(ny, x, diagonalKernel.getPixel(x, y));
						ny--;
					} 
				} 
				lastDeg = deg;
				aktSobelIp = getSobel(aktIp, sobelKernel);
				return aktSobelIp;
			} 
		} 
		return aktSobelIp;
	} 


	/**
	 * Hilfsmethode. Wird zur Berechnung des Kompassgradienten gebraucht.
	 * Wendet auf eine uebergebenes Bild einen ebenfalls uebergebenen Kernel an.
	 * @param org Das Eingabebild als ImageJ-ShortProcessor-Datentyp
	 * @param kernel Der Kernel der Filterung, ebenfalls ImageJ-ShortProcessor-Datentyp
	 * @return ShortProcessor das gefilterte Bild.
	 */
	private ShortProcessor getSobel(ShortProcessor org, ShortProcessor kernel) {
		int			x, y, i, j, val;
		short[] res = new short[org.getWidth() * org.getHeight()];

		// result.setColorModel(org.getColorModel());
		for (x = 0; x < org.getWidth() - 5; x++) {
			for (y = 0; y < org.getHeight() - 5; y++) {
				val = 0;
				for (i = 0; i < 5; i++) {
					for (j = 0; j < 5; j++) {
						val = val + (int) (aktPixel[((y + j) * org.getWidth()) + x + i] /* org.getPixel(x+i,y+i) */ * kernel.getPixel(i, j));
					} 
				} 
				if (val < 0) {
					val = 0;	// Math.abs(val);

					// val = val;
				} 
				res[((y + 2) * org.getWidth()) + x + 2] = (short) val;
			} 
		} 
		String	k = "\n";

		for (y = 0; y < 5; y++) {
			for (x = 0; x < 5; x++) {
				val = kernel.getPixel(x, y);

				// if (val == 0) k = k + "0";
				// if (val == 0) k = k + "*";
				if (val >= 0) {
					k = k + val + " ";
				} 
				if (val == -1) {
					k = k + "# ";
				} 
			} 
			k = k + "\n";
		} 

		// IJ.write( "\nKernel:\n " + k);

		ShortProcessor	result = new ShortProcessor(org.getWidth(), org.getHeight(), res, org.getColorModel(), false);

		// ImagePlus newPlus = new ImagePlus( "Tester", result);
		// newPlus.show();
		return result;
	} 

}


/**
 * Hilfsklasse, wird zur Segmentierung benoetigt. Sie untersucht die Intensitaetsprofile
 * entlang zweier Scanlines und berechnet anhand vorgegebener Intensitaesuebergaenge
 * den bestmoeglichen Konturpunkt.
 * @see jm.segmentation.KidneyExtraction
 * @version  1.0, 14/02/2000
 * @author   Jens Martin
 */
class DoubleLine {
	double[]	line1;
	double[]	line2;


	/**
	 * Der Konstruktor.
	 * @param    l1    double-Array, Intensitaetsprofil der Scanline 1
	 * @param    l2    double-Array, Intensitaetsprofil der Scanline 2
	 */
	public DoubleLine(double[] l1, double[] l2) {
		copyFrom(l1, l2);
	}


	/**
	 * Hilfsmethode. Kopiert die Werte der uebergebenen Intensitaetsprofile in Klassenvaraiblen.
	 * @param    l1    double-Array, Intensitaetsprofil der Scanline 1
	 * @param    l2    double-Array, Intensitaetsprofil der Scanline 2
	 */
	public void copyFrom(double[] l1, double[] l2) {
		line1 = new double[l1.length];
		line2 = new double[l2.length];
		for (int i = 0; i < l1.length; i++) {
			line1[i] = l1[i];
		} 
		for (int i = 0; i < l2.length; i++) {
			line2[i] = l2[i];
		} 
	} 


	/**
	 * Untersucht die beiden Intensitaetsprofile und liefert den gefundenen Konturpunkt
	 * zurueck. Die Kriterien zur Detektion des Punktes sind in dieser Methode
	 * fest kodiert und wurden durch Tests ermittelt.
	 * @param    idxStart  Arrayindex, ab dem das Intensitaetsprofil der Scanlines untersucht werden soll.
	 * @param    idxEnd    Arrayindex, bis zu dem das Intensitaetsprofil der Scanlines untersucht werden soll.
	 * @return   int       Position im Array, die dem gefundenen Kontupunkt entspricht.
	 */
	public int findBestPoint(int idxStart, int idxEnd) {
		int			idxGoodMax = -1;
		int			idx = idxStart;
		double	line1Max = 0.0;

		while ((idx != idxEnd) && (idx < line1.length)) {
			if (line1[idx] > line1Max) {
				line1Max = line1[idx];
			} 
			idx = next(idxStart, idxEnd, idx);
		} 

		idx = idxStart;
		double	absMax = Double.MIN_VALUE;

		while ((idx != idxEnd) && (idx < line1.length) && (idx < line2.length)) {
			double	val = (line2[idx] / 2.0) + (line1Max - line1[idx]);

			// val = val * (double)idx / (double)line1.length;
			if (val > absMax) {
				absMax = val;
				idxGoodMax = idx;
			} 

			// if ( val > 100) break;
			idx = next(idxStart, idxEnd, idx);
		} 
		return idxGoodMax;
	} 


	/**
	 * Liest eine Intensitaetswert an der uebergebenen Position innerhalb des Intensitaetprofils aus.
	 * @param    idxPos    Arrayindex
	 * @return   double    Intensitaeswert
	 */
	public double getValueAt(int idxPos) {
		int			idx = 0;
		int			idxStart = 0;
		int			idxEnd = line1.length - 1;
		double	line1Max = 0.0;

		while (idx != idxEnd) {
			if (line1[idx] > line1Max) {
				line1Max = line1[idx];
			} 
			idx = next(idxStart, idxEnd, idx);
		} 

		double	val = (line2[idxPos] / 2.0) + (line1Max - line1[idx]);

		// IJ.write("Val: " + val + "/" + ( line2[idxPos] / 2.0 ) + "/" + ( line1Max - line1[idx] ));
		return val;
	} 


	/**
	 * Hilfsmethode. Glaettet die uebergebenen Intensitaets-Kennlinien durch Mittelwert-Bildung aus den
	 * Nachbarpixeln.
	 * @param    range     legt fest, wieviele Nachbarwerte mit eingezogen werden sollen
	 */
	public void smoothLines(int range) {
		double	floatingMid;
		int			i, j, pos;

		for (i = 0; i < line1.length; i++) {
			floatingMid = 0.0;
			for (j = 0; j < range; j++) {
				pos = i - (range / 2) + j;
				if (pos < 0) {
					pos = pos + range;
				} 
				if (pos >= line1.length) {
					pos = pos - range;
				} 
				floatingMid = floatingMid + line1[pos];

			} 
			floatingMid = floatingMid / (double) range;

			// IJ.write("FM1: " + floatingMid);
			line1[i] = floatingMid;
		} 

		for (i = 0; i < line2.length; i++) {
			floatingMid = 0.0;
			for (j = 0; j < range; j++) {
				pos = i - (range / 2) + j;
				if (pos < 0) {
					pos = pos + range;
				} 
				if (pos >= line2.length) {
					pos = pos - range;
				} 
				floatingMid = floatingMid + line2[pos];

			} 
			floatingMid = floatingMid / (double) range;

			// IJ.write("FM2: " + floatingMid);
			line2[i] = floatingMid;
		} 

	} 


	/**
	 * Hilfsmethode. Liefert den naechsten Index entsprechned der Durchlaufrichtung
	 * durch ein Array.
	 * @param    idxStart     Startposition
	 * @param    idxEnd       Endposition (muss nicht groesser als die Startposition sein !!!)
	 * @param    idx          aktueller Index
	 * @return   int naechste Position
	 */
	public int next(int idxStart, int idxEnd, int idx) {
		if (idxStart < idxEnd) {
			return ++idx;
		} else if (idxStart > idxEnd) {
			return --idx;
		} 
		return idx;
	} 

}


/*--- formatting done in "My Own Convention" style on 08-04-2000 ---*/

