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

/*
 * The package jm.sit is part of the Renal Function Project
 * for analysis of dynamic contrast medium evalutions 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.sit;

import java.awt.*;
import java.lang.Math.*;
import java.awt.image.*;

import ij.*;
import ij.io.*;
import ij.gui.*;
import ij.process.*;

import jm.kidney.*;
import jm.cepstrum.*;
import jm.sit.*;
import jm.extension.*;
import jm.util.*;


/**
 * Klasse zur Berechnung der Signal-Intensitaets-Zeit-Kurven.
 * Realisiert die Pixelklassifizierung fuer Neirenrinde und Nierenmark.
 * Wird von der Klasse AnalysisDialog aufgerufen und benoetigt die Klasse
 * SITResults zur Darstellung der berechneten Ergenisse in einm Dialogfenster
 * @see jm.kidney.AnalysisDialog
 * @see jm.sit.SITResults
 * @version  4.0, 14/02/2000
 * @author   Jens Martin
 */
public class SITAnalysis extends Thread {
	private double					grayAdd = 30.0;
	public double						devScaleFactor = 0.001;

	private double					adjustmentMin = 0.5;
	private double					adjustmentMax = 1.5;

	private int							maxImage = 0;
	private int							maxTime = 0;
	public AnalysisDialog		analysisDialog = null;

	public SITResults				sitResults = null;
	public ImagePlus				rightRefImage = null;
	public ImagePlus				leftRefImage = null;

	private Color						myYellow = new Color(243, 239, 118);
	private Color						myRed = new Color(227, 75, 75);

	private ImagePlus				imagePlus = null;
	private ImageStack			imageStack = null;
	private ExImageCanvas		imageCanvas = null;

	// private ImageStack newStack = null;
	private ImageProcessor	imageProcessor = null;
	private int							width = 0;
	private int							height = 0;
	private Roi							rRoi = null;
	private Roi							lRoi = null;
	private Rectangle				rRect = null;
	private Rectangle				lRect = null;
	private double[][]			stdDeviation = null;
	private double[][]			averageGray = null;
	private double[][]			stdDeviationMedulla = null;
	private double[][]			averageGrayMedulla = null;
	private double[][]			stdDeviationCortex = null;
	private double[][]			averageGrayCortex = null;
	private double[][]			secondDerivation = null;

	public int[]						pixelCountKidney = null;
	public int[]						pixelCountMedulla = null;
	public int[]						pixelCountCortex = null;

	public int[]						surfaceKidney = null;
	public int[]						surfaceMedulla = null;
	public int[]						surfaceCortex = null;

	private double[]				pixStdDeviation = null;
	private double[]				pixAverageGray = null;
	private double[]				pixGrayAscent = null;
	private short[]					pixClassification = null;
	private short[]					pixDistancePercent = null;

	private short[]					pix = null;

	private int[][]					centerGravX = null;
	private int[][]					centerGravY = null;

	private int[]						maxGray = null;
	private int[]						imageNrTreshhold = null;
	private double[]				maxDeviation = null;
	private double[]				maxAverageGray = null;
	private double[]				maxGrayAscent = null;
	public double[]					acquisitionTimes = null;

	public ImagePlus				ipLeftKidney = null;
	public ImagePlus				ipRightKidney = null;
	private ColorProcessor	cpClassifyHelper = null;
	public boolean[]				autoClassification = null;
	public boolean[]				noRois = null;
	private boolean[][]			enabledForSIT = null;

	public boolean					stopRunning = false;
	public boolean					finished = false;
	public boolean					error = false;

	private boolean					verbose = true;
	private int							startImage = 1;

	public double[]					classTresholdPelvis = null;
	public double[]					classTresholdMedulla = null;
	public double[]					classTresholdCortex = null;

	private Roi[]						leftRoi = null;
	private Roi[]						rightRoi = null;
	private Roi[][]					leftRoiCortex = null;
	private Roi[][]					rightRoiCortex = null;
	private Roi[][]					leftRoiMedulla = null;
	private Roi[][]					rightRoiMedulla = null;


	/**
	 * Der Konstruktor. Er bekommt eine Referenz auf das Hauptfenster der renalen
	 * Funktionsanalyse uebergeben, um Zugriff auf dortige Funktionen und Variablen
	 * zu haben.
	 * @param dlg Die Referenz auf das Plugin-Hauptfenster
	 */
	public SITAnalysis(AnalysisDialog dlg)	// , ImagePlus img, int[] imageTimes, Roi leftRoi, Roi rightRoi)
	 {
		try {
			analysisDialog = dlg;
			imagePlus = (ImagePlus) analysisDialog.imagePlus;
			imageStack = imagePlus.getStack();
			maxImage = imagePlus.getStackSize();
			imageCanvas = analysisDialog.imageCanvas;
			startImage = analysisDialog.startImage;

			height = imagePlus.getHeight();
			width = imagePlus.getWidth();

			leftRoi = imageCanvas.getRoi(JMConst.LEFT);
			rightRoi = imageCanvas.getRoi(JMConst.RIGHT);
			leftRoiCortex = imageCanvas.getRoi(JMConst.LEFT, JMConst.CORTEX);
			rightRoiCortex = imageCanvas.getRoi(JMConst.RIGHT, JMConst.CORTEX);
			leftRoiMedulla = imageCanvas.getRoi(JMConst.LEFT, JMConst.MEDULLA);
			rightRoiMedulla = imageCanvas.getRoi(JMConst.RIGHT, JMConst.MEDULLA);
			enabledForSIT = imageCanvas.getEnabledForSIT();

			autoClassification = new boolean[2];
			noRois = new boolean[2];

			// autoClassification[JMConst.LEFT] = true;
			// autoClassification[JMConst.RIGHT] = true;

			for (int half = JMConst.LEFT; half <= JMConst.RIGHT; half++) {
				autoClassification[half] = true;
				noRois[half] = false;
				Roi[]		roi = leftRoi;
				Roi[][] cortexRoi = leftRoiCortex;
				Roi[][] medullaRoi = leftRoiMedulla;

				if (half == JMConst.RIGHT) {
					roi = rightRoi;
					cortexRoi = rightRoiCortex;
					medullaRoi = rightRoiMedulla;
				}

				int lastX = -1;
				int lastY = -1;
				int x = -1;
				int y = -1;
				int counter = 0;

				for (int i = startImage - 1; i < maxImage; i++) {
					if (roi != null && roi[i] != null && enabledForSIT[half][i]) {
						counter++;
						x = roi[i].getBoundingRect().x;
						y = roi[i].getBoundingRect().y;
						if (lastX < 0) {
							lastX = x;
						}
						if (lastY < 0) {
							lastY = y;
						}
						if (lastX != x || lastY != y || (cortexRoi != null && cortexRoi[i][0] != null) || (medullaRoi != null && medullaRoi[i][0] != null)) {
							autoClassification[half] = false;
							i = maxImage;
						}
					}
					lastX = x;
					lastY = y;
				}
				if (autoClassification[half] && counter == 0) {
					autoClassification[half] = false;

					// Check if there are any Rois
					for (int i = startImage - 1; i < maxImage; i++) {
						if ((cortexRoi != null && cortexRoi[i][0] != null) || (medullaRoi != null && medullaRoi[i][0] != null)) {
							counter++;
						}
					}
					if (counter == 0) {
						noRois[half] = true;
					}
				}
			}

			if (verbose) {
				IJ.write("Checking both ROIs...");
			}
			if (autoClassification[JMConst.LEFT]) {
				int i = 0;

				while (lRoi == null) {
					lRoi = leftRoi[i++];
				}
				lRect = lRoi.getBoundingRect();
				if (verbose) {
					IJ.write("Left ROI x:" + lRect.x + " y: " + lRect.y + " w: " + lRect.width + " h: " + lRect.height);
				}
			} else {
				lRect = new Rectangle(Integer.MAX_VALUE, Integer.MAX_VALUE, Integer.MIN_VALUE, Integer.MIN_VALUE);
				Rectangle rect;

				for (int i = startImage - 1; i < maxImage; i++) {
					if (enabledForSIT[JMConst.LEFT][i]) {
						if (leftRoi != null && leftRoi[i] != null) {
							rect = leftRoi[i].getBoundingRect();
							if (rect.x < lRect.x) {
								lRect.x = rect.x;
							}
							if (rect.y < lRect.y) {
								lRect.y = rect.y;
							}
							if (rect.width > lRect.width) {
								lRect.width = rect.width;
							}
							if (rect.height > lRect.height) {
								lRect.height = rect.height;
							}
						}
						int j = 0;

						while (j < imageCanvas.MAXROIS && leftRoiCortex != null && leftRoiCortex[i][j] != null) {
							rect = leftRoiCortex[i][j].getBoundingRect();
							if (rect.x < lRect.x) {
								lRect.x = rect.x;
							}
							if (rect.y < lRect.y) {
								lRect.y = rect.y;
							}
							if (rect.width > lRect.width) {
								lRect.width = rect.width;
							}
							if (rect.height > lRect.height) {
								lRect.height = rect.height;
							}
							j++;
						}
						j = 0;
						while (j < imageCanvas.MAXROIS && leftRoiMedulla != null && leftRoiMedulla[i][j] != null) {
							rect = leftRoiMedulla[i][j].getBoundingRect();
							if (rect.x < lRect.x) {
								lRect.x = rect.x;
							}
							if (rect.y < lRect.y) {
								lRect.y = rect.y;
							}
							if (rect.width > lRect.width) {
								lRect.width = rect.width;
							}
							if (rect.height > lRect.height) {
								lRect.height = rect.height;
							}
							j++;
						}
					}
				}
				if (lRect.width < 0) {
					lRect = null;
				}
			}

			if (autoClassification[JMConst.RIGHT]) {
				int i = 0;

				while (rRoi == null) {
					rRoi = rightRoi[i++];
				}
				rRect = rRoi.getBoundingRect();
				IJ.write("Right ROI x:" + rRect.x + " y: " + rRect.y + " w: " + rRect.width + " h: " + rRect.height);
			} else {
				rRect = new Rectangle(Integer.MAX_VALUE, Integer.MAX_VALUE, Integer.MIN_VALUE, Integer.MIN_VALUE);
				Rectangle rect;

				for (int i = startImage - 1; i < maxImage; i++) {
					if (enabledForSIT[JMConst.RIGHT][i]) {
						if (rightRoi != null && rightRoi[i] != null) {
							rect = rightRoi[i].getBoundingRect();
							if (rect.x < rRect.x) {
								rRect.x = rect.x;
							}
							if (rect.y < rRect.y) {
								rRect.y = rect.y;
							}
							if (rect.width > rRect.width) {
								rRect.width = rect.width;
							}
							if (rect.height > rRect.height) {
								rRect.height = rect.height;
							}
						}
						int j = 0;

						while (j < imageCanvas.MAXROIS && rightRoiCortex != null && rightRoiCortex[i][j] != null) {
							rect = rightRoiCortex[i][j].getBoundingRect();
							if (rect.x < rRect.x) {
								rRect.x = rect.x;
							}
							if (rect.y < rRect.y) {
								rRect.y = rect.y;
							}
							if (rect.width > rRect.width) {
								rRect.width = rect.width;
							}
							if (rect.height > rRect.height) {
								rRect.height = rect.height;
							}
							j++;
						}
						j = 0;
						while (j < imageCanvas.MAXROIS && rightRoiMedulla != null && rightRoiMedulla[i][j] != null) {
							rect = rightRoiMedulla[i][j].getBoundingRect();
							if (rect.x < rRect.x) {
								rRect.x = rect.x;
							}
							if (rect.y < rRect.y) {
								rRect.y = rect.y;
							}
							if (rect.width > rRect.width) {
								rRect.width = rect.width;
							}
							if (rect.height > rRect.height) {
								rRect.height = rect.height;
							}
							j++;
						}
					}
				}
				if (rRect.width < 0) {
					rRect = null;
				}
			}

			int[] imageTimes = analysisDialog.acquisitionTime;

			if (imageTimes != null) {
				int arrayLen = imageTimes.length;

				maxTime = imageTimes[arrayLen - 1] / 1000;
				if (maxTime == 0) {
					maxTime = arrayLen * 2;
				}
				acquisitionTimes = new double[arrayLen];
				for (int i = 0; i < arrayLen; i++) {
					acquisitionTimes[i] = (double) imageTimes[i] / 1000.0;
				}
			}

			// initImages(false);
		} catch (Exception e) {
			IJ.write("Can't initialise Result Window");
			error = true;
			e.printStackTrace();
		}
	}


	/**
	 * Hilfsmethode. Legt fest, bei welchem Bild der Sequenz mit den Berechnungen
	 * begonnen werden soll. Wird vor dem Start der Berechnungen mit dem im
	 * Voreinstllungs-Dailog angegebenen Wert aufgerufen.
	 * @param n die Bildnummer innerhalb der Sequenz.
	 */
	public void setStartImage(int n) {
		if (n < maxImage && n >= 1) {
			startImage = n;
		}
	}


	/**
	 * Initialisiert der beiden eingefaerbten Klassifizierungsbilder fuer linke und rechte Niere
	 * sowie zwei Hilfsbilder, die zur Klassifizierung benoetigt werden. Dieser werden mehr im Sinne
	 * eines Arrays gebraucht.
	 * @param showImages true, wenn die Klassifizierungsbilder und Hilfsbilder zur Laufzeit angezeigt werden sollen.
	 */
	private void initImages(boolean showImages) {

		// imageNrTreshhold[JMConst.LEFT]
		short[] pixLeft = (short[]) imageStack.getPixels(imageNrTreshhold[JMConst.LEFT]);
		short[] pixRight = (short[]) imageStack.getPixels(imageNrTreshhold[JMConst.RIGHT]);

		pix = new short[width * height];
		for (int y = 0; y < height; y++) {
			for (int x = 0; x < width / 2; x++) {
				pix[(y * width) + x] = pixLeft[(y * width) + x];
			}
			for (int x = width / 2; x < width; x++) {
				pix[(y * width) + x] = pixRight[(y * width) + x];
			}
		}

		// pix = (short[])imageStack.getPixels(startImage);
		cpClassifyHelper = new ColorProcessor(width, height);

		for (int x = 0; x < width; x++) {
			for (int y = 0; y < height; y++) {
				double	gray = (double) pix[(y * width) + x];

				gray = gray + grayAdd;
				int r = (int) gray;
				int g = (int) gray;
				int b = (int) gray;

				if (r > 255) {
					r = 255;
				}
				if (g > 255) {
					g = 255;
				}
				if (b > 255) {
					b = 255;
				}
				int val = 0xff000000 | ((r & 0xff) << 16) | ((g & 0xff) << 8) | b & 0xff;

				cpClassifyHelper.putPixel(x, y, val);
			}
		}
		Rectangle Rect = null;

		for (int half = JMConst.LEFT; half <= JMConst.RIGHT; half++) {
			if (half == JMConst.LEFT) {
				Rect = lRect;
			}
			if (half == JMConst.RIGHT) {
				Rect = rRect;
			}
			ColorProcessor	cp = null;

			if (Rect != null) {
				cp = new ColorProcessor(Rect.width, Rect.height);
				for (int x = 0; x < Rect.width; x++) {
					for (int y = 0; y < Rect.height; y++) {
						double	gray = (double) pix[((y + Rect.y) * width) + (x + Rect.x)];

						gray = gray + grayAdd;
						int r = (int) gray;
						int g = (int) gray;
						int b = (int) gray;

						if (r > 255) {
							r = 255;
						}
						if (g > 255) {
							g = 255;
						}
						if (b > 255) {
							b = 255;

						}
						int val = 0xff000000 | ((r & 0xff) << 16) | ((g & 0xff) << 8) | b & 0xff;

						cp.putPixel(x, y, val);
					}
				}
			} else {

				// generate a black Dummy-Image
				cp = new ColorProcessor(20, 20);
				for (int x = 0; x < 20; x++) {
					for (int y = 0; y < 20; y++) {
						int val = 0x00000000;

						cp.putPixel(x, y, val);
					}
				}
			}
			if (half == JMConst.LEFT) {
				ipLeftKidney = new ImagePlus("right Kidney", cp);
				if (showImages) {
					ipLeftKidney.show();
				}
			} else {
				ipRightKidney = new ImagePlus("left Kidney", cp);
				if (showImages) {
					ipRightKidney.show();
				}
			}
		}		// end for
 	 }


	/**
	 * Haupt-Berechnungsmethode der Klasse. Wird aufgerufen, wenn der Thread gestartet
	 * wird. Implementiert die abstrakte run-Method der Klasse Thread.
	 */
	public void run() {
		int i = 0;
		int x = 0;
		int y = 0;

		stopRunning = false;

		// ColorModel cm = imageStack.getColorModel();
		// if ( !(lRect != null || rRect != null) )
		if (noRois[JMConst.LEFT] && noRois[JMConst.RIGHT]) {
			this.finished = false;	// set finish flag
			this.stopRunning = true;
			return;
		}

		try {
			int[]		mask = null;													// mask of the roi
			short[] pix = null;														// pixel array of actual slice

			IJ.showProgress(0.001);
			int			maxProgress = 0;
			int			counter = 0;
			double	sum;
			int			grayVal, xSum, ySum, graySum;

			// arrays to store series properties in
			stdDeviation = new double[2][maxImage /* +1 */];
			averageGray = new double[2][maxImage /* +1 */];
			stdDeviationMedulla = new double[2][maxImage /* +1 */];
			averageGrayMedulla = new double[2][maxImage /* +1 */];
			stdDeviationCortex = new double[2][maxImage /* +1 */];
			averageGrayCortex = new double[2][maxImage /* +1 */];
			secondDerivation = new double[2][maxImage /* +1 */];
			centerGravX = new int[2][maxImage /* +1 */];
			centerGravY = new int[2][maxImage /* +1 */];
			maxGray = new int[2];
			imageNrTreshhold = new int[2];
			maxDeviation = new double[2];
			maxAverageGray = new double[2];
			maxGrayAscent = new double[2];

			pixelCountKidney = new int[2];
			pixelCountMedulla = new int[2];
			pixelCountCortex = new int[2];

			surfaceKidney = new int[2];
			surfaceMedulla = new int[2];
			surfaceCortex = new int[2];

			classTresholdPelvis = new double[2];
			classTresholdMedulla = new double[2];
			classTresholdCortex = new double[2];

			// arrays to store pixel properties in
			pixStdDeviation = new double[width * height];
			pixAverageGray = new double[width * height];
			pixGrayAscent = new double[width * height];
			pixClassification = new short[width * height];
			pixDistancePercent = new short[width * height];

			if (lRect != null) {
				maxProgress++;
			}
			if (rRect != null) {
				maxProgress++;
			}
			maxProgress = maxProgress * 2 * maxImage;

			// for-next from 0 to 1 for LEFT and RIGHT kidney
			int half;

			for (half = 0; half <= 1; half++) {
				pixelCountKidney[half] = 0;

				classTresholdPelvis[half] = 1.0;
				classTresholdMedulla[half] = 1.0;
				classTresholdCortex[half] = 1.0;

				Rectangle rect = lRect;
				Roi				roi = lRoi;

				if (half == JMConst.RIGHT)									// half == 1
				 {
					rect = rRect;
					roi = rRoi;
				}
				i = startImage;
				while (i <= maxImage) {
					double[]	roiProperties = new double[6];

					for (int j = 0; j < 6; j++) {
						roiProperties[j] = 0.0;
					}
					if (half == JMConst.LEFT && enabledForSIT[JMConst.LEFT][i - 1] && leftRoi != null) {
						roiProperties = getRoiProperties(i, leftRoi[i - 1]);
					}
					if (half == JMConst.RIGHT && enabledForSIT[JMConst.RIGHT][i - 1] && rightRoi != null) {
						roiProperties = getRoiProperties(i, rightRoi[i - 1]);

						// thread: short break for other processes
					}
					try {
						this.sleep(5);
					} catch (InterruptedException ie) {
						IJ.write("Interrupted in SIT-calcualtion");
					}

					averageGray[half][i - 1] = roiProperties[0];
					stdDeviation[half][i - 1] = roiProperties[1];
					centerGravX[half][i - 1] = (int) roiProperties[2];
					centerGravY[half][i - 1] = (int) roiProperties[3];
					maxGray[half] = (int) roiProperties[4];
					pixelCountKidney[half] = (int) roiProperties[5];

					if (verbose) {
						IJ.write("Slice: " + i + " Half:" + half + " - Average Gray: " + averageGray[half][i - 1] + " StandardDeviation: " + stdDeviation[half][i - 1]);

						// increament i for while
					}
					i++;
				}																						// while end
 			 }																						// for end

			// calculate the turn over point
			imageNrTreshhold[JMConst.LEFT] = 1;
			imageNrTreshhold[JMConst.RIGHT] = 1;

			if (autoClassification[JMConst.LEFT]) {
				imageNrTreshhold[JMConst.LEFT] = getMaxArrayPos(derivation(floatingAverage(floatingAverage(averageGray[JMConst.LEFT], 3), 3))) + 3;		// + 1;
				if (verbose) {
					IJ.write("Turn Over-Point for Right Kidney is Image: " + imageNrTreshhold[JMConst.LEFT]);
				}
			}
			if (autoClassification[JMConst.RIGHT]) {
				imageNrTreshhold[JMConst.RIGHT] = getMaxArrayPos(derivation(floatingAverage(floatingAverage(averageGray[JMConst.RIGHT], 3), 3))) + 3;		// + 1;
				if (verbose) {
					IJ.write("Turn Over-Point for Left Kidney is Image: " + imageNrTreshhold[JMConst.RIGHT]);
				}
			}


			if (imageNrTreshhold[JMConst.LEFT] > maxImage) {
				imageNrTreshhold[JMConst.LEFT] = maxImage;
			}
			if (imageNrTreshhold[JMConst.RIGHT] > maxImage) {
				imageNrTreshhold[JMConst.RIGHT] = maxImage;

			}
			initImages(false);

			// secondDerivation[JMConst.LEFT] = derivation(averageGray[JMConst.LEFT]);
			// secondDerivation[JMConst.RIGHT] = derivation(averageGray[JMConst.RIGHT]);
			// imagePlus.setSlice(imageNrTreshhold[JMConst.LEFT]);
			// imagePlus.updateAndDraw();
			if ((autoClassification[JMConst.LEFT]) || (autoClassification[JMConst.RIGHT])) {
				for (x = 0; x < width; x++) {
					for (y = 0; y < height; y++) {
						pixDistancePercent[(y * width) + x] = 100;
					}
				}

				for (half = 0; half <= 1; half++) {
					maxAverageGray[half] = 0.0;
					maxDeviation[half] = 0.0;
					maxGrayAscent[half] = 0.0;

					Rectangle rect = lRect;
					Roi				roi = lRoi;

					if (half == JMConst.RIGHT) {
						rect = rRect;
						roi = rRoi;
					}
					if (rect != null && autoClassification[half]) {
						mask = roi.getMask();
						if (mask != null) {

							// calc average gray for each pixel of the rois over all images
							for (x = 0; x < rect.width; x++) {
								for (y = 0; y < rect.height; y++) {
									if (mask[(y * (rect.width)) + x] != -1) {
										counter = 0;
										sum = 0.0;
										for (i = startImage; i <= maxImage; i++) {
											if (enabledForSIT[half][i - 1]) {
												pix = (short[]) imageStack.getPixels(i);
												sum = sum + (double) pix[((y + rect.y) * width) + (x + rect.x)];
												counter++;
											}
										}
										sum = sum / (double) counter;		// maxImage;
										pixAverageGray[((y + rect.y) * width) + (x + rect.x)] = sum;
										if (maxAverageGray[half] < sum) {
											maxAverageGray[half] = sum;
										}
									}
								}
							}

							// calc standard deviation for each Pixel of the left Roi over all Images
							for (x = 0; x < rect.width; x++) {
								for (y = 0; y < rect.height; y++) {
									if (mask[(y * (rect.width)) + x] != -1) {
										counter = 0;
										sum = 0.0;
										for (i = startImage; i <= maxImage; i++) {
											if (enabledForSIT[half][i - 1]) {
												pix = (short[]) imageStack.getPixels(i);
												double	d = ((double) pix[((y + rect.y) * width) + (x + rect.x)]) - pixAverageGray[((y + rect.y) * width) + (x + rect.x)];

												sum = sum + (d * d);
												counter++;
											}
										}
										double	dev = Math.sqrt((1.0 / ((double) counter)) * sum);

										pixStdDeviation[((y + rect.y) * width) + (x + rect.x)] = dev;
										if (maxDeviation[half] < dev) {
											maxDeviation[half] = dev;
										}
									}
								}
							}

							// IJ.write("maxmax:" +maxAverageGray[half] + " / " +maxDeviation[half] );

							if (verbose) {
								IJ.write("Calculating Help-Parameters...");

								// calculate the distance of each pixel to the center of gravity in percent

							}
							for (x = 0; x < rect.width; x++) {
								for (y = 0; y < rect.height; y++) {
									if (mask[(y * (rect.width)) + x] != -1) {
										int xa = -1;
										int ya = -1;
										int pos = startImage - 1;

										while (xa < 0 || ya < 0) {
											xa = centerGravX[half][pos] - rect.x;
											ya = centerGravY[half][pos] - rect.y;
											pos++;
										}
										int			xb = x;
										int			yb = y;
										double	dist = calcDistance(xa, ya, xb, yb);
										int			dx, dy;
										int			vx, vy, px, py, fa, fy, fx, fxy;

										// /////
										if (xb > xa) {
											dx = 1;
											vx = yb - ya;
										} else {
											dx = -1;
											vx = ya - yb;
										}
										if (yb > ya) {
											dy = 1;
											vy = xa - xb;
										} else {
											dy = -1;
											vy = xb - xa;
										}

										fa = 0;
										px = xa;
										py = ya;

										while ((px != xb) || (py != yb)) {
											if (((py * (rect.width)) + px < mask.length) && (mask[(py * (rect.width)) + px] != -1)) {

												// todo:
												// Winkel mit den Vektoren / sin-Fkt bestimmten
												// entfernung in Prozent ausrechnen
												double	d = calcDistance(xa, ya, px, py);
												short		percent = 0;

												if (dist > 0.0) {
													percent = (short) ((d / dist) * 100.0);
												}
												short oldVal = pixDistancePercent[((py + rect.y) * width) + (px + rect.x)];

												if (percent < oldVal) {

													// ins Array zeichnen
													pixDistancePercent[((py + rect.y) * width) + (px + rect.x)] = (short) percent;
												}
											}

											// ip.putPixel(px,py,254);
											fx = fa + vx;
											fy = fa + vy;
											fxy = fx + fy - fa;
											if (Math.abs(fx) <= Math.abs(fy) && Math.abs(fx) <= Math.abs(fxy)) {
												px += dx;
												fa = fx;
											} else if (Math.abs(fy) < Math.abs(fx) && Math.abs(fy) < Math.abs(fxy)) {
												fa = fy;
												py += dy;
											} else {
												fa = fxy;
												px += dx;
												py += dy;
											}
										}
									}
								}
							}

							// IJ.write("Now Gray ascent");
							// calc the middle-gray-ascent for each pixel of the rois from start to trunoverpoint
							for (x = 0; x < rect.width; x++) {
								for (y = 0; y < rect.height; y++) {
									if (mask[(y * (rect.width)) + x] != -1) {
										counter = 0;
										short val1, val2;

										sum = 0.0;
										for (i = startImage + 1; i <= imageNrTreshhold[half]; i++) {
											if (enabledForSIT[half][i - 1]) {
												pix = (short[]) imageStack.getPixels(i - 1);
												val1 = pix[((y + rect.y) * width) + (x + rect.x)];
												pix = (short[]) imageStack.getPixels(i);
												val2 = pix[((y + rect.y) * width) + (x + rect.x)];
												if (val2 < val1) {
													val2 = val1;
												}
												sum = sum + (double) (val2 - val1);
												counter++;
											}
										}
										sum = sum / (double) counter;		// imageNrTreshhold[half];
										sum = sum * 10;
										pixGrayAscent[((y + rect.y) * width) + (x + rect.x)] = sum;
										if (maxGrayAscent[half] < sum) {
											maxGrayAscent[half] = sum;
										}
									}
								}
							}

							// IJ.write("MGA:" + maxGrayAscent[half]);
						}
					}
				}
			}

			// Testzwecke: Bilder anzeigen
			boolean			showImages = false;

			ColorModel	cm = imagePlus.getProcessor().getColorModel();

			pix = new short[width * height];
			for (i = 0; i < width * height; i++) {
				pix[i] = (short) pixAverageGray[i];
			}
			ShortProcessor	sp1 = new ShortProcessor(width, height, pix, cm, false);
			ImagePlus				ip1 = new ImagePlus("Test1", sp1);

			if (showImages) {
				ip1.show();

				// ///////////
			}
			for (i = 0; i < width * height; i++) {
				pix[i] = (short) pixStdDeviation[i];
			}
			ShortProcessor	sp2 = new ShortProcessor(width, height, pix, cm, false);
			ImagePlus				ip2 = new ImagePlus("Test2", sp2);

			if (showImages) {
				ip2.show();

				// ///////////
			}
			ShortProcessor	sp3 = new ShortProcessor(width, height, pixDistancePercent, cm, false);
			ImagePlus				ip3 = new ImagePlus("Test3", sp3);

			if (showImages) {
				ip3.show();

				// ///////////
			}
			for (i = 0; i < width * height; i++) {
				pix[i] = (short) pixGrayAscent[i];
			}
			ShortProcessor	sp4 = new ShortProcessor(width, height, pix, cm, false);
			ImagePlus				ip4 = new ImagePlus("Test4", sp4);

			if (showImages) {
				ip4.show();

				// /////

				// Zeige die Wendepunkte der beiden SIT-Kurven
				// IJ.error("!" + imageNrTreshhold[JMConst.LEFT]+"\n"+imageNrTreshhold[JMConst.RIGHT]);
			}
			int j;

			if (lRect != null) {
				pix = (short[]) imageStack.getPixels(imageNrTreshhold[JMConst.LEFT]);
				short[] lPix = new short[lRect.width * lRect.height];

				for (i = 0; i < lRect.width; i++) {
					for (j = 0; j < lRect.height; j++) {
						lPix[(j * lRect.width) + i] = (short) pix[((lRect.y + j) * width) + lRect.x + i];
					}
				}
				ShortProcessor	spTemp = new ShortProcessor(lRect.width, lRect.height, lPix, cm, false);

				leftRefImage = new ImagePlus("Temp", spTemp);

				// leftRefImage = ipTemp.getImage();
			} else {

				// create a dummy Image
				short[] lPix = new short[20 * 20];

				for (i = 0; i < 20 * 20; i++) {
					lPix[i] = 0;
				}
				ShortProcessor	spTemp = new ShortProcessor(20, 20, lPix, cm, false);

				leftRefImage = new ImagePlus("Temp", spTemp);
			}

			if (rRect != null) {
				pix = (short[]) imageStack.getPixels(imageNrTreshhold[JMConst.RIGHT]);
				short[] rPix = new short[rRect.width * rRect.height];

				for (i = 0; i < rRect.width; i++) {
					for (j = 0; j < rRect.height; j++) {
						rPix[(j * rRect.width) + i] = (short) pix[((rRect.y + j) * width) + rRect.x + i];

						// (short)imageStack.getProcessor().getPixel(rRect.x+i,rRect.y+j);
					}
				}
				ShortProcessor	spTemp = new ShortProcessor(rRect.width, rRect.height, rPix, cm, false);

				rightRefImage = new ImagePlus("Temp", spTemp);

				// rightRefImage = ipTemp.getImage();
			} else {

				// create a dummy Image
				short[] rPix = new short[20 * 20];

				for (i = 0; i < 20 * 20; i++) {
					rPix[i] = 0;
				}
				ShortProcessor	spTemp = new ShortProcessor(20, 20, rPix, cm, false);

				leftRefImage = new ImagePlus("Temp", spTemp);
			}

			if (verbose) {
				IJ.write("Wendepunkt links: " + imageNrTreshhold[JMConst.LEFT] + "\nWendepunkt rechts: " + imageNrTreshhold[JMConst.RIGHT]);


			}
			for (half = JMConst.LEFT; half <= JMConst.RIGHT; half++) {
				int otherHalf = Math.abs(half - 1);

				if (verbose) {
					IJ.write("Classifying Pixels for Half: " + half);
				}
				if (autoClassification[half]) {
					classifyPixels(half, classTresholdPelvis[half], classTresholdMedulla[half], classTresholdCortex[half]);

					// IJ.write("classify Pixels hat geklappt");
					calcPixelCurve(half, JMConst.MEDULLA);

					// IJ.write("calc Pixel Curve medulla hat geklappt");
					calcPixelCurve(half, JMConst.CORTEX);

					// IJ.write("calc Pixel Curve cortex hat geklappt");

					if (analysisDialog.interpolateNullValues) {
						averageGray[half] = interpolateNullValues(averageGray[half]);
						stdDeviation[half] = interpolateNullValues(stdDeviation[half]);

						// interpolateNullValues(stdDeviationCortex[half]);
						// interpolateNullValues(stdDeviationMedulla[half]);
					}
					if (startImage > 1) {
						averageGray[half] = fillMissingValues(averageGray[half]);

						// averageGrayCortex[half] = fillMissingValues(averageGrayCortex[half]);
						// averageGrayMedulla[half] = fillMissingValues(averageGrayMedulla[half]);
						stdDeviation[half] = fillMissingValues(stdDeviation[half]);

						// stdDeviationCortex[half] = fillMissingValues(stdDeviationCortex[half]);
						// stdDeviationMedulla[half] = fillMissingValues(stdDeviationMedulla[half]);
					}
					initResultWindow(half);

				} else {
					Roi[][] roiMedulla = leftRoiMedulla;
					Roi[][] roiCortex = leftRoiCortex;

					if (half == JMConst.RIGHT) {
						roiMedulla = rightRoiMedulla;
						roiCortex = rightRoiCortex;
					}
					pixelCountCortex[half] = 0;
					pixelCountMedulla[half] = 0;
					int pixelCountCor = 0;
					int pixelCountMed = 0;
					int counterCountCor = 0;
					int counterCountMed = 0;

					for (i = startImage - 1; i < maxImage; i++) {

						// IJ.write("Bin in Slice: " + i + " Half: " + half);
						if (enabledForSIT[half][i]) {

							// IJ.write("Enabled for SIT Slice: " + i + " Half: " + half);

							averageGrayCortex[half][i] = 0.0;
							averageGrayMedulla[half][i] = 0.0;

							counter = 0;
							j = 0;
							while (j < imageCanvas.MAXROIS && roiCortex != null && roiCortex[i][j] != null) {

								// IJ.write("JOTT Cortex: " + j);
								double[]	roiProperties = new double[6];

								for (int p = 0; p < 6; p++) {
									roiProperties[p] = 0.0;
								}
								roiProperties = getRoiProperties(i + 1, roiCortex[i][j]);
								averageGrayCortex[half][i] += roiProperties[0];
								stdDeviationCortex[half][i] += roiProperties[1];
								pixelCountCor += roiProperties[5];
								j++;
								counter++;
								counterCountCor++;
							}
							if (counter > 1) {
								averageGrayCortex[half][i] = averageGrayCortex[half][i] / (double) counter;
								stdDeviationCortex[half][i] = stdDeviationCortex[half][i] / (double) counter;
							}
							counter = 0;
							j = 0;
							while (j < imageCanvas.MAXROIS && roiMedulla != null && roiMedulla[i][j] != null) {

								// IJ.write("JOTT Medulla: " + j);
								double[]	roiProperties = new double[6];

								for (int p = 0; p < 6; p++) {
									roiProperties[p] = 0.0;
								}
								roiProperties = getRoiProperties(i + 1, roiMedulla[i][j]);
								averageGrayMedulla[half][i] += roiProperties[0];
								stdDeviationMedulla[half][i] += roiProperties[1];
								pixelCountMed += roiProperties[5];
								j++;
								counter++;
								counterCountMed++;
							}
							if (counter > 1) {
								averageGrayMedulla[half][i] = averageGrayMedulla[half][i] / (double) counter;
								stdDeviationMedulla[half][i] = stdDeviationMedulla[half][i] / (double) counter;
							}
						}
					}
					if (counterCountCor > 0) {
						pixelCountCortex[half] = pixelCountCor / counterCountCor;
					}
					if (counterCountMed > 0) {
						pixelCountMedulla[half] = pixelCountMed / counterCountMed;

					}
					if (analysisDialog.interpolateNullValues) {
						averageGray[half] = interpolateNullValues(averageGray[half]);
						averageGrayCortex[half] = interpolateNullValues(averageGrayCortex[half]);
						averageGrayMedulla[half] = interpolateNullValues(averageGrayMedulla[half]);
						stdDeviation[half] = interpolateNullValues(stdDeviation[half]);
						stdDeviationCortex[half] = interpolateNullValues(stdDeviationCortex[half]);
						stdDeviationMedulla[half] = interpolateNullValues(stdDeviationMedulla[half]);
					}
					if (startImage > 1) {
						averageGray[half] = fillMissingValues(averageGray[half]);
						averageGrayCortex[half] = fillMissingValues(averageGrayCortex[half]);
						averageGrayMedulla[half] = fillMissingValues(averageGrayMedulla[half]);
						stdDeviation[half] = fillMissingValues(stdDeviation[half]);
						stdDeviationCortex[half] = fillMissingValues(stdDeviationCortex[half]);
						stdDeviationMedulla[half] = fillMissingValues(stdDeviationMedulla[half]);
					}

					initResultWindow(half);
				}
				if (verbose) {
					IJ.write("Done.");

				}
			}

			setMinMaxAxisValues(JMConst.MAIN);
			setMinMaxAxisValues(JMConst.MEDULLA);
			setMinMaxAxisValues(JMConst.CORTEX);

			IJ.showProgress(1.0);
			IJ.showStatus("");

			this.finished = true;

			// leftAWTImage = ipLeftKidney.getImage();
			// rightAWTImage = ipRightKidney.getImage();
			// dlg.updateImage(JMConst.LEFT,leftAWTImage);
			// dlg.updateImage(JMConst.RIGHT,rightAWTImage);

			sitResults.updateImageLabel(JMConst.LEFT, ipRightKidney);
			sitResults.updateImageLabel(JMConst.RIGHT, ipLeftKidney);

			imagePlus.unlock();
			this.finished = true;													// set finish flag
			this.stopRunning = true;

		} catch (Exception e) {
			IJ.write("Error while calculation SIT-Graph !");
			error = true;
			e.printStackTrace();
			this.stopRunning = true;
			//this.stop();
		}
		if (analysisDialog != null) {
			analysisDialog.sitFinished();
		}
		//this.stop();
	}


	/**
	 * Berechnet der minimalen und maximalen Intensitaetswerte der SIT-Kurven. Diese
	 * werden dem Dialog zur Darstellung der Kurven uebergeben und die dortigen
	 * Koordinatensysteme entsprechend angepasst. So koennen die Kurven immer optimal
	 * zur die Bildschirmanzeige dargestellt werden.
	 * @param tissue die Gewebeart (0=ganze Niere | 1=Medulla | 2=Cortex)
	 */
	public void setMinMaxAxisValues(int tissue) {
		int min = 0;
		int max = 100;

		if (tissue == JMConst.MAIN) {
			min = Integer.MAX_VALUE;
			min = Math.min(getMinArrayVal(substArrays(averageGray[JMConst.LEFT], scaleArray(stdDeviation[JMConst.LEFT], devScaleFactor))), min);
			min = Math.min(getMinArrayVal(substArrays(averageGray[JMConst.RIGHT], scaleArray(stdDeviation[JMConst.RIGHT], devScaleFactor))), min);
			max = Integer.MIN_VALUE;
			max = Math.max(getMaxArrayVal(addArrays(averageGray[JMConst.LEFT], scaleArray(stdDeviation[JMConst.LEFT], devScaleFactor))), max);
			max = Math.max(getMaxArrayVal(addArrays(averageGray[JMConst.RIGHT], scaleArray(stdDeviation[JMConst.RIGHT], devScaleFactor))), max);

			// IJ.write("minY / maxY " + min + " / " + max);
			// sitResults.setYScale(JMConst.MAIN,min,max);
		}
		if (tissue == JMConst.MEDULLA) {
			min = Integer.MAX_VALUE;
			min = Math.min(getMinArrayVal(substArrays(averageGrayMedulla[JMConst.LEFT], scaleArray(stdDeviationMedulla[JMConst.LEFT], devScaleFactor))), min);
			min = Math.min(getMinArrayVal(substArrays(averageGrayMedulla[JMConst.RIGHT], scaleArray(stdDeviationMedulla[JMConst.RIGHT], devScaleFactor))), min);
			max = Integer.MIN_VALUE;
			max = Math.max(getMaxArrayVal(addArrays(averageGrayMedulla[JMConst.LEFT], scaleArray(stdDeviationMedulla[JMConst.LEFT], devScaleFactor))), max);
			max = Math.max(getMaxArrayVal(addArrays(averageGrayMedulla[JMConst.RIGHT], scaleArray(stdDeviationMedulla[JMConst.RIGHT], devScaleFactor))), max);

			// IJ.write("minY / maxY " + min + " / " + max);
			// sitResults.setYScale(JMConst.MEDULLA,min,max);
		}
		if (tissue == JMConst.CORTEX) {
			min = Integer.MAX_VALUE;
			min = Math.min(getMinArrayVal(substArrays(averageGrayCortex[JMConst.LEFT], scaleArray(stdDeviationCortex[JMConst.LEFT], devScaleFactor))), min);
			min = Math.min(getMinArrayVal(substArrays(averageGrayCortex[JMConst.RIGHT], scaleArray(stdDeviationCortex[JMConst.RIGHT], devScaleFactor))), min);
			max = Integer.MIN_VALUE;
			max = Math.max(getMaxArrayVal(addArrays(averageGrayCortex[JMConst.LEFT], scaleArray(stdDeviationCortex[JMConst.LEFT], devScaleFactor))), max);
			max = Math.max(getMaxArrayVal(addArrays(averageGrayCortex[JMConst.RIGHT], scaleArray(stdDeviationCortex[JMConst.RIGHT], devScaleFactor))), max);

			// IJ.write("minY / maxY " + min + " / " + max);
			// sitResults.setYScale(JMConst.CORTEX,min,max);
		}
		sitResults.setYScale(tissue, min, max);
		sitResults.setDevScaleFactor(devScaleFactor, tissue);
	}


	/**
	 * Untersucht alle Pixel innerhalb der uebergebenen ROI auf deren Eigenschaften. Diese
	 * werden in einem Double-Array gespeichert und zurueckgegeben:
	 * <BR> result[0] = mittlere Grauwert
	 * <BR> result[1] = Standardabweichung
	 * <BR> result[2] = Position des Schwerpunktes der ROI (x-Koordinate)
	 * <BR> result[3] = Position des Schwerpunktes der ROI (y-Koordinate)
	 * <BR> result[4] = maximale Grauwert aller Pixel innerhalb der Roi
	 * <BR> result[5] = Anzahl der Pixel innehalb der ROI
	 * @param slice die Bildnummer innerhalb der Sequenz
	 * @param roi die ROI
	 */
	private double[] getRoiProperties(int slice, Roi roi) {
		Rectangle rect = null;

		if (roi != null) {
			rect = roi.getBoundingRect();
		}
		double	averageGray = 0.0;
		double	stdDeviation = 0.0;
		double	centerGravX = 0.0;
		double	centerGravY = 0.0;
		boolean hasMask = false;
		int			maxGray = 0;
		int			pixelCount = 0;
		int			counter, x, y, grayVal;
		int			xSum, ySum;

		short[] pix = (short[]) imageStack.getPixels(slice);

		if (rect != null) {
			int[] mask = roi.getMask();

			if (mask != null) {
				hasMask = true;

				// IJ.write(">>>" + rect.width + "/" + rect.height + "/" + rect.x + "/" + rect.y + "/" + mask.length);
			}
			averageGray = 0.0;
			stdDeviation = 0.0;
			counter = 0;
			xSum = 0;
			ySum = 0;
			maxGray = Integer.MIN_VALUE;
			pixelCount = 0;

			// calculate the average gray of each roi for each slice
			// IJ.write(" Calcualte the average Gray");
			for (x = 0; x < rect.width; x++) {
				for (y = 0; y < rect.height; y++) {
					if ((!hasMask) || ((hasMask) && (mask[(y * (rect.width)) + x] != -1))) {

						// the pixel is part of the roi
						counter++;						// count the No. of pixels in the roi
						grayVal = pix[((y + rect.y) * width) + (x + rect.x)];
						if (grayVal > maxGray) {
							maxGray = grayVal;	// store max gray value
						}
						averageGray = averageGray + (double) grayVal;
						xSum = xSum + ((x + rect.x) * grayVal);
						ySum = ySum + ((y + rect.y) * grayVal);
					}
				}
			}

			// the center of gravity as a co-result (needed later)
			centerGravX = xSum / (int) averageGray;
			centerGravY = ySum / (int) averageGray;
			averageGray = averageGray / (double) counter;
			pixelCount = counter;

			double	sum = 0.0;

			// calculate the standard deviation of each roi of each slice
			// IJ.write(" calculate the standard deviation");
			for (x = 0; x < rect.width; x++) {
				for (y = 0; y < rect.height; y++) {
					if ((!hasMask) || ((hasMask) && (mask[(y * (rect.width)) + x] != -1))) {
						double	d = ((double) pix[((y + rect.y) * width) + (x + rect.x)]) - averageGray;

						sum = sum + (d * d);
					}
				}
			}
			stdDeviation = Math.sqrt((1.0 / ((double) counter)) * sum);
		}

		double[]	result = new double[6];

		result[0] = averageGray;
		result[1] = stdDeviation;
		result[2] = centerGravX;
		result[3] = centerGravY;
		result[4] = maxGray;
		result[5] = pixelCount;
		return result;
	}


	/**
	 * Initialisiert das Dialog-Fenster zur Anzeige der berechneten Kurven.
	 * Die Berechneten Kurvenwerte werden als Arrays uebergeben, ausserdem
	 * werden die Dimensionen der Koordinatensysteme zur Anzeige gesetzt.
	 * @param half die Bildhaelfte, entsprechend der Niere (0=links | 1=rechts)
	 */
	public void initResultWindow(int half) {
		if (sitResults == null) {
			sitResults = new SITResults(this, ipRightKidney, ipLeftKidney);
		}
		sitResults.setAdjustmentRange(adjustmentMin, adjustmentMax);
		sitResults.setAverageGray(half, JMConst.MAIN, averageGray[half]);

		// dlg.setAllAverageGray(JMConst.RIGHT,averageGray[RIGHT]);
		sitResults.setStandardDeviation(half, JMConst.MAIN, stdDeviation[half]);

		// dlg.setAllStandardDeviation(JMConst.RIGHT,stdDeviation[RIGHT]);
		sitResults.setAverageGray(half, JMConst.CORTEX, averageGrayCortex[half]);
		sitResults.setStandardDeviation(half, JMConst.CORTEX, stdDeviationCortex[half]);
		sitResults.setAverageGray(half, JMConst.MEDULLA, averageGrayMedulla[half]);
		sitResults.setStandardDeviation(half, JMConst.MEDULLA, stdDeviationMedulla[half]);

		sitResults.setXScale(JMConst.MAIN, 0, maxTime);
		sitResults.setXScale(JMConst.MEDULLA, 0, maxTime);
		sitResults.setXScale(JMConst.CORTEX, 0, maxTime);

		int marks = averageGray[half].length;

		while (marks > 30) {
			marks = marks / 2;
		}
		sitResults.setXMarks(JMConst.MAIN, marks);
		sitResults.setXMarks(JMConst.MEDULLA, marks);
		sitResults.setXMarks(JMConst.CORTEX, marks);

	}


	/**
	 * Fuehrt die Klassifizeirung der Pixel durch. Benoetigt zur Berechnung zwei
	 * ImagePlus-Klassifizierungsbilder, in denen die Pixel entsprechend der
	 * Klassifizierung eingefaerbt werden sowie zwei Hilfsbilder, die die Entfernung
	 * der Pixel von Schwerpunkt dier Nieren sowie die Werte des mittleren Grauwert-
	 * Anstiegs enthalten.
	 * @param half             die Bildhaelfte, entsprechend der Niere (0=links | 1=rechts)
	 * @param thresholdPelvis  der Schwellenwert zur Einstellung der Gewebegrenze zwischen Nierenbecken und Nierenmark
	 * @param thresholdMedulla der Schwellenwert zur Einstellung der Gewebegrenze zwischen Nierenmark und Nierenrinde
	 * @param thresholdCortex  der Schwellenwert zur Einstellung der Gewebegrenze zwischen Nierenrinde und dem die Niere imgebenden Gewebe
	 */
	public void classifyPixels(int half, double thresholdPelvis, double thresholdMedulla, double thresholdCortex) {
		int[]						mask = null;
		Rectangle				rect = null;
		Roi							roi = null;
		ImageProcessor	ip = null;
		double					dev, dist, asc;
		boolean					classified;
		int							rectWidth, rectHeight, rectX, rectY, val;

		if (half == JMConst.LEFT) {
			rect = lRect;
			roi = lRoi;
			if (ipLeftKidney != null) {
				ip = ipLeftKidney.getProcessor();
			}
		}
		if (half == JMConst.RIGHT) {
			rect = rRect;
			roi = rRoi;
			if (ipRightKidney != null) {
				ip = ipRightKidney.getProcessor();
			}
		}

		ip = (ColorProcessor) ip;
		rectWidth = rect.width;
		rectHeight = rect.height;
		rectX = rect.x;
		rectY = rect.y;

		// iporg = (ColorProcessor)iporg;
		double	inverseThresholdCortex = Math.abs(thresholdCortex - 2.0);
		double	weakThresholdCortex = ((thresholdCortex - 1.0) / 2.0) + 1.0;

		int			lumination = (int) grayAdd;
		int			lumMinus = 0xff101010;
		int			yellowAdd = ((lumination & 0xff) << 16) | ((lumination & 0xff) << 8);

		int			redAdd = ((lumination & 0xff) << 16);
		int			greenAdd = ((lumination & 0xff) << 8);
		int			blueAdd = lumination & 0xff;


		if (rect != null) {
			mask = roi.getMask();

			// ip = new ColorProcessor(rect.width,rect.height);
			if (mask != null) {
				pixelCountCortex[half] = 0;
				pixelCountMedulla[half] = 0;

				// IJ.write("Nue gaehts los");
				for (int x = 0; x < rectWidth; x++) {
					for (int y = 0; y < rectHeight; y++) {
						if (mask[(y * (rectWidth)) + x] != -1) {
							classified = false;
							dev = pixStdDeviation[((y + rectY) * width) + (x + rectX)];
							dist = pixDistancePercent[((y + rectY) * width) + (x + rectX)];
							asc = pixGrayAscent[((y + rectY) * width) + (x + rectX)];

							// dev = dev * 0.65 * thresholdPelvis;
							// dist = dist * 0.40 * thresholdPelvis;
							if ((dev < (maxDeviation[half] * 0.40 * thresholdPelvis)) && (dist < (65.0 * thresholdPelvis))) {
								pixClassification[((y + rectY) * width) + (x + rectX)] = JMConst.PELVIS;
								val = cpClassifyHelper.getPixel((x + rectX), (y + rectY));

								// int val = ip.getPixel(x,y);
								// int gray = iporg.getPixel((x+rect.x),(y+rect.y)) + (int)grayAdd;
								// int val = 0xff000000 | ((gray&0xff)<<16) | ((gray&0xff)<<8) | gray&0xff;
								// val = (val & 0xff00ff00) + greenAdd;
								val = val - lumMinus;
								ip.putPixel(x, y, val);
								classified = true;
							}
							if ((dev < (maxDeviation[half] * 0.75 * thresholdMedulla)) && (dist < (95.0 * thresholdMedulla)) && (asc < (maxGrayAscent[half] * 0.60 * thresholdMedulla)) && (!classified)) {
								pixClassification[((y + rectY) * width) + (x + rectX)] = JMConst.MEDULLA;
								pixelCountMedulla[half] = pixelCountMedulla[half] + 1;
								val = cpClassifyHelper.getPixel((x + rectX), (y + rectY));

								// int val = ip.getPixel(x,y);
								// int gray = iporg.getPixel((x+rect.x),(y+rect.y)) + (int)grayAdd;
								// int val = 0xff000000 | ((gray&0xff)<<16) | ((gray&0xff)<<8) | gray&0xff;
								val = val - lumMinus;
								val = val + redAdd;

								// val = (val & 0xff0000ff) + blueAdd;
								ip.putPixel(x, y, val);
								classified = true;
							}


							/*
							 * if (((dev > (maxDeviation[half] * 0.65 * thresholdCortex)) ||
							 * (dist >= (95.0 * thresholdCortex)) ) &&
							 * (!classified))
							 * {
							 * pixClassification[((y+rectY)*width) + (x+rectX)] = PELVIS;
							 * //pixelCountCortex[half] = pixelCountCortex[half] + 1;
							 * val = cpClassifyHelper.getPixel((x+rectX),(y+rectY));
							 * //int val = ip.getPixel(x,y);
							 * //int gray = iporg.getPixel((x+rect.x),(y+rect.y)) + (int)grayAdd;
							 * //int val = 0xff000000 | ((gray&0xff)<<16) | ((gray&0xff)<<8) | gray&0xff;
							 * //val = val + yellowAdd;
							 * //val = (val & 0xffff0000) + redAdd;
							 * ip.putPixel(x,y,val);
							 * //ip.putPixel(x,y,0xffff0000);
							 * classified = true;
							 * }
							 */

							if (!classified) {
								pixClassification[((y + rectY) * width) + (x + rectX)] = JMConst.CORTEX;
								pixelCountCortex[half] = pixelCountCortex[half] + 1;
								val = cpClassifyHelper.getPixel((x + rectX), (y + rectY));

								// int val = ip.getPixel(x,y);
								// int gray = iporg.getPixel((x+rect.x),(y+rect.y)) + (int)grayAdd;
								// int val = 0xff000000 | ((gray&0xff)<<16) | ((gray&0xff)<<8) | gray&0xff;
								val = val - lumMinus;
								val = val + yellowAdd;

								// val = (val & 0xffff0000) + redAdd;

								// val = 0xffffff00;
								ip.putPixel(x, y, val);
								classified = true;
							}


							if ((dist > (90.0 * weakThresholdCortex)) && (asc < (maxGrayAscent[half] * 0.55 * inverseThresholdCortex)) && (dev < (maxDeviation[half] * 0.50 * inverseThresholdCortex))) {
								pixClassification[((y + rectY) * width) + (x + rectX)] = JMConst.PELVIS;

								// pixelCountCortex[half] = pixelCountCortex[half] + 1;
								val = cpClassifyHelper.getPixel((x + rectX), (y + rectY));

								// int val = ip.getPixel(x,y);
								// int gray = iporg.getPixel((x+rect.x),(y+rect.y)) + (int)grayAdd;
								// int val = 0xff000000 | ((gray&0xff)<<16) | ((gray&0xff)<<8) | gray&0xff;
								// val = val + yellowAdd;
								// val = (val & 0xffff0000) + redAdd;
								val = val - lumMinus;
								ip.putPixel(x, y, val);

								// ip.putPixel(x,y,0xffff0000);
								classified = true;
							}
						}
					}
				}
			}

			if (half == JMConst.LEFT) {

				// ipLeftKidney.setProcessor(ipLeftKidney.getTitle(),ip);
				ipLeftKidney.updateImage();
			}
			if (half == JMConst.RIGHT) {

				// ipRightKidney.setProcessor(ipRightKidney.getTitle(),ip);
				// ipRightKidney = new ImagePlus("RIGHT",ip);
				ipRightKidney.updateImage();
			}

		}
	}


	/**
	 * Berechnet der SIT-Kurven fuer Nierenmark und Nierenrinden entsprechend
	 * der Klassifizierung der Pixel. Da die Ergebnisse der Klassifizierung in einem
	 * zweidimensionalen Array als Klassenvariable angespeichert sind, braucht sie
	 * nicht als Parameter mit uebergeben werden.
	 * @param half             die Bildhaelfte, entsprechend der Niere (0=links | 1=rechts)
	 * @param int              der Gewebetyp
	 */
	public void calcPixelCurve(int half, int tissue) {
		Rectangle rect = null;

		if (half == JMConst.LEFT) {
			rect = lRect;
		}
		if (half == JMConst.RIGHT) {
			rect = rRect;

		}
		if ((rect != null) && (pixClassification != null)) {
			for (int i = startImage; i <= maxImage; i++) {
				if (tissue == JMConst.MEDULLA) {
					averageGrayMedulla[half][i - 1] = 0.0;
					stdDeviationMedulla[half][i - 1] = 0.0;
				}
				if (tissue == JMConst.CORTEX) {
					averageGrayCortex[half][i - 1] = 0.0;
					stdDeviationCortex[half][i - 1] = 0.0;
				}

				if (enabledForSIT[half][i - 1]) {
					pix = (short[]) imageStack.getPixels(i);
					int			counter = 0;
					double	sum = 0;
					double	d = 0;
					double	avGray = 0;

					// calculate average gray
					for (int x = 0; x < rect.width; x++) {
						for (int y = 0; y < rect.height; y++) {
							short pixClass = pixClassification[((y + rect.y) * width) + (x + rect.x)];

							if (pixClass == tissue) {
								counter++;
								short grayVal = pix[((y + rect.y) * width) + (x + rect.x)];

								if (tissue == JMConst.MEDULLA) {
									averageGrayMedulla[half][i - 1] = averageGrayMedulla[half][i - 1] + (double) grayVal;
								}
								if (tissue == JMConst.CORTEX) {
									averageGrayCortex[half][i - 1] = averageGrayCortex[half][i - 1] + (double) grayVal;
								}
							}
						}
					}
					if (tissue == JMConst.MEDULLA) {
						averageGrayMedulla[half][i - 1] = averageGrayMedulla[half][i - 1] / (double) counter;
						avGray = averageGrayMedulla[half][i - 1];
					}
					if (tissue == JMConst.CORTEX) {
						averageGrayCortex[half][i - 1] = averageGrayCortex[half][i - 1] / (double) counter;
						avGray = averageGrayCortex[half][i - 1];
					}

					// calculate standard deviation
					for (int x = 0; x < rect.width; x++) {
						for (int y = 0; y < rect.height; y++) {
							short pixClass = pixClassification[((y + rect.y) * width) + (x + rect.x)];

							if (pixClass == tissue) {
								d = ((double) pix[((y + rect.y) * width) + (x + rect.x)]) - avGray;
								sum = sum + (d * d);
							}
						}
					}
					if (tissue == JMConst.MEDULLA) {
						stdDeviationMedulla[half][i - 1] = Math.sqrt((1.0 / ((double) counter)) * sum);
					}
					if (tissue == JMConst.CORTEX) {
						stdDeviationCortex[half][i - 1] = Math.sqrt((1.0 / ((double) counter)) * sum);
					}
				}		// for
 			 }
		}
		if (analysisDialog.interpolateNullValues) {
			if (tissue == JMConst.MEDULLA) {
				averageGrayMedulla[half] = interpolateNullValues(averageGrayMedulla[half]);
				stdDeviationMedulla[half] = interpolateNullValues(stdDeviationMedulla[half]);
			}
			if (tissue == JMConst.CORTEX) {
				averageGrayCortex[half] = interpolateNullValues(averageGrayCortex[half]);
				stdDeviationCortex[half] = interpolateNullValues(stdDeviationCortex[half]);
			}
		}
		if (startImage > 1) {
			if (tissue == JMConst.MEDULLA) {
				averageGrayMedulla[half] = fillMissingValues(averageGrayMedulla[half]);
				stdDeviationMedulla[half] = fillMissingValues(stdDeviationMedulla[half]);
			}
			if (tissue == JMConst.CORTEX) {
				averageGrayCortex[half] = fillMissingValues(averageGrayCortex[half]);
				stdDeviationCortex[half] = fillMissingValues(stdDeviationCortex[half]);
			}
		}
	}


	/**
	 * Hilfsmethode. Liefert den minimalen Wert innerhalb des uebergebenen Double-Arrays
	 * zurueck, zur naechsten ganzen Zahl (Integer) abgerundet.
	 * @param data Array von double-Werten
	 * @return int der minimale Wert
	 */
	private int getMinArrayVal(double[] data) {
		double	dmin = Double.MAX_VALUE;

		for (int i = 0; i < data.length; i++) {
			if (data[i] < dmin) {
				dmin = data[i];
			}
		}
		return (int) Math.floor(dmin);
	}


	/**
	 * Hilfsmethode. Liefert den maximalen Wert innerhalb des uebergebenen Double-Arrays
	 * zurueck, zur naechsten ganzen Zahl (Integer) aufgerundet.
	 * @param data Array von double-Werten
	 * @return int der maximale Wert
	 */
	private int getMaxArrayVal(double[] data) {
		double	dmax = Double.MIN_VALUE;

		for (int i = 0; i < data.length; i++) {
			if (data[i] > dmax) {
				dmax = data[i];
			}
		}
		return (int) Math.ceil(dmax);
	}


	/**
	 * Hilfsmethode. Liefert die Position innerhalb des uebergebenen Double-Arrays
	 * zurueck, die den maximalen Wert enthaelt.
	 * @param data Array von double-Werten
	 * @return int Postion (Index) des groessten Arrayeintrags
	 */
	private int getMaxArrayPos(double[] data) {
		double	dmax = Double.MIN_VALUE;
		int			pos = 0;

		for (int i = 0; i < data.length; i++) {
			if (data[i] > dmax) {
				dmax = data[i];
				pos = i;
			}
		}
		return pos;
	}


	/**
	 * Hilfsmethode. Berechnet den gleitenden Mittelwert fuer ein Double-Array,
	 * dies entspricht einer Glaettung der durch das Array dargestellten Kurve.
	 * @param data Array von double-Werten
	 * @param range legt fest, wieviele der Nachbarwerte eines betrachteten Arrayeintrags fuer die Mittelwertbildung einbezogen werden sollen.
	 * @return double das geglaettete Array
	 */
	private double[] floatingAverage(double[] data, int range) {
		double[]	result = new double[data.length];
		double		floatingMid;
		int				i, j, pos;

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

			}
			floatingMid = floatingMid / (double) range;
			result[i] = floatingMid;
		}
		return result;
	}


	/**
	 * Addiert zwei Double-Array stellenweise.
	 * Die Arrays sollten die gleiche Laenge haben.
	 * @param data1 Array von double-Werten
	 * @param data2 Array von double-Werten
	 * @return double das Ergebnis der Addition
	 */
	private double[] addArrays(double[] data1, double[] data2) {
		double[]	result = new double[data1.length];

		for (int i = 0; i < data1.length; i++) {
			result[i] = data1[i] + data2[i];
		}
		return result;
	}


	/**
	 * Subtrahiert zwei Double-Array stellenweise.
	 * Die Arrays sollten die gleiche Laenge haben.
	 * @param data1 Array von double-Werten
	 * @param data2 Array von double-Werten
	 * @return double das Ergebnis der Subtraktion.
	 */
	private double[] substArrays(double[] data1, double[] data2) {
		double[]	result = new double[data1.length];

		for (int i = 0; i < data1.length; i++) {
			result[i] = data1[i] - data2[i];
			if (result[i] < 0) {
				result[i] = 0;

				// IJ.write("data1 / data2 / result " + data1[i] + " / " + data2[i] + " / " + result[i]);
			}
		}
		return result;
	}


	/**
	 * Multipliziert alle Eintraege eines Double-Arrays
	 * mit einem angegebenen Faktor.
	 * @param data Array von double-Werten
	 * @param scaler der Multiplikator
	 * @return double das Ergebnis der Berechnung
	 */
	private double[] scaleArray(double[] data, double scaler) {
		double[]	result = new double[data.length];

		for (int i = 0; i < data.length; i++) {
			result[i] = data[i] * scaler;
			if (result[i] < 0) {
				result[i] = 0;
			}
		}
		return result;
	}


	/**
	 * Fuellt Nullwerte am Anfang und Ende eines Double-Arrays mit sinnvollen
	 * Werten auf. Diese werden aus der jeweils ersten oder letzten belegten
	 * Arrayposition entnommen: Beispiel: 0, 0, 0, 12.4 -> 12.4, 12.4, 12.4, 12.4
	 * @param data Array von double-Werten
	 * @return double das Ergebnis der Berechnung
	 */
	private double[] fillMissingValues(double[] data) {
		double[]	result = new double[data.length];
		double		firstVal, lastVal;
		int				i, j, pos;

		for (i = 0; i < data.length; i++) {
			result[i] = data[i];
		}
		i = 0;
		while (i < data.length && data[i] < 0.001) {
			i++;
		}
		if (i == data.length) {
			return result;
		}
		if (i > 0) {
			firstVal = data[i];
			for (pos = i - 1; pos >= 0; pos--) {
				result[pos] = firstVal;
			}
		}
		i = 0;
		while (i < data.length && result[i] > 0.001) {
			i++;
		}
		if (i == data.length) {
			return result;
		}
		lastVal = data[i - 1];
		for (pos = i; pos < data.length; pos++) {
			result[pos] = lastVal;
		}
		return result;
	}


	/**
	 * Berechnet die (diskrete) erste Ableitung eines Double-Arrays.
	 * @param data Array von double-Werten
	 * @return double das Ergebnis der Berechnung
	 */
	private double[] derivation(double[] data) {
		double[]	result = new double[data.length];
		double		gradePoint;
		double		gradePrev;
		double		gradeNext;

		int				i, j, pos;

		for (i = 1; i < (data.length - 1); i++) {
			gradePrev = (data[i] - data[i - 1]);
			gradeNext = (data[i + 1] - data[i]);
			gradePoint = (gradePrev + gradeNext) / 2.0;
			result[i] = gradePoint;
		}
		return result;

	}


	/**
	 * Fuellt Nullwerte innerhalb eines Double-Arrays mit sinnvollen
	 * Werten auf. Diese werden durch Interpolation aus den Nachbareintraegen
	 * erzeugt. Beispiel: 12.7, 0, 16.1 -> 12.7, 14.4, 16.1
	 * @param data Array von double-Werten
	 * @return double das Ergebnis der Berechnung
	 */
	private double[] interpolateNullValues(double[] data) {
		double[]	result = new double[data.length];
		double		ascent;
		int				i, j, pos;

		for (i = startImage - 1; i < data.length; i++) {
			result[i] = data[i];
			if (data[i] < 0.001) {
				j = i + 1;
				while (j < data.length && data[j] < 0.001) {
					j++;
				}
				if (j < data.length) {
					ascent = (data[j] - data[i - 1]) / (double) (j - (i - 1));

					// IJ.write ("ascent: " + ascent + " -> from: " + data[j] + "-" + data[i-1] + " / " + j + "-" + (i-1));
					double	dcounter = 1.0;

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

						// IJ.write ("pos: " + pos + " !!! " + (ascent * dcounter) + "222: " + data[pos-1]);
						result[pos] = data[pos - 1] + (ascent * dcounter);
						dcounter = dcounter + 1.0;
					}
				}
				i = j - 1;
			}
		}
		return result;
	}


	/**
	 * 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)));
	}


	/**
	 * Setzt den Wertebereich, innerhalb dem die Schwellenwerte veraendert werden koennen.
	 * Die Schwellenwerte werden  mit 1,0 intialisiert. Der min- und max-Wert wird durch
	 * die Einstellungen im Voreinstellungs-Dialog (jm.kidney.PresetsDialog) festgelegt.
	 * @param min   der minimal erlaubte Wert
	 * @param max   der maximal erlaubte Wert
	 */
	public void setAdjustmentRange(double min, double max) {
		adjustmentMin = min;
		adjustmentMax = max;
	}


	/**
	 * Hilfsmethode. Legt fest, ob Zwischenergebnisse in textueller Form im
	 * ImageJ-Hauptfenster ausgegeben werden sollen.
	 * @param v     true, wenn Text ausgegeben werden soll.
	 */
	public void setVerbose(boolean v) {
		verbose = v;
	}

}


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

