import ij.*;
import ij.plugin.filter.PlugInFilter;
import ij.process.*;
import ij.gui.*;
import ij.measure.Measurements;
import ij.measure.ResultsTable;
import java.awt.*;
import java.util.*;
import ij.plugin.frame.RoiManager;

/** Bob Dougherty 4/27/2026
 Compute the Generalized Contrast-to-Noise Ratio (gCNR) for two ROIs on a 32-bit image.
 Reference:
 Rodriguez-Molares, A., Rindal, O. M. H., D’hooge, J., Måsøy, S. E., Austeng, A., Bell, M. A. L., & Torp, H. (2020).
 The generalized contrast-to-noise ratio: a formal definition for lesion detectability. IEEE Transactions on
 Ultrasonics, Ferroelectrics, and Frequency Control, 67(4), 745–759. doi.org
 */

public class GCNR_ implements PlugInFilter {

	ImagePlus imp;

	public int setup(String arg, ImagePlus imp) {
		if (IJ.versionLessThan("1.17j"))
			return DONE;
		this.imp = imp;
		return DOES_32;
	}
	public void run(ImageProcessor ip) {
		RoiManager rm = RoiManager.getInstance();
    	Roi[] rois = rm.getRoisAsArray();
		if (rois == null){
			IJ.error("There needs to be exactly two rois in the RoiManager. Use T to add them.");
			return;
		}
		if (rois.length != 2){
			IJ.error("There needs to be exactly two rois in the RoiManager. Use T to add them.");
			return;
		}

  		int nBins = (int)Prefs.get("gcnr.nBins",256); 
  		double dr = Prefs.get("gcnr.dr",60); 
   		boolean showPlot = Prefs.get("gcnr.showPlot",true); 
  		boolean zeroTop = Prefs.get("gcnr.zeroTop",false); 
 		boolean showResults = Prefs.get("gcnr.showResults",true); 
 		String name = Prefs.get("gcnr.name",""); 
		GenericDialog gd = new GenericDialog("gCNR...", IJ.getInstance());
		gd.addStringField("label", name,20);
		gd.addNumericField("Dynamic range for plotting the reference line (0 to skip)", dr, 1);
		gd.addNumericField("nBins (256)", nBins, 0);
		gd.addCheckbox("Show_plot", showPlot);
		gd.addCheckbox("Shift the top bin to 0 on in the plot", zeroTop);
		gd.addCheckbox("Show_value in Results Table", showResults);
		gd.showDialog();
		if (gd.wasCanceled()) {
			return;
		}
		name = gd.getNextString();
		dr = gd.getNextNumber();
		nBins = (int)gd.getNextNumber();
		showPlot = gd.getNextBoolean();
		zeroTop = gd.getNextBoolean();
		showResults = gd.getNextBoolean();
		Prefs.set("gcnr.name",name);
		Prefs.set("gcnr.nBins",nBins);
 		Prefs.set("gcnr.dr",dr);
  		Prefs.set("gcnr.showPlot",showPlot);
  		Prefs.set("gcnr.zeroTop",zeroTop);
  		Prefs.set("gcnr.showResults",showResults);
		
		ImageStatistics[] stats = new ImageStatistics[2];
  		int options = Measurements.AREA + Measurements.MEAN + Measurements.MIN_MAX;
		for(int iRoi = 0; iRoi < 2; iRoi++){
			imp.setRoi(rois[iRoi],false);
  			stats[iRoi] = imp.getStatistics(options, nBins);

		}
		double globalMin = Math.min(stats[0].min,stats[1].min);
		double globalMax = Math.max(stats[0].max,stats[1].max);			
		
		for(int iRoi = 0; iRoi < 2; iRoi++){
			imp.setRoi(rois[iRoi],false);
 			stats[iRoi] = imp.getStatistics(options, nBins, globalMin, globalMax);
		}
		imp.setRoi(null,false);
		long[][] histograms = new long[2][];
		for(int iRoi = 0; iRoi < 2; iRoi++){
			histograms[iRoi] = stats[iRoi].getHistogram();		
		}
        int n = histograms[0].length;
		double[] x = new double[n];
		double[][] p = new double[2][histograms[0].length];
		double[][] pNorm = new double[2][histograms[0].length];
		double binWidth = (globalMax - globalMin)/(nBins - 1);
		for(int i = 0; i < n; i++){
			for(int iRoi = 0; iRoi < 2; iRoi++)p[iRoi][i] = histograms[iRoi][i];
			x[i] = globalMin + binWidth*i;
		}
		double pNormMax = 0;
		for(int iRoi = 0; iRoi < 2; iRoi++){
			double tot = 0;
       		for(int i = 0; i < n; i++)tot += p[iRoi][i];
       		if(tot > 0){
       			for(int i = 0; i < n; i++){
       				pNorm[iRoi][i]= p[iRoi][i]/tot;
       				if(pNorm[iRoi][i] > pNormMax)pNormMax = pNorm[iRoi][i];
       			}
       		}
		}
		double OVL = 0;
		for(int i = 0; i < n; i++){
			OVL += Math.min(pNorm[0][i],pNorm[1][i]);
		}
		double gCNR = 1 - OVL;
		if(showResults){
			ResultsTable rt = ResultsTable.getResultsTable();
			rt.incrementCounter();
			rt.addValue("gCNR", gCNR);
			rt.show("Results");
		}
		if(showPlot){
			if(zeroTop)for(int i = 0; i < n; i++)x[i] -= x[n-1];
			Plot plot = new Plot("gCNR from "+imp.getTitle(), "value", "Probability Density");
			plot.setFontSize(18);
			plot.setColor(Color.BLACK);
			plot.addPoints(x, pNorm[0],Plot.CIRCLE);
			plot.setColor(Color.RED);
			plot.addPoints(x, pNorm[1],Plot.CROSS);
			plot.setLimits(x[0],x[n-1],0,pNormMax);
			plot.setColor(Color.BLACK);
			if(dr > 0){
				double xdr = x[n-1] - dr;
				double[] xL = new double[]{xdr,xdr};
				double[] yL = new double[]{0,0.9*pNormMax};
				plot.setLineWidth(2);
				plot.addPoints(xL, yL,Plot.LINE);	
			}
			String label = "gCNR = ";
			if(name.length() > 0)label = name+" "+label;
			plot.addLabel(0.02,0.08,label+IJ.d2s(gCNR,3));
			plot.show();
		}
	}
}