// Card_Profile
// (c) GC
// 14/08/2002

// Modifica 27/01/03    Threshold selection (default: 60% of maximum value).
// Modifica 16/12/03    Esclusione del centro della ROI ellittica nella ricerca
//                      dell'"hotspot" (la presenza del centro provocava un bug).

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

/*
 *   Card_Profile  takes the image region bounded by an oval region and samples the oval
 *   at N equal angles around the oval. The program generates a ProfilePlot of either:
 *
 *    0) maximum intensity values along N radii from the oval center to oval points.
 *      This measure is called the circumferential profile and is used in Nuclear Medicine
 *      cardiac blood perfusion studies, in which short axis slices of the left
 *      ventricle are roughly doughnut shaped.
 *    1) Sum of pixel intensities along each radius, from the center to the oval
 *      boundary.
 *    2) pixel intensities at sample points along the oval.

 *   Card_Profile has been developed (as an enhancement of "Oval Profile Plot" by Bill 
 *   O'Connell) to meet a specific request from cardiologists in Nuclear Medicine.
 *
 */

public class Card_Profile  implements PlugInFilter {

	protected ImagePlus imp;

	protected static int npoints = 60;
   protected static int nBHole = 0;
   protected static double scFact = 1.25;
   protected static double percThreshold = 60.0;
	private static String[] modes = {"Maximum Intensity", "Radial Sums",
                                          "Along Oval"};
	public static int CIRCUMFERENTIAL = 0;
	protected static int mode = CIRCUMFERENTIAL;
	protected static boolean hotspots = true;

	public int setup(String args, ImagePlus imp) {
		if(args.equals("about")) {
			showAbout();
			return DONE;
		}
		if (IJ.versionLessThan("1.27e")) {
         IJ.showMessage ("ImageJ 1.27e or more required.");
			return DONE;
      }
		if(imp != null) {
			this.imp = imp;
		}
		return DOES_8G+DOES_16;
	}

	public void run(ImageProcessor ip) {

		Roi roi = imp.getRoi();
		if (roi==null || roi.getType()!=Roi.OVAL) {
			IJ.error("Oval selection required.");
			return;
		}
		if (!getParameters())
			return;;
      Calibration oldCal = imp.getCalibration();
      Calibration cal = oldCal.copy();
      cal.setUnit("Degrees");
      cal.pixelWidth = 360.0 / npoints;
      cal.pixelHeight = 360.0 / npoints;
      imp.setCalibration(cal);
	   CardProfilePlot p = new CardProfilePlot(imp, npoints, 0, hotspots, 
                                              nBHole, scFact, percThreshold);
      p.createWindow();
      imp.setCalibration(oldCal);
	}

	boolean getParameters() {

		GenericDialog param = new GenericDialog("Circumferential Profile", IJ.getInstance());
		param.addNumericField("Number of Points", npoints, 0);
      param.addNumericField("Threshold (% max value)", percThreshold, 0);
      param.addNumericField("Black hole (% mean radius)", nBHole, 0);
      param.addNumericField("Scale factor (mm per pixel)", scFact, 2);
//   	param.addChoice("Analysis mode", modes, modes[mode]);
		param.addCheckbox("Show Hotspots", hotspots);
		param.showDialog();

		if(!param.wasCanceled()) {
			npoints = (int)param.getNextNumber();
         percThreshold = param.getNextNumber();
         if (percThreshold < 10.0 || percThreshold > 100.0)
         {
            IJ.error("Allowed values between 10.0 and 100.0.");
            return false;
         }
         nBHole = (int)param.getNextNumber();
         if (nBHole > 60) {
            IJ.error("Maximum value for black hole: 60.");
            return false;
         }
//			mode = param.getNextChoiceIndex();
         scFact = param.getNextNumber();
         if (scFact < 0.01) {
            IJ.error("Minimum scale factor: 0.01.");
            return false;
         }
			hotspots = param.getNextBoolean();
			return true;

		} else
			return false;

	}
      
	void showAbout() {
      IJ.showMessage(
         "Card_Profile  takes the image region bounded by an oval region and samples the oval" +
         "at N equal angles around the oval. The program generates a ProfilePlot of maximum " +
         "intensity values along N radii from the oval center to oval points."  +
         "This measure is called the circumferential profile and is used in Nuclear Medicine" +
         "cardiac blood perfusion studies, in which short axis slices of the left" +
         "ventricle are roughly doughnut shaped."
      );
	}
}

class CardProfilePlot extends ProfilePlot {

	int npoints = 60;
	int mode = 0;  // Circ. profile
   boolean hotspots = false; // Show graphed points in image
   int nBHole = 0;
   double scFact = 1.25;
   double percThreshold = 60.0;

	ImageProcessor ip2;

	public CardProfilePlot(ImagePlus imp, int npoints, int mode, boolean hotspots, 
                          int nBHole, double scFact, double percThreshold) {

		super();

		this.npoints = npoints;
		this.imp = imp;
		this.mode = mode;
		this.hotspots = hotspots;
      this.nBHole = nBHole;
      this.scFact = scFact;
      this.percThreshold = percThreshold;

		Roi roi = imp.getRoi();
		if (roi==null) {
			IJ.error("Selection required.");
			return;
		}
		int roiType = roi.getType();
		if(roiType != Roi.OVAL) {
			IJ.error("Oval selection required.");
			return;
		}
		// RoiType is Oval
		imp.setColor(Color.white);
		Calibration cal = imp.getCalibration();
		pixelSize = cal.pixelWidth;
		units = cal.getUnits();
		yLabel = cal.getValueUnit();
		ImageProcessor ip = imp.getProcessor();
		if (hotspots) {
			ip2 = ip.duplicate();
			ip2.setColor(Color.white);
			ip2.snapshot();
		}
		ip.setCalibrationTable(cal.getCTable());
		profile = getOvalProfile(roi, ip, imp);

		ip.setCalibrationTable(null);
		ImageWindow win = imp.getWindow();
		if (win!=null)
			magnification = win.getCanvas().getMagnification();
		else
			magnification = 1.0;
	}

   double[] getOvalProfile(Roi roi, ImageProcessor ip, ImagePlus img) {

		Rectangle b = roi.getBoundingRect();
		double width = b.width;
		double height = b.height;

		// get radii from oval center
		double cx = b.x + width / 2.0;
		double cy = b.y + height / 2.0;
		// double theta0 = -Math.PI/2;
		double theta0 = 0;
		double tink = 2 * Math.PI / npoints;
		double[] profile = new double[npoints];
      int[] hotXVect = new int[npoints];
      int[] hotYVect = new int[npoints];

      double dist = 0.0;

      double medRadius = ((width + height) / 2.0) / 2.0;
      double minDist = medRadius * nBHole / 100.0;

      double minRadius = 0.0;
      double maxRadius = 0.0;

      if (width >= height) {
         minRadius = height / 2.0;
         maxRadius = width / 2.0;
      }
      else
      {
         minRadius = width / 2.0;
         maxRadius = height / 2.0;
      }

      IJ.write("CARD PROFILE PARAMETERS");
      IJ.write(" ");

      IJ.write("Image: " + img.getTitle());
      IJ.write(" ");
      IJ.write("Scale factor: " + fmt(scFact) + " mm per pixel");
      IJ.write(" ");
      IJ.write("Mean radius: " + fmt(medRadius * scFact) + " mm (" + fmt(medRadius) + " pixels)");
      IJ.write("Max radius: " + fmt(maxRadius * scFact) + " mm (" + fmt(maxRadius) + " pixels)");
      IJ.write("Min radius: " + fmt(minRadius * scFact) + " mm (" + fmt(minRadius) + " pixels)");
      IJ.write(" ");
      IJ.write("Black hole radius (" + nBHole + "%): " + fmt(minDist * scFact) + " mm (" + fmt(minDist) + " pixels)");
      
      IJ.write(" ");

		int i = 0;

      double minHotval = Double.MAX_VALUE;
      double minTheta = Double.MAX_VALUE;
      double maxHotval = Double.MIN_VALUE;
      double maxTheta = Double.MIN_VALUE;

      int maxPoint = 0;
      int minPoint = 0;

      double hotDist = 0.0;

     	for(double theta = theta0; (theta < theta0 + 2 * Math.PI) && (i < npoints); theta += tink) {

			double dx = Math.cos(theta);
			double dy = Math.sin(theta);

			double x = cx;
			double y = cy;
         double sum = 0;

			double hotval = Double.MIN_VALUE;
			double hotx=0, hoty=0;

			while(roi.contains((int)x, (int)y)) {

				double val = ip.getInterpolatedPixel(x, y);
            dist = Math.sqrt (Math.pow(x-cx, 2) + Math.pow(y-cy, 2));

            if (dist > minDist) {  // era dist >= minDist (16/12/03)
   
               if (val > hotval) { // circumferential

   			      hotval = val;
                  hotDist = dist;
	   		      hotx = x;
				      hoty = y;
               }
            }

				x += dx;
				y += dy;
			}

         if (hotval > maxHotval) {
            maxHotval = hotval;
            maxTheta = theta;
            maxPoint = i;
         }

         if (hotval < minHotval) {
            minHotval = hotval;
            minTheta = theta;
            minPoint = i;
         }

         profile[i] = hotval;
         hotXVect[i] = (int)hotx;
         hotYVect[i] = (int)hoty;

			i++;
		}
            
      if (hotspots) {
         for (i = 0; i < npoints; ++i) {
            ip2.drawDot2(hotXVect[i], hotYVect[i]);
         }       
      }         

      int pointsUnderThreshold = 0;

      double percPoints = 0.0;
      double thValue = maxHotval * percThreshold / 100.0;

      for (i = 0; i < npoints; ++i)
         if (profile[i] < thValue)
            pointsUnderThreshold++;

      percPoints = (double)pointsUnderThreshold * 100.0 / (double)npoints;
         
//    IJ.write (" ");
      IJ.write ("Max value: " + fmt(maxHotval) + " (at angle: " + fmt(maxTheta*360/(2*Math.PI)) + ", step " + maxPoint + ")");
      IJ.write ("Min value: " + fmt(minHotval) + " (at angle: " + fmt(minTheta*360/(2*Math.PI)) + ", step " + minPoint + ")");
      IJ.write ("Min/max ratio: " + fmt(minHotval / maxHotval));
      IJ.write (" ");
      IJ.write ("Threshold " + fmt(percThreshold) + "%: " + fmt(thValue));
      IJ.write ("Points under threshold: " + pointsUnderThreshold + " (" + fmt(percPoints) + "% of " + npoints + " points)");

      IJ.write (" ");
      IJ.write (" ");

		return profile;
	}

	public void createWindow() {
		super.createWindow();
		if (hotspots && ip2!=null) {
			ImagePlus imp2 = new ImagePlus("Hotspots of "+imp.getShortTitle(), ip2);
			imp2.setRoi((Roi)imp.getRoi().clone());
			imp2.show();
		}
	}

   public double fmt(double n) {
      return (double)((int)((n + 0.005) * 100.0) / 100.0);
   }
}

// *** EOF ***