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

/*
 * The package jm.cepstrum 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.cepstrum;

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

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

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


/**
 * Realisiert den Algothimus zur Kepstrum-Filterung der Nierensequenzen.
 * Wird von der Klasse AnalysisDialog aufgerufen und benoetigt die Klasse
 * Cepstrum zur Durchfuehrung der Berechnungen.
 * @see jm.kidney.AnalysisDialog
 * @see jm.cepstrum.Cepstrum
 * @version  3.6, 13/05/2005, tha: Kleine Aenderungen in run und applyShift.
 * @author   Jens Martin
 */
public class CepstrumAnalysis extends Thread {
	private int							maxImage = 0;
	private int							startImage = 1;
	private boolean					verbose = true;
	private boolean					secondChance = true;
	private boolean					writeLog = true;
	private boolean					newWindow = false;
	private int							minDiffSwitch = 1000;
	private java.awt.Image	referenceImage = null;
	private java.awt.Image	lastImage = null;
	private java.awt.Image	actualImage = null;
	private AnalysisDialog	analysisDialog = null;
	private ImagePlus				imagePlus = null;
	private ImageStack			imageStack = null;
	private ImageStack			newStack = 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					analysisAccuracy = 1.0;
	private int							min_x_left;
	private int							min_y_left;
	private int							max_x_left;
	private int							max_y_left;
	private int							min_x_right;
	private int							min_y_right;
	private int							max_x_right;
	private int							max_y_right;
	private int							lastChangeLeft = 0;
	private int							lastChangeRight = 0;

	private double					minImageDiff;

	public String						seriesDate = "";
	public String						studyDate = "";
	public String						patientName = "";
	public String						patientBirth = "";
	public String						patientWeight = "";

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

	public int[]						cepstrumQuality = 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
	 * @param rightRoi Die rechte Roi
	 * @param leftRoi Die linke Roi
	 */
	public CepstrumAnalysis( /* ImagePlus */AnalysisDialog dlg, Roi rightRoi, Roi leftRoi) {
		analysisDialog = dlg;
		imagePlus = analysisDialog.imagePlus;

		// imagePlus = (ImagePlus)img;
		imageStack = imagePlus.getStack();
		maxImage = imagePlus.getStackSize();
		height = imagePlus.getHeight();
		width = imagePlus.getWidth();
		min_y_left = height + 1;
		min_y_right = height + 1;
		max_y_left = -1;
		max_y_right = -1;
		analysisAccuracy = 1.0;
		lastChangeLeft = 1;
		lastChangeRight = 1;
		startImage = analysisDialog.startImage;

		cepstrumQuality = new int[maxImage];

		if (leftRoi != null) {
			lRoi = leftRoi;
			lRect = lRoi.getBoundingRect();
			if (verbose) {
				IJ.write("left ROI ---  x:" + lRect.x + " y: " + lRect.y + " w: " + lRect.width + " h: " + lRect.height);
			} 
		} 
		if (rightRoi != null) {
			rRoi = rightRoi;
			rRect = rRoi.getBoundingRect();
			if (verbose) {
				IJ.write("right ROI ---  x:" + rRect.x + " y: " + rRect.y + " w: " + rRect.width + " h: " + rRect.height);
			} 
		} 


	}


	/**
	 * Hilfsmethode. Setzt die Genauigkeit der Kepstrum-Filterung.
	 * @param accuracy double Wert fuer die Genauigkeit. Zulaessiger Bereich: 1.0-4.0
	 */
	public void setAccuracy(double accuracy) {
		analysisAccuracy = accuracy;
	} 


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


	/**
	 * Hilfsmethode. Legt fest, ob die korrigierte Nierensequenz wieder
	 * in das alte Fenster zurueckgeschrieben werden soll (und damit
	 * die alte Sequenz verworfen wird) oder ob ein neues Fenster zur Anzeige
	 * verwendet werden soll.
	 * @param w     true, wenn ein neues Fenster gewuenscht ist.
	 */
	public void newWindow(boolean w) {
		newWindow = w;
	} 


	/**
	 * Hilfsmethode. Legt fest, ob der Algorithmus zur Filterung der Sequenz
	 * fehlertolerant arbeiten soll.
	 * @param sc     true = Fehlertoleranz einschalten
	 */
	public void setSecondChance(boolean sc) {
		secondChance = sc;
	} 


	/**
	 * Hilfsmethode. Legt fest, ob Zwischenergebnisse textuell in eine Protokolldatei
	 * (LogFile) geschrieben werden sollen.
	 * @param log     true, wenn Text in die Dateiausgegeben werden soll.
	 */
	public void setWriteLog(boolean log) {
		writeLog = log;
	} 


	/**
	 * Hilfsmethode. Legt fest, nach wie vielen Iterationen des Algorithmus ohne
	 * Aenderung der maximalen Verschiebungsachse auf den minimalen Differenzfilter
	 * umgeschaltet werden soll.
	 * @param s     gibt die Anzahl der Iterationen an. true, wenn Text ausgegeben werden soll.
	 */
	public void setMinDiff(int s) {
		minDiffSwitch = s;
	} 


	/**
	 * 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, y;
		boolean			success;
		ColorModel	cm = imageStack.getColorModel();

		newStack = new ImageStack(width, height, cm);

		// mark: alle ersten Images bis zum Startimage einfach kopieren
		for (i = 1; i <= startImage; i++) {
			newStack.addSlice(imageStack.getSliceLabel(i), imageStack.getPixels(i));
		} 
		java.awt.Point	shift = new java.awt.Point(0, 0);
		java.awt.Point	firstShift = new java.awt.Point(0, 0);
		java.awt.Point	secondShift = new java.awt.Point(0, 0);
		double					firstMinDiff = 4.0;
		double					secondMinDiff = 4.0;
		Cepstrum				cepstrumFilter = new Cepstrum();

		cepstrumFilter.showImages = false;
		cepstrumFilter.setVerbose(verbose);
		cepstrumFilter.setWriteLog(writeLog);

		ByteProcessor leftIP_ref = null;
		ByteProcessor leftIP_act = null;
		ByteProcessor rightIP_ref = null;
		ByteProcessor rightIP_act = null;

		analysisDialog.setCepstrumResult(JMConst.LEFT, startImage, 0);
		analysisDialog.setCepstrumResult(JMConst.RIGHT, startImage, 0);

		try {
      // tha 2002.5.13: In der aktuellen Version von ImageJ (1.26) werden keine 
      // unsigned Images mehr unterstuetzt.
			// boolean unsigned = imagePlus.getProcessor().isUnsigned();

			if (lRect != null) {
				int lrw = (int) Math.ceil((double) lRect.width * analysisAccuracy);
				int lrh = (int) Math.ceil((double) lRect.height * analysisAccuracy);

				leftIP_ref = new ByteProcessor(lrw, lrh, new byte[lrw * lrh], cm);
				leftIP_act = new ByteProcessor(lrw, lrh, new byte[lrw * lrh], cm);
			} 
			if (rRect != null) {
				int rrw = (int) Math.ceil((double) rRect.width * analysisAccuracy);
				int rrh = (int) Math.ceil((double) rRect.height * analysisAccuracy);

				rightIP_ref = new ByteProcessor(rrw, rrh, new byte[rrw * rrh], cm);
				rightIP_act = new ByteProcessor(rrw, rrh, new byte[rrw * rrh], cm);
			} 
			if (!(lRect != null || rRect != null)) {
				return;

			} 
			IJ.showProgress(0.001);
			int maxProgress = 0;
			int counter = 0;

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

			// Grosse While-Schleife zum durchlaufen der Bilder;
			i = startImage + 1;
			while (!(stopRunning == true || i > maxImage)) {
				for (int half = JMConst.LEFT; half <= JMConst.RIGHT; half++) {
					Rectangle			rect;
					ByteProcessor halfIP_ref;
					ByteProcessor halfIP_act;
					String				strHalf;
					int						lastChangeHalf;

					if (half == JMConst.LEFT) {
						rect = lRect;
						strHalf = "left";
						halfIP_ref = leftIP_ref;
						halfIP_act = leftIP_act;
						lastChangeHalf = lastChangeLeft;
					} else {
						rect = rRect;
						strHalf = "right";
						halfIP_ref = rightIP_ref;
						halfIP_act = rightIP_act;
						lastChangeHalf = lastChangeRight;
					} 

					// Thread: damit auch andere mal drankommen
					try {
						this.sleep(5);
					} catch (InterruptedException ie) {
						IJ.write("Interrupted in Ceptrum-Filter !");
					} 

					// here starts the true work
					if (rect != null) {

						// IJ.showStatus(" Image " + i + " / " + maxImage + " - Cepstrum filter left kidney");
						if (writeLog) {
							ResultLog logFile = new ResultLog("cepstrum.log");

							logFile.append("" + patientName + "|" + patientBirth + "|" + patientWeight + "|" + studyDate + "|" + seriesDate + "|" + i + "|" + strHalf + "|" + analysisAccuracy + "|");
						} 
						if (verbose) {
							IJ.write("__________________________________________");
						} 
						if (verbose) {
							IJ.write("Movement Correction for Image " + i + " - " + strHalf);

							// Prepare the Images for Cepstrum-Filter or linearTranslationCorrection
							// todo: equalize histograms of the both images
						} 
						boolean firstTry = true;
						boolean secondTry = false;
						int			tryCounter = 0;

						while (firstTry || secondTry) {
							tryCounter++;
							int lastGoodResult = i - 1;

							while ((cepstrumQuality[lastGoodResult] > 4) && (lastGoodResult > 1)) {
								lastGoodResult--;
							} 
							if ((secondTry) && (lastGoodResult > 1)) {
								lastGoodResult--;
								while ((cepstrumQuality[lastGoodResult] > 4) && (lastGoodResult > 1)) {
									lastGoodResult--;
								} 
							} 
							prepareImages(half, halfIP_ref, lastGoodResult, halfIP_act, i);

							// Histogramm-Angleich:
							halfIP_act = histogramEquilisation(halfIP_ref, halfIP_act);
							firstTry = false;

							// Can we switch to linear ???
							if ((i - lastChangeHalf) < minDiffSwitch) {
								IJ.showStatus(" Image " + i + " / " + maxImage + " - Cepstrum-Filter " + strHalf + " kidney");

								// Create java.awt.Images for the Cepstrum Class
								referenceImage = halfIP_ref.createImage();
								actualImage = halfIP_act.createImage();

								// Start the cepstrum filter
								success = cepstrumFilter.calculateTranslation(referenceImage, actualImage, shift);
								if (writeLog) {
									ResultLog logFile = new ResultLog("cepstrum.log");

									logFile.append("" + shift.x + "|" + shift.y + "|" + ((double) shift.x / analysisAccuracy) + "|" + ((double) shift.y / analysisAccuracy) + "|");
								} 

								// Check the result..
								int range = (int) Math.round(analysisAccuracy * 2.0);

								if (range > 7) {
									range = 7;

									// minImageDiff = checkCepstrumResult(filterKBWindow(filterSobel(halfIP_ref),true),filterKBWindow(filterSobel(halfIP_act),true),range,2,shift);
								} 
								minImageDiff = checkCepstrumResult(halfIP_ref, halfIP_act, range, 2, shift);
							} else {
								IJ.showStatus(" Image " + i + " / " + maxImage + " - Minimum-Difference filter " + strHalf + " kidney");
								linearTranslationCorrection(halfIP_ref, halfIP_act, half, shift);
							} 

							// Thread: others' turn
							try {
								this.sleep(5);
							} catch (InterruptedException ie) {
								IJ.write("Interrupted in Ceptrum-Filter !");
							} 

							// Set Progress bar
							counter++;
							IJ.showProgress((double) counter / (double) maxProgress);
							if (secondTry) {
								secondMinDiff = minImageDiff;
								secondShift.x = shift.x;
								secondShift.y = shift.y;
							} else {
								firstMinDiff = minImageDiff;
								firstShift.x = shift.x;
								firstShift.y = shift.y;
							} 
							cepstrumQuality[i - 1] = 5;
							if ((cepstrumFilter.peakQuality <= 20.0) && (cepstrumFilter.peakAboveAverage >= 150.0) && (minImageDiff <= 21.0)) {
								cepstrumQuality[i - 1] = 4;
							} 
							if ((cepstrumFilter.peakQuality <= 13.0) && (cepstrumFilter.peakAboveAverage >= 200.0) && (minImageDiff <= 14.0)) {
								cepstrumQuality[i - 1] = 3;
							} 
							if ((cepstrumFilter.peakQuality <= 10.0) && (cepstrumFilter.peakAboveAverage >= 200.0) && (minImageDiff <= 11.0)) {
								cepstrumQuality[i - 1] = 2;
							} 
							if ((cepstrumFilter.peakQuality <= 8.0) && (cepstrumFilter.peakAboveAverage >= 300.0) && (minImageDiff <= 7.0)) {
								cepstrumQuality[i - 1] = 1;
							} 
							analysisDialog.setCepstrumResult(half, i, cepstrumQuality[i - 1]);

							if ((cepstrumQuality[i - 1] == 5) && (secondTry == false)) {
								secondTry = true;
							} else {
								secondTry = false;
							} 
						} 
						if (tryCounter == 2) {

							// pruefen, welcher der beiden Punkte weniger schlecht ist...
							IJ.write("DRIN1");
							Point temp = compareShifts(firstShift, firstMinDiff, secondShift, secondMinDiff);

							shift.x = temp.x;
							shift.y = temp.y;
						} 
						applyShift(i, half, shift);
						if (writeLog) {
							ResultLog logFile = new ResultLog("cepstrum.log");

							logFile.append("" + shift.x + "|" + shift.y + "|" + ((double) shift.x / analysisAccuracy) + "|" + ((double) shift.y / analysisAccuracy) + "\r\n");
						} 
						counter++;
						IJ.showProgress((double) counter / (double) maxProgress);
					} 
				} 
				i++;
			} 
		} catch (Exception e) {
			error = true;
			e.printStackTrace();
			IJ.write("Error in Movement Correction !");
		} 
		if (i > maxImage) {
			this.finished = true;		// set finish flag

			// delete old ImageStack
			for (i = 1; i <= maxImage; i++) {
				imageStack.deleteLastSlice();
			} 

			// create new Image Stack
			for (i = 1; i <= maxImage; i++) {
				imageStack.addSlice(newStack.getSliceLabel(i), newStack.getPixels(i));
			} 
		} 

		// imagePlus.updateImage();
		imagePlus.updateAndDraw();

		if (verbose) {
			IJ.write("Movement Correction finished");
		} 
		IJ.showProgress(1.0);
		IJ.showStatus("");
		this.stopRunning = true;

		if (cepstrumFilter.showImages) {
			cepstrumFilter.fftStack.setColorModel(cm);
			cepstrumFilter.fft2Stack.setColorModel(cm);
			cepstrumFilter.psStack.setColorModel(cm);
			cepstrumFilter.ps2Stack.setColorModel(cm);
			cepstrumFilter.orgStack.setColorModel(cm);

			ImagePlus test1 = new ImagePlus("fftImage", cepstrumFilter.fftStack);
			ImagePlus test2 = new ImagePlus("fft2Image", cepstrumFilter.fft2Stack);
			ImagePlus test3 = new ImagePlus("psImage", cepstrumFilter.psStack);
			ImagePlus test4 = new ImagePlus("ps2Image", cepstrumFilter.ps2Stack);
			ImagePlus test5 = new ImagePlus("orgImage", cepstrumFilter.orgStack);

			test1.show();
			IJ.wait(500);
			test2.show();
			IJ.wait(500);
			test3.show();
			IJ.wait(500);
			test4.show();
			IJ.wait(500);
			test5.show();
			IJ.wait(500);
		} 

		// c.reshape( 100,100,910,650);
		// c.show();
		// cepstrumFilter = null;
		newStack = null;
		analysisDialog.cepstrumFinished();
		imagePlus.unlock();
	} 


	/**
	 * Hilfsmethode. Liefert die bessere von zwei moeglichen, berechneten Verschiebungen
	 * zurueck. Wird von der Methode run benoetigt, um das SecondChance-Konzept
	 * zur fehlertoleranten Kepstrum-Filterung zu realisierten.
	 * @param p1 Verschiebungsvektor 1
	 * @param m1 mittlerer Grauwert des Differenzbildes 1
	 * @param p2 Verschiebungsvektor 2
	 * @param m1 mittlerer Grauwert des Differenzbildes  1
	 * @return Point der bessere Verschiebungsvektor
	 */
	private Point compareShifts(Point p1, double m1, Point p2, double m2) {

		// IJ.write("DRIN2");
		double	tr = 22.0;
		double	e1 = Math.abs((double) p1.x / analysisAccuracy) + Math.abs((double) p1.y / analysisAccuracy);
		double	e2 = Math.abs((double) p2.x / analysisAccuracy) + Math.abs((double) p2.y / analysisAccuracy);

		// IJ.write("WERTE:" + e1 + "!" + e2);
		if (e1 < tr && e2 > tr) {
			return p1;
		} 
		if (e1 > tr && e2 < tr) {
			return p2;
		} 
		if (e1 < e2) {
			return p1;
		} else {
			return p2;
		}
	} 


	/**
	 * Verschiebt ein angegebenes Halbbild der Sequenz. Es wird der ganzzahlige
	 * Verschiebungsvektor uebergeben, der noch durch die Skalierungs geteilt wird.
	 * Ist das Ergebnis nicht mehr ganzzahlig, wird automatisch bikubische Interpolation
	 * zur Verschiebung angewendet.
	 * @param slice die Schicht der Bilderstapels, entsp. der Bildnummer der Sequenz
	 * @param half  die Bildhaelfte
	 * @param shift der Verschiebungsvektor
	 */
	private void applyShift(int slice, int half, java.awt.Point shift) {
    // tha 2002.5.13: In der aktuellen Version von ImageJ (1.26) wird  
    // ImageStack.setSlice(int slice) nicht mehr unterstuetzt.
		// imageStack.setSlice(slice);

		// IJ.write ("Setze Slice: " + slice);
		// IJ.write ("Lese  Slice: " +  imageStack.getCurrentSlice() );

		short pix_org[] = new short[width * height];

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

		short pix[] = new short[width * height];

		System.arraycopy(pix_org, 0, pix, 0, width * height);

		// short pix[][] = new short [width*height*2][];
		// pix[0] = (short[]) imageStack.getPixels(slice);
		// pix[1] = (short[]) imageStack.getPixels(slice);
		// short sobel[] = new short [width*height];
		int			x, y;


		// int shift_x = shift.x / accuracy;
		// int shift_y = shift.y / accuracy;

		double	shift_x = (double) shift.x / analysisAccuracy;
		double	shift_y = (double) shift.y / analysisAccuracy;

		if (half == JMConst.LEFT) {
			boolean change = false;

			if (shift.y < min_y_left) {
				min_y_left = shift.y;
				min_x_left = shift.x;
				change = true;
			} 
			if (shift.y > max_y_left) {
				max_y_left = shift.y;
				max_x_left = shift.x;
				change = true;
			} 
			if (change == true) {
				lastChangeLeft++;
			} 
		} else {
			boolean change = false;

			if (shift.y < min_y_right) {
				min_y_right = shift.y;
				min_x_right = shift.x;
				change = true;
			} 
			if (shift.y > max_y_right) {
				max_y_right = shift.y;
				max_x_right = shift.x;
				change = true;
			} 
			if (change == true) {
				lastChangeRight++;
			} 
		} 

		if (verbose) {
			IJ.write("shifting Images, shift: x=" + shift_x + " y=" + shift_y);

			// int halfShift_x = shift.x - (shift_x * accuracy);
			// int halfShift_y = shift.y - (shift_y * accuracy);

		} 
		int start_x, stop_x, adder_x, start_y, stop_y, adder_y, start_half, stop_half;

		if (half == JMConst.LEFT) {
			start_half = 0;
			stop_half = width / 2;
		} else {
			start_half = width / 2;
			stop_half = width;
		} 

		start_x = start_half;
		stop_x = stop_half;
		adder_x = 1;
		start_y = 0;
		stop_y = height;
		adder_y = 1;

		if (shift.x > 0) {
			start_x = stop_half - 1;	// (width / 2) - 1;
			stop_x = start_half - 1;	// (-1);
			adder_x = (-1);
		} 
		if (shift.y > 0) {
			start_y = height - 1;
			stop_y = (-1);
			adder_y = (-1);
		} 

		// IJ.write("ApplyShift: Vorbereitung fertig!");

		if (!(analysisAccuracy > 1.0)) {
			y = start_y;
			while (y != stop_y) {
				x = start_x;
				while (x != stop_x) {
					short val = pix[(y * width) + x];
					int		px = x + (int) Math.round(shift_x);
					int		py = y + (int) Math.round(shift_y);		// - 1;

					if (!(px < start_half || py < 0 || px >= stop_half || py >= height)) {
						pix[(py * width) + px] = val;
					} 
					x = x + adder_x;
				} 
				y = y + adder_y;
			} 
		} else {
			BicubicInterpolation	bic = new BicubicInterpolation(pix, width, height);

			y = start_y;
			while (y != stop_y) {
				x = start_x;
				while (x != stop_x) {
					int val = bic.getPixel((double) x - shift_x, (double) y - shift_y);		// [0][(y*256)+x];

					pix[(y * width) + x] = (short) val;
					x = x + adder_x;
				} 
				y = y + adder_y;
			} 
		} 

		if (newStack.getSize() == slice) {
			newStack.deleteLastSlice();
		} 
		newStack.addSlice(imageStack.getSliceLabel(slice), pix);

    // tha 2002.5.13: In der aktuellen Version von ImageJ (1.26) wird  
    // ImageStack.setSlice(int slice) nicht mehr unterstuetzt.
		// imageStack.setSlice(slice);
    // tha 2002.5.13: In der aktuellen Version von ImageJ (1.26) wird  
    // ImageStack.setPixel(Object pix) nicht mehr unterstuetzt. Statt dessen
    // muss jetzt eine Slice-Nummer mit angegeben werden.
		// imageStack.setPixels(pix);
		imageStack.setPixels(pix, slice);
	} 


	/**
	 * Ueberprueft das Ergebnis der Kepstrum-Filterung numerisch und fuehrt
	 * eine Korrektur durch, falls notwendig.
	 * @param ip_ref  der Imageprocessor des Referenzbildes
	 * @param ip_act  der Imageprocessor des Folgebildes
	 * @param range   die Groesse des tu testetenden Korrekturfensters
	 * @param accuracy   die Genaugigkeit
	 * @param shift   die vom Kepstrum berechnete Verschiebung
	 * @return double den mittleren Grauwert des Differenzbildes zwischen Refezen- und Folgebild an der besten Verschiebungsposition
	 */
	double checkCepstrumResult(ImageProcessor ip_ref, ImageProcessor ip_act, int range, int accuracy, java.awt.Point shift) {

		// boolean result = false;
		// check all the input-parameters
		// if (half > 1 || half < 0) return result;
		if (range < 1 || accuracy < 1) {
			return -1.0;
		} 
		if (verbose) {
			IJ.write("Checking result...");

			// imageStack.setSlice(slice);

		} 
		short pix[][] = new short[ip_ref.getWidth() * ip_ref.getHeight() * 2][];

		// IJ.write ("Lese Pixel des Ref-IP");
		// pix[0] = (short[]) ip_ref.getPixelsCopy();
		// IJ.write ("Lese Pixel des Acz-IP");
		// pix[1] = (short[]) ip_act.getPixelsCopy();
		// IJ.write ("Rechne mir nun einen...");
		int		startx = range;
		int		endx = ip_ref.getWidth() - range;
		int		starty = range;
		int		endy = ip_ref.getHeight() - range;

		if (shift.x < 0) {
			startx = startx + Math.abs(shift.x);
		} else {
			endx = endx - shift.x;
		} 
		if (shift.y < 0) {
			starty = starty + Math.abs(shift.y);
		} else {
			endy = endy - shift.y;
		} 

		int			x, y, i, j, counter;
		double	add[][] = new double[range * 2 + 1][range * 2 + 1];

		for (j = range * (-1); j <= range; j++) {
			for (i = range * (-1); i <= range; i++) {
				counter = 0;
				add[i + range][j + range] = 0;
				for (y = starty; y < endy; y++) {
					int offset = y % accuracy;

					for (x = startx + offset; x < endx; x = x + accuracy) {
						add[i + range][j + range] = add[i + range][j + range] + (double) Math.abs(ip_ref.getPixel(x, y) - ip_act.getPixel((x - shift.x - i), (y - shift.y - j)));

						// (double)Math.abs( pix[1][(y*width)+x] - pix[0][((y+shift.y+j)*width)+x+shift.x+i] );
						counter++;
					} 
				} 
				add[i + range][j + range] = add[i + range][j + range] / (double) counter;

				// IJ.write ("New val for addon i: " + i + " j: " + j + " -> " + add[i+range][j+range]);

			} 

		} 
		int			min_i = -1;
		int			min_j = -1;
		double	min_val = Double.MAX_VALUE;

		for (j = range * (-1); j <= range; j++) {
			for (i = range * (-1); i <= range; i++) {
				if (add[i + range][j + range] < min_val) {
					min_val = add[i + range][j + range];
					min_i = i;
					min_j = j;
				} 
			} 
		} 
		shift.x = shift.x + min_i;
		shift.y = shift.y + min_j;
		if (verbose) {
			IJ.write("Found MinDiffValue " + add[min_i + range][min_j + range] + " at i: " + min_i + " j: " + min_j + " -> correcting Shift to x:" + shift.x + " y: " + shift.y);
		} 
		ResultLog log07 = new ResultLog("cepstrum.log");

		log07.append("" + add[min_i + range][min_j + range] + "|");

		// return result;
		return add[min_i + range][min_j + range];
	} 


	/**
	 * Berechnet der Deformationsmatrix fuer die Rastertransformation. Wurde aus Zeitgruenden
	 * noch nicht fertig implementiert. Wird z.Zt. nicht aufgerufen.
	 * eine Korrektur durch, falls notwendig.
	 * @param slice   die Bildnummer der Sequenz
	 * @param half    die Bildhaelfte
	 * @param ip_act  der Imageprocessor des Folgebildes
	 * @param shift   die berechnete Verschiebung
	 */
	private void getDeformationMatrix(int slice, int half, ImageProcessor ip_act, java.awt.Point shift) {
		short pix[] = new short[width * height];

		pix = (short[]) imageStack.getPixels(slice);
		BicubicInterpolation	bi_ref = new BicubicInterpolation(pix, width, height);
		int										pix_ref[] = new int[9 * 9];

		int										i, j, x, y, rect_width, rect_height, rect_x, rect_y;
		int										sx, sy, counterX, counterY, matrixX, matrixY;
		double								dx, dy;
		short									val_ref, val_act, val_diff;

		if (half == JMConst.LEFT) {
			rect_width = lRect.width;
			rect_height = lRect.height;
			rect_x = lRect.x;
			rect_y = lRect.y;
		} else {
			rect_width = rRect.width;
			rect_height = rRect.height;
			rect_x = rRect.x;
			rect_y = rRect.y;
		} 
		matrixX = rect_width / 5;
		matrixY = rect_height / 5;

		// if (rect_width % 5 != 0) matrixX++;
		// if (rect_height % 5 != 0) matrixY++;
		double	factor = 1.0 / analysisAccuracy;

		for (j = 0; j < matrixY; j++) {
			for (i = 0; i < matrixX; i++) {

				// pix_ref fuellen
				counterX = 0;
				for (sx = (-2); sx < 7; sx++) {
					counterY = 0;
					for (sy = (-2); sy < 7; sy++) {

						// 
						dx = ((double) (rect_x + (i * 5))) + (((double) (sx)) * factor);
						dy = ((double) (rect_y + (j * 5))) + (((double) (sy)) * factor);
						pix_ref[(counterY * 9) + counterX] = bi_ref.getPixel(dx, dy);
						counterY++;
					} 
					counterX++;
				} 

				// vergleichen
				for (x = 0; x < 4; x++) {
					for (y = 0; y < 4; y++) {
						sx = (i * 5) + x;
						sy = (j * 5) + y;
					} 
				} 
			} 
		} 


	} 


	/**
	 * Bereitet die Eingabebilder fuer die Kepstrum-Filterung vor.
	 * noch nicht fertig implementiert. Wird z.Zt. nicht aufgerufen.
	 * eine Korrektur durch, falls notwendig.
	 * @param half    die Bildhaelfte
	 * @param ip_ref  der Imageprocessor des Referenzbildes
	 * @param slice_ref   die Referenz-Bildnummer der Sequenz
	 * @param ip_act  der Imageprocessor des Folgebildes
	 * @param slice_act   die Folge-Bildnummer der Sequenz
	 */
	private void prepareImages(int half, ByteProcessor ip_ref, int slice_ref, ByteProcessor ip_act, int slice_act) {
		if (verbose) {
			IJ.write("Preparing Images for Cepstrum Filter. Accuracy is " + analysisAccuracy);

			// imageStack.setSlice(slice);
		} 
		short pix[][] = new short[2][];

		pix[0] = (short[]) imageStack.getPixels(slice_ref);
		pix[1] = (short[]) imageStack.getPixels(slice_act);
		int			x, y, rect_width, rect_height, rect_x, rect_y;
		double	dx, dy;
		short		val_ref, val_act;

		if (half == JMConst.LEFT) {
			rect_width = lRect.width;
			rect_height = lRect.height;
			rect_x = lRect.x;
			rect_y = lRect.y;
		} else {
			rect_width = rRect.width;
			rect_height = rRect.height;
			rect_x = rRect.x;
			rect_y = rRect.y;
		} 

		if (!(analysisAccuracy > 1.0)) {
			for (x = 0; x < rect_width; x++) {
				for (y = 0; y < rect_height; y++) {
					val_ref = pix[0][((y + rect_y) * width) + x + rect_x];
					val_act = pix[1][((y + rect_y) * width) + x + rect_x];
					ip_ref.putPixel(x, y, val_ref);
					ip_act.putPixel(x, y, val_act);

					// todo: besser nur mit Arrays arbeiten, geht schneller!
				} 
			} 
		} else	// upscale the images with factor given by analysisAccuracy
		 {

			// rect_width = rect_width / 2;
			// rect_height = rect_height / 2;
			if (rect_x + rect_width >= width) {
				rect_width = width - rect_x - 1;
			} 
			if (rect_y + rect_height >= height) {
				rect_height = height - rect_y - 1;
			} 
			double								adder = 1.0 / analysisAccuracy;

			// IJ.write("So, nun geht's los: x/y/w/h: " + rect_x + " / " + rect_y + " / " + rect_width + " / " + rect_height);
			// IJ.write("Und die ImageProcessor Werte " + rect_x + " / " + rect_y + " / " + ip_ref.getWidth() + " / " + ip_ref.getHeight() );
			double								px, py;
			BicubicInterpolation	bi_ref = new BicubicInterpolation(pix[0], width, height);
			BicubicInterpolation	bi_act = new BicubicInterpolation(pix[1], width, height);

			x = 0;
			for (dx = 0.0; dx < (double) rect_width; dx = dx + adder) {
				y = 0;
				for (dy = 0.0; dy < (double) rect_height; dy = dy + adder) {
					px = dx + (double) rect_x;
					py = dy + (double) rect_y;
					val_ref = (short) bi_ref.getPixel(px, py);
					val_act = (short) bi_act.getPixel(px, py);
					ip_ref.putPixel(x, y, val_ref);
					ip_act.putPixel(x, y, val_act);
					y++;
				} 
				x++;
			} 
		} 

		/*
		 * ImagePlus res1 = new ImagePlus("refImage",ip_ref);
		 * IJ.wait(350);
		 * res1.show();
		 * IJ.wait(350);
		 * ImagePlus res = new ImagePlus("actImage",ip_act);
		 * IJ.wait(350);
		 * res.show();
		 * IJ.wait(350);
		 */
	} 


	/**
	 * Berechnet die lineare Bewegungsfilterung.
	 * Hier werden entlang der protokolllierten Verschiebungachse alle Differenzbilder
	 * berechnet und das mit dem minimalen mittl. Grauwert wird als das am besten
	 * geeignetste Zurueckgeliefert.
	 * @param ip_ref  der Imageprocessor des Referenzbildes
	 * @param ip_act  der Imageprocessor des Folgebildes
	 * @param shift   hier wird die Berechnete Verschiebung gespeichert
	 * @return int Zeigt an , ob die Filterung erfolgreich durchgefuehrt werden konnte. 1 = ok, -1 = Fehler
	 */
	int linearTranslationCorrection(ByteProcessor ip_ref, ByteProcessor ip_act, int half, java.awt.Point shift) {
		int							result = -1;
		int							xa, ya, xb, yb;
		int							vx, vy, fa, fxy, dx, dy, fx, fy, x, y;
		int							offsetx = 0;
		int							offsety = 0;
		int							counter = 0;
		java.awt.Point	actPoint = new Point(0, 0);
		java.awt.Point	bestPoint = new Point(0, 0);
		double					minDiffVal = Double.MAX_VALUE;
		double					actDiffVal = 0.0;

		if (half == JMConst.LEFT) {

			// if (min_y_left < 0) offsety = Math.abs(min_y_left) + 1;
			// int min = Math.min(min_x_left,max_x_left);
			// if (min < 0) offsetx = Math.abs(min) + 1;
			xa = min_x_left + offsetx;
			ya = min_y_left + offsety;
			xb = max_x_left + offsetx;
			yb = max_y_left + offsety;

		} else {

			// if (min_y_right < 0) offsety = Math.abs(min_y_right) + 1;
			// int min = Math.min(min_x_right,max_x_right);
			// if (min < 0) offsetx = Math.abs(min) + 1;
			xa = min_x_right + offsetx;
			ya = min_y_right + offsety;
			xb = max_x_right + offsetx;
			yb = max_y_right + offsety;
		} 

		// Linienalgo nach A. Trzewik

		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;
		x = xa;
		y = ya;

		// ImageProcessor ip = newStack.getProcessor();

		while ((x != xb) || (y != yb)) {

			// IJ.write("Punkt auf Linie:" + (x-offsetx) + "/" + (y-offsety) );// Punkt: x,y !!! m_ImageOut->SetPixel(x, y, TRUE);
			// ip.putPixel(x-offsetx,y-offsety,254);
			if (((counter++) % 2) == 0) {

				// Punkt untersuchen
				// IJ.write("Jetzt");
				actPoint.x = x;
				actPoint.y = y;
				actDiffVal = findMinDiffValue(ip_ref, ip_act, actPoint);
				if (actDiffVal < minDiffVal) {
					minDiffVal = actDiffVal;
					bestPoint.x = actPoint.x;
					bestPoint.y = actPoint.y;
				} 
			} 
			fx = fa + vx;
			fy = fa + vy;
			fxy = fx + fy - fa;
			if (Math.abs(fx) <= Math.abs(fy) && Math.abs(fx) <= Math.abs(fxy)) {
				x += dx;
				fa = fx;
			} else if (Math.abs(fy) < Math.abs(fx) && Math.abs(fy) < Math.abs(fxy)) {
				fa = fy;
				y += dy;
			} else {
				fa = fxy;
				x += dx;
				y += dy;
			} 
		} 

		// test also the last Point on the Line
		actPoint.x = x;
		actPoint.y = y;
		actDiffVal = findMinDiffValue(ip_ref, ip_act, actPoint);
		if (actDiffVal < minDiffVal) {
			minDiffVal = actDiffVal;
			bestPoint.x = actPoint.x;
			bestPoint.y = actPoint.y;
		} 
		shift.x = bestPoint.x;
		shift.y = bestPoint.y;

		return result;
	} 


	/**
	 * Hilfsfunktion zur Ueberpruefung der Kepstrum-Ergebnisses.
	 * Ruft checkCepstrumResult af.
	 * @param ip_ref  der Imageprocessor des Referenzbildes
	 * @param ip_act  der Imageprocessor des Folgebildes
	 * @param shift   die von der Kepstrum-Filterung berechnete Verschiebung
	 * @return double den mittleren Grauwert des Differenzbildes an der besten Position
	 */
	private double findMinDiffValue(ByteProcessor ip_ref, ByteProcessor ip_act, java.awt.Point shift) {
		int r = (int) Math.round(analysisAccuracy * 2.0);

		if (r > 7) {
			r = 7;
		} 
		return checkCepstrumResult(ip_ref, ip_act, r, 2, shift);
	} 


	/**
	 * Fuehrt die Sobel-Filterung des eingegebenen Bildes durch.
	 * @param shift   ip der ImageProcessor des Eingabebildes.
	 * @return ByteProcessor das Sobel-Bild.
	 */
	public ByteProcessor filterSobel(ByteProcessor ip) {
		byte[]	pixel = (byte[]) ip.getPixels();
		int			width = ip.getWidth();
		int			height = ip.getHeight();
		byte[]	sobel = new byte[width * height];

		int			x, y;
		int			p1 = 0, p2 = 0, p3 = 0, p4 = 0, p5 = 0, p6 = 0, p7 = 0, p8 = 0, p9 = 0;
		int			sobelResult = 0;
		int			sobelX = 0;
		int			sobelY = 0;

		for (y = 1; y < height - 1; y++) {
			for (x = 1; x < width - 1; x++) {
				if (x > 0 /* first_row */) {
					p1 = pixel[((y - 1) * width) + (x - 1)];	// img.get(xx-1,yy-1);
					p2 = pixel[((y - 1) * width) + (x)];			// img.get(xx  ,yy-1);
					p3 = pixel[((y - 1) * width) + (x + 1)];	// img.get(xx+1,yy-1);
					p4 = pixel[((y) * width) + (x - 1)];			// img.get(xx-1,yy  );
					p5 = pixel[((y) * width) + (x)];					// img.get(xx  ,yy  );
					p6 = pixel[((y) * width) + (x + 1)];			// img.get(xx+1,yy  );
					p7 = pixel[((y + 1) * width) + (x - 1)];	// img.get(xx-1,yy+1);
					p8 = pixel[((y + 1) * width) + (x)];			// img.get(xx  ,yy+1);
					p9 = pixel[((y + 1) * width) + (x + 1)];	// img.get(xx+1,yy+1);
				} else {
					p1 = p2;
					p2 = p3;
					p3 = pixel[((y - 0) * width) + (x - 0)];	// img.get(xx+1,yy-1);
					p4 = p5;
					p5 = p6;
					p6 = pixel[((y - 0) * width) + (x - 0)];	// img.get(xx+1,yy  );
					p7 = p8;
					p8 = p9;
					p9 = pixel[((y - 0) * width) + (x - 0)];	// img.get(xx+1,yy+1);
				} 

				// Calculate horizontal Sobel
				sobelX = p7 + 2 * p8 + p9 - p1 - 2 * p2 - p3;

				// Calculate vertical Sobel
				sobelY = p3 + 2 * p6 + p9 - p1 - 2 * p4 - p7;

				// Correct Overflows
				sobelX = Math.abs(sobelX);
				sobelY = Math.abs(sobelY);

				// Put to resulting Image
				sobelResult = (sobelX + sobelY) / 2;				// / 2 / 4;
				if (sobelResult > 255) {
					sobelResult = 255;

				} 
				sobel[(y * width) + x] = (byte) sobelResult;
			} 
		} 
    // tha 2002.5.13: In der aktuellen Version von ImageJ (1.26) werden keine 
    // unsigned Images mehr unterstuetzt.
		// boolean unsigned = imagePlus.getProcessor().isUnsigned();
		ByteProcessor result = new ByteProcessor(width, height, sobel, ip.getColorModel());

		// show result in Window:
		// IJ.wait(350);
		// ImagePlus res = new ImagePlus("sobel",result);
		// IJ.wait(350);
		// res.show();
		// IJ.wait(350);
		return result;
	} 


	/**
	 * Fuehrt einen Histogramm-Angleich durch.
	 * @param refIP der ImageProcessor des Referenzbildes
	 * @param aktIP der ImageProcessor des Folgebildes
	 * @return ByteProcessor das durch Anwendung des Histogramm-Angleichs korrigierte FOlgebild
	 */
	private ByteProcessor histogramEquilisation(ImageProcessor refIP, ImageProcessor aktIP) {
		if (refIP == null && aktIP == null) {
			return null;

			// if (refIP != null && aktIP == null) return refIP;
			// if (refIP == null && aktIP != null) return aktIP;

			// fuer Testzwecke: Anzeige der Zweischenergebnisse

			/*
			 * ImagePlus res = new ImagePlus("non-adjusted",aktIP);
			 * IJ.wait(350);
			 * res.show();
			 * IJ.wait(350);
			 */

		} 
		byte[]					orgPixel = (byte[]) aktIP.getPixels();
		int							width = refIP.getWidth();
		int							height = refIP.getHeight();
		int[]						aktPixel = new int[width * height];

		ImagePlus				refImp = new ImagePlus("", refIP);
		ImagePlus				aktImp = new ImagePlus("", aktIP);

		ImageStatistics refStat = refImp.getStatistics();
		ImageStatistics aktStat = aktImp.getStatistics();

		double					refMin = refStat.min;
		double					refMax = refStat.max;

		double					aktMin = aktStat.min;
		double					aktMax = aktStat.max;
		double					c0 = refMin - aktMin;

		// IJ.write("HA-> " + refMin + " ! " + refMax + " ! "+ aktMin + " ! "+ aktMax + " ! " + c0);

		aktMin = aktMin + c0;
		aktMax = aktMax + c0;
		int i;

		for (i = 0; i < width * height; i++) {
			aktPixel[i] = ((int) (orgPixel[i])) + (int) c0;
		} 

		// for (i = 0; i < width*height; i++)
		// {
		// aktPixel[i] = aktPixel[i] + c0;
		// }

		for (i = 0; i < width * height; i++) {
			double	val = ((double) (aktPixel[i]));
			double	d = (val - aktMin) / (aktMax - aktMin);
			double	e = refMin + (d * (refMax - refMin));

			aktPixel[i] = ((int) e);
		} 
		for (i = 0; i < width * height; i++) {
			orgPixel[i] = ((byte) (aktPixel[i]));
		} 

		// boolean unsigned = imagePlus.getProcessor().isUnsigned();
		ByteProcessor result = new ByteProcessor(width, height, orgPixel, refIP.getColorModel());

		// fuer Testzwecke: Anzeige der Zweischenergebnisse

		/*
		 * IJ.wait(350);
		 * ImagePlus res2 = new ImagePlus("adjusted",result);
		 * IJ.wait(350);
		 * res2.show();
		 * IJ.wait(350);
		 * ImagePlus ref = new ImagePlus("origin",refIP);
		 * IJ.wait(350);
		 * ref.show();
		 * IJ.wait(350);
		 */
		return result;
	} 

}


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

