import ij.*;
import ij.measure.*;
import ij.plugin.*;
import ij.plugin.filter.*;
import ij.process.*;
import ij.gui.*;
import java.awt.*;

/**
* Bob Dougherty 6/27/2004
* Adapt Wayne Rasband's Custom Particle Analyzer to the case of Measure Roi

* Version 1 2/26/2007  Update with using the new version of the Custom Particle Analyzer by
*   Greg Joss and Wayne Rasband

*/
/*	License:
	Copyright (c) 2004, 2005, OptiNav, Inc.
	All rights reserved.

	Redistribution and use in source and binary forms, with or without
	modification, are permitted provided that the following conditions
	are met:

		Redistributions of source code must retain the above copyright
	notice, this list of conditions and the following disclaimer.
		Redistributions in binary form must reproduce the above copyright
	notice, this list of conditions and the following disclaimer in the
	documentation and/or other materials provided with the distribution.
		Neither the name of OptiNav, Inc. nor the names of its contributors
	may be used to endorse or promote products derived from this software
	without specific prior written permission.

	THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
	"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
	LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
	A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
	CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
	EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
	PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
	PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
	LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
	NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
	SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

public class Measure_Roi_Curve implements PlugIn {

	public void run(String arg) {
		if (IJ.versionLessThan("1.26i"))
			return;
		ImagePlus imp = WindowManager.getCurrentImage();
		analyzeStackParticles(imp);
	}

	public void analyzeStackParticles(ImagePlus imp) {
		if (imp.getBitDepth()==24)
			{IJ.error("Grayscale image required"); return;}
		MeasureCurveAnalyzer pa = new MeasureCurveAnalyzer();
		int flags = pa.setup("", imp);
		if (flags==PlugInFilter.DONE)
			return;
		if ((flags&PlugInFilter.DOES_STACKS)!=0) {
			for (int i=1; i<=imp.getStackSize(); i++) {
				imp.setSlice(i);
				pa.run(imp.getProcessor());
			}
		} else
			pa.run(imp.getProcessor());
	}

}

class MeasureCurveAnalyzer extends ParticleAnalyzer {
	int NCURVE = 128;
	double[] xL = new double[NCURVE];
	double[] yL = new double[NCURVE];
	double[] xR = new double[NCURVE];
	double[] yR = new double[NCURVE];
	double[] xC = new double[NCURVE];
	double[] yC = new double[NCURVE];
	int sliceC, feretC, orthoC;
	int w,h,nBorder;
	double wp,hp;
	byte[] m;
	double[] x,y;

	// Overrides method with the same in AnalyzeParticles that's called once for each particle
	protected void saveResults(ImageStatistics stats, Roi r) {
 		Calibration cal = imp.getCalibration();
		wp = cal.pixelWidth;
		hp = cal.pixelHeight;
		double feret = 0;
		double ortho = 0;
		if(r instanceof PolygonRoi){
			PolygonRoi pr = (PolygonRoi)r;
			nBorder = pr.getNCoordinates();
			int[] intCoords = pr.getXCoordinates();
			x = new double[nBorder];
			for (int kp = 0; kp < nBorder; kp++)
				x[kp] = wp*intCoords[kp];
			intCoords = pr.getYCoordinates();
			y = new double[nBorder];
			for (int kp = 0; kp < nBorder; kp++)
				y[kp] = hp*intCoords[kp];

			//find the length, or Feret's diameter
			double cos = 1;
			double sin = 0;
			int kTop = 0;
			int kBottom = 0;
			for (int k1 = 0; k1 < nBorder; k1++){
				for (int k2 = k1; k2 < nBorder; k2++){
					double dx = x[k2] - x[k1];
					double dy = y[k2] - y[k1];
					double d = Math.sqrt(dx*dx + dy*dy);
					if(d > feret){
						feret = d;
						cos = dx/d;
						sin = dy/d;
						kTop = k2;
						kBottom = k1;
					}
				}
			}
			//kTop > kBottom by constuction
			//Visualize the ROI as clockwise
			//Pull out left and right side curves, as originally sampled.
			int nLeft = kTop - kBottom + 1;
			double[] xLeft = new double[nLeft];
			double[] yLeft = new double[nLeft];
			for (int k = 0; k < nLeft; k++){
				xLeft[k] = x[kBottom+k];
				yLeft[k] = y[kBottom+k];
			}
			int nRight = (kBottom + nBorder) - kTop + 1;
			double[] xRight = new double[nRight];
			double[] yRight = new double[nRight];
			for (int k = 0; k < nRight; k++){
				xRight[k] = x[(kBottom + nBorder - k)%nBorder];
				yRight[k] = y[(kBottom + nBorder - k)%nBorder];
			}
			//Resample the curves on an equal arc length basis
			resample(xLeft,yLeft,xL,yL);
			resample(xRight,yRight,xR,yR);
			//Define the center curve
			for (int k = 0; k < NCURVE; k++){
				xC[k] = (xL[k] + xR[k])/2;
				yC[k] = (yL[k] + yR[k])/2;
			}
			//Measure the curved length of the centerline and replace the
			//feret length with with result.

			//Option: comment out the next 5 lines to retain the original feret length.
			feret = 0;
			for (int k = 0; k < (NCURVE-1); k++){
				feret += Math.sqrt((xC[k+1]-xC[k])*(xC[k+1]-xC[k]) +
						(yC[k+1]-yC[k])*(yC[k+1]-yC[k]));
			}
			//Find the orthogonal width.  Method: for each k between 0 and NCURVE-1,
			//find the distance between (xL[k],yL[k]) and (xR[k],yR[k]).
			//More generally, with k
			//fixed, attempt to make sure this is a perendicular distance by minimizing
			//as a function of i the distance beween
			//(xL[k+i],yL[k+i]) and (xR[k-i],yR[k-i]).  Then maximixe this result over k.

			ortho = 0;
			for (int k = 0; k < NCURVE; k++){
				double minDist = Double.MAX_VALUE;
				double iDist;
				int iRange = (int)Math.min(k,NCURVE-1-k);
				for(int i = -iRange; i <= iRange; i++){
					iDist = Math.sqrt((xR[k+i]-xL[k-i])*(xR[k+i]-xL[k-i]) +
						(yR[k+i]-yL[k-i])*(yR[k+i]-yL[k-i]));
					if(iDist < minDist){
						minDist = iDist;
					}
				}
				if(minDist > ortho){
					ortho = minDist;
				}
			}
		} else if (r instanceof OvalRoi){
			java.awt.Rectangle rect = r.getBoundingRect();
			double wr = wp*rect.getWidth();
			double hr = hp*rect.getHeight();
			if(wr > hr){
				feret = wr;
				ortho = hr;
			}else{
				feret = hr;
				ortho = wr;
			}

		} else {
			//Roi must be a rectangle
			java.awt.Rectangle rect = r.getBoundingRect();
			double wr = wp*rect.getWidth();
			double hr = hp*rect.getHeight();
			feret = Math.sqrt(wr*wr + hr*hr);
			if (feret == 0){
				ortho = 0;
			}else{
				ortho = 2*wr*hr/feret;
			}
		}

		analyzer.saveResults(stats, r);

		rt.addValue("Curve_Length", feret);
		rt.addValue("Curve_Width", ortho);
		if (showResults)
			analyzer.displayResults();
	}
	//Parameterize a curve as a function of arc length.
	void resample(double[] xIn,double[] yIn,double[] xOut,double[] yOut){
		int nIn = xIn.length;
		int nOut = xOut.length;
		if(nIn == 0){
			for(int i = 0; i < nOut; i++){
				xOut[i] = 0;
				yOut[i] = 0;
			}
			return;
		}
		if(nIn == 1){
			for(int i = 0; i < nOut; i++){
				xOut[i] = xIn[0];
				yOut[i] = yIn[0];
			}
			return;
		}
		double[] arcLength = new double[nIn];
		arcLength[0] = 0;
		for (int i = 1; i < nIn; i++){
			arcLength[i] = arcLength[i-1] +
				Math.sqrt((xIn[i] - xIn[i-1])*(xIn[i] - xIn[i-1]) +
				(yIn[i] - yIn[i-1])*(yIn[i] - yIn[i-1]));
		}
		//Scale so that the total arc length is nOut-1;
		double factor = arcLength[nIn-1]/(nOut - 1);
		for (int i = 0; i < nIn; i++){
			arcLength[i] /= factor;
		}
		arcLength[nIn-1] = nOut - 1;//Avoid any round-off error for this one.
		int iInNext = 1;//The first input index that has arc length >=.iOut
		for (int iOut = 0; iOut < nOut; iOut++){
			//Adjust iInNext to point to the first input that is beyond iOut in arc length.
			while(arcLength[iInNext] < iOut){
				iInNext++;
			}
			//Use iOut as the arc length to perform linear interpolation on
			//the input index.
			double iIn = iInNext;
			double denom = arcLength[iInNext] - arcLength[iInNext-1];
			if (denom > 0){
				iIn = ((iInNext - 1)*(arcLength[iInNext] - iOut) +
						iInNext*(iOut - arcLength[iInNext-1]))/denom;
			}
			int iInInt = (int)iIn;
			xOut[iOut] = xIn[iInInt];
			yOut[iOut] = yIn[iInInt];
			if(iInInt < (nIn - 1)){
				double fract = iIn - iInInt;
				xOut[iOut] += fract*(xIn[iInInt+1] - xIn[iInInt]);
				yOut[iOut] += fract*(yIn[iInInt+1] - yIn[iInInt]);
			}
		}
		return;
	}

}