content: import java.awt.Rectangle;
  
  import edu.emory.mathcs.jtransforms.fft.DoubleFFT_1D;
  import ij.IJ;
  import ij.ImagePlus;
  import ij.gui.Plot;
  import ij.gui.Roi;
  import ij.plugin.PlugIn;
  
  public class Measure_Slanted_Edge implements PlugIn {
  
  	ImagePlus image;
  	Roi roi;
  	
  	/**
  	 * Computes the absolute value of the complex number at position pos in the array
  	 * @param pos the position
  	 * @param array the array which contains the values
  	 * @return the absolute value
  	 */
  	public static double abs (int pos, double[] array){
  		return Math.sqrt(Math.pow(array[pos *2], 2) + Math.pow(array[(2*pos)+1], 2));
  	}
  
  	
  	public static int getNextPowerOfTwo(int value){
  		int i = 2;
  		while (i < value) {
  			i *= 2;
  		}
  		return i;
  	}
  	
  	public static double []  fft(double []  array){
  		DoubleFFT_1D fft = new DoubleFFT_1D(getNextPowerOfTwo(array.length));
  		double [] test = new double [getNextPowerOfTwo(array.length) * 2];
  		for (int i = 0; i < array.length; i++){
  			test[i*2] = array[i];
  		}
  		fft.complexForward(test);
  		return test;
  	}
  	
  	/* Perform linear regression on a set of points
  	 * 
  	 * @param x		Array of x-coordinates of points
  	 * @param y		Array of y-coordinates of points
  	 * @param n		Number of points
  	 * @return		An array containing three values: the slope of the best-fit line, the y-intercept for the best-fit line, and the Pearson correlation
  	 */
  	private double [] regress(double x[] , double y [], int n)
  	{
  		double sumx, sumy, sumxy, sumx2, sumy2;
  		double slope, intercept, corr;
  		double ret[]=new double[3];
  
  		sumx = sumy = sumxy = sumx2 = sumy2 = 0.0;
  		for (int i=0; i<n; ++i){
  			sumx = sumx + x[i];
  			sumy = sumy + y[i];
  			sumxy = sumxy + x[i] * y[i];
  			sumx2 = sumx2 + x[i] * x[i];
  			sumy2 = sumy2 + y[i] * y[i];
  		}
  
  		/* Linear least-squares regression */
  		slope = (n * sumxy - sumx * sumy ) / (n * sumx2 -  (sumx*sumx) );
  		intercept = (sumy - slope * sumx) / n;
  
  		/* sample Pearson correlation */
  		corr = (n * sumxy - sumx * sumy) / Math.sqrt((n * sumx2 - ( sumx*sumx )) * (n * sumy2 - ( sumy*sumy )) );
  
  		ret[0]=slope;
  		ret[1]=intercept;
  		ret[2]=corr;
  		return ret;
  
  	}
  
  	public void evaluate() {
  		int thispix[], hit[];
  		int xsize,ysize,center_x, lt_limit, rt_limit,st_line,end_line;
  		int midpt,tmp,sl_size, os_factor,dim,offset;
  		double midval,M,B;
  		double inflect[], xpos[], sledge[],tempedge[],tdata[],rdata[];
  		double tdata_plot[],rdata_plot[];
  		double fcenter_x,pixsize;
  		double regarr[], slope, intercept;
  
  		Rectangle bounds=roi.getBounds();
  
  		xsize=bounds.width;
  		ysize=bounds.height;
  
  		// 2D array of the pixels of the roi
  		int pixels[][]=new int[ysize][xsize];
  
  		os_factor=(int)IJ.getNumber("Oversample Factor: ", 2);
  		pixsize=IJ.getNumber("Pixel Size [mm]: ",0.5618);
  		dim=(int)IJ.getNumber("Length of FFT array: ", 512);
  
  		// Get the pixels specified by the bounding box
  		for (int i=0; i<ysize;i++)
  		{
  			for(int j=0; j<xsize;j++)
  			{
  				thispix=image.getPixel(j+bounds.x,i+bounds.y);
  				pixels[i][j]=thispix[0];
  			}
  		}
  
  		center_x=xsize/2;
  
  		// Columns of the image to include in the calculation
  		// (Don't start quite at the edges)
  		lt_limit = center_x - xsize/4;
  		rt_limit = center_x + xsize/4;
  
  		// Rows of the image to include in the calculation
  		st_line = 0;
  		end_line = ysize;
  
  		// Array that will contain the "column" of the point of inflection
  		// for each row of the selected roi.  I say "column" in quotes because
  		// it will, in general, not be an integer value.
  		inflect=new double[ysize];
  
  		// Array that will contain the row numbers of the selected roi
  		// (st_line, st_line + 1, st_line + 2, ..., end_line).
  		xpos=new double[ysize];
  
  		// Iterate over each row of the roi
  		for (int i=st_line; i<end_line; i++)
  		{	   
  			double lt_ave, rt_ave;
  			int l;
  			lt_ave = 0.0; rt_ave = 0.0;
  
  			// Find the average intensity of the left and right sides of the row
  			for (int j=lt_limit;j<(lt_limit+20);j++)
  			{
  				lt_ave += pixels[i][j];
  			}
  
  			for (int j=rt_limit-20;j<rt_limit;j++)
  			{
  				rt_ave += pixels[i][j];
  			}
  			lt_ave/=20.;
  			rt_ave/=20.;
  
  
  			if (lt_ave<rt_ave)
  			{
  				// First, find an initial estimate of the midpoint (inflection point) location
  				midval =  (rt_ave - lt_ave)/2 + lt_ave;
  				l=lt_limit;
  				while (pixels[i][l] < midval) l++;
  				midpt = --l;
  
  				// Now, fine-tune that estimate
  				lt_ave =  rt_ave = 0.0;
  				for (int k=midpt-50;k<midpt-30;k++)
  				{
  					lt_ave += pixels[i][k];
  				}
  				for (int k=midpt+30;k<midpt+50;k++)
  				{
  					rt_ave += pixels[i][k];
  				}
  				lt_ave/=20.;
  				rt_ave/=20.;
  				midval =  (rt_ave - lt_ave) /2.+ lt_ave;
  				l=midpt-30;
  				while (pixels[i][l] < midval) l++;
  				midpt = --l;
  			}
  
  			else
  			{
  				// First, find an initial estimate of the midpoint (inflection point) location
  				midval =  (lt_ave - rt_ave)/2 + rt_ave;
  				l=lt_limit;
  				while (pixels[i][l] > midval) l++;
  				midpt = --l;
  
  				// Now, fine-tune that estimate
  				lt_ave = rt_ave = 0.0;
  				for (int k=midpt-50;k<midpt-30;k++)
  				{
  					lt_ave += pixels[i][k];
  				}
  				for (int k=midpt+30;k<midpt+50;k++)
  				{
  					rt_ave += pixels[i][k];
  				}
  				lt_ave/=20.;
  				rt_ave/=20.;
  
  				midval =  (lt_ave - rt_ave)/2. + rt_ave;
  				l=midpt-30;
  				while (pixels[i][l] > midval) l++;
  				midpt = --l;
  			}
  
  			// Now, the idea is to find the "column" of the calculated inflection point with sub-integer accuracy.
  			// (Again, I put "column" in quotes because of the sub-integer accuracy involved).
  			// The calculated inflection point has intensity 'midval' and is assumed to be on a line between
  			// column 'midpt' and 'midpt'+1.
  			M = pixels[i][midpt+1]-pixels[i][midpt];
  			B = pixels[i][midpt] - M*midpt;
  
  			// inflect[i] is the i'th row's "column" index for the midpoint (inflection point)
  			inflect[i] = (midval- B)/M;
  			xpos[i]=i;
  
  		}
  
  		// Perform linear regression on the midpoint location to find the
  		// best-fit estimate of where the edge is within the roi
  		regarr=regress(xpos,inflect,end_line-st_line);
  		slope=regarr[0];
  		intercept=regarr[1];
  
  		// Find the "column" (with sub-integer accuracy) of 'midval' (inflection point) for the center row of the roi
  		fcenter_x = slope*(st_line+(end_line+st_line)/2)+intercept;
  
  		// Prepare to deal with some array index offset business later
  		tmp = Math.max((int)(fcenter_x-xsize/2.),0);
  
  		// Size of the over-sampled edge response
  		sl_size = xsize*os_factor;
  
  		// Scale variables appropriately according to how much we over-sample
  		fcenter_x = fcenter_x*os_factor;
  		pixsize/=os_factor;
  
  		// 'sledge' will contain the over-sampled edge response function
  		sledge=new double[dim];
  		hit=new int[dim];
  
  		for (int j=0;j<dim;j++)
  		{
  			sledge[j] = 0.0;
  			hit[j] = 0;
  		}
  
  		// The idea for the following 'for' loop is to upsample each row
  		// and add a shifted version of it to 'sledge'.
  		// The shift ('delta' below) is how far that specific row's
  		// "column" (again, sub-pixel accuracy) index for the inflection point
  		// is from that of the center row's 
  		for (int i=st_line;i<end_line;i++)
  		{
  			int delta;
  			tempedge=new double[sl_size];
  			for (int j=0;j<sl_size;j++)
  			{
  				tempedge[j] = 0;
  			}
  
  			for (int j=0;j<sl_size;j+=os_factor)
  			{
  				tempedge[j] = pixels[i][j/os_factor+tmp];
  			}
  
  			delta = (int)(inflect[i] *os_factor - fcenter_x);
  			if (delta > 0 )
  				for (int j=0;j<sl_size-delta;j++)
  				{
  					sledge[j] += tempedge[j+delta];
  					if (tempedge[j+delta]!=0)  hit[j]++;
  				}
  			if (delta <= 0)
  				for (int j=sl_size-1;j>=-delta;j--)
  				{
  					sledge[j] += tempedge[j+delta];
  					if (tempedge[j+delta]!=0)  hit[j]++;
  				}
  		}
  
  		for (int j=0;j<dim;j++)
  		{
  			if (hit[j] > 0) sledge[j]/=(double)hit[j];
  		}
  
  
  		// Pad for the array for the FFT
  
  		if (sl_size < dim )
  		{
  			int pad = (dim-sl_size)/2;
  			double lt_ave=0.0, rt_ave = 0.0;
  			for (int k=0;k<20;k++)
  				lt_ave += sledge[k];
  			for (int k=sl_size-20;k<sl_size;k++)
  				rt_ave += sledge[k];
  			lt_ave/=20.;
  			rt_ave/=20.;
  			for (int k=(dim-pad);k<dim;k++)
  				sledge[k] = rt_ave;
  			for (int k=(dim-pad-1);k>pad;k--)
  				sledge[k]= sledge[k-pad];
  			for (int k=0;k<=pad; k++)
  				sledge[k] = lt_ave;
  			sl_size = dim;
  		}
  
  		offset = (sl_size - dim)/2;
  		tdata=new double[dim];
  		rdata=new double[dim+1];
  
  		// Differentiate the ERF to get the LSF
  		for (int j=0;j<dim-1;j++)
  		{
  			tdata[j+1] = sledge[j+1+offset] - sledge[j+offset];
  		}
  		tdata[0]=0;
  
  		// Shift the LSF
  		// (I don't know why this is necessary)
  		for (int i=1;i<dim/2;i++)
  			rdata[i] = tdata[i+dim/2];
  		for (int i=dim/2+1;i<dim+1;i++)
  			rdata[i] = tdata[i-dim/2];
  
  		// 'tdata' and 'rdata' have some funky indexing, so we want to make
  		// versions of them that are plot-able (i.e., have good indexing)
  		tdata_plot=new double[dim-1];
  		rdata_plot=new double[dim];
  
  		for(int i=0;i<dim-1;i++) tdata_plot[i]=tdata[i+1];
  		for(int i=0;i<dim;i++) rdata_plot[i]=rdata[i+1];
  
  		// Plot ERF and LSF
  		createPlot("ERF", sledge).show();
  		createPlot("LSF", tdata_plot).show();
  
  		double MTF1[]= fft(rdata_plot);				   
  		// This is from VisualizationUtil.java
  		// I should really modify/add methods in that file
  
  		// Calculate magnitude of MTF
  		double [] absValues = new double [MTF1.length / 2];
  		for (int i = 0; i < absValues.length; i ++){
  			absValues[i] = abs(i, MTF1);
  		}
  
  		// Normalize MTF to have max value of 1
  		for(int i = absValues.length-1; i >= 0; i--){
  			absValues[i]=absValues[i]/absValues[0];
  		}
  		double yValues[]=absValues;
  
  		double xValues[]=new double[yValues.length];
  
  		// Create the frequency axis, scaled appropriately
  		for (int i=0; i < xValues.length; i++) xValues[i]=i*(1.0/xValues.length)*(1.0/pixsize);
  
  		// Plot MTF
  		createPlot(xValues, yValues, "MTF", "Frequency [mm^(-1)]", "").show();
  	}
  	
  	public static Plot createPlot(double [] xValues, double [] yValues, String title, String xLabel, String yLabel){
  		double miny = Double.MAX_VALUE;
  		double maxy = -Double.MAX_VALUE;
  		double minx = Double.MAX_VALUE;
  		double maxx = -Double.MAX_VALUE;
  		for (int i = 0; i < xValues.length; i ++){
  			miny = (yValues[i] < miny) ? yValues[i] : miny;
  			maxy = (yValues[i] > maxy) ? yValues[i] : maxy;
  			minx = (xValues[i] < minx) ? xValues[i] : minx;
  			maxx = (xValues[i] > maxx) ? xValues[i] : maxx;
  
  		}
  		if (miny == maxy){
  			maxy++;
  		}
  		if (minx == maxx){
  			maxx++;
  		}
  		Plot plot = new Plot(title, xLabel, yLabel, xValues, yValues, Plot.DEFAULT_FLAGS);
  		plot.setLimits(minx, maxx, miny, maxy);
  		return plot;
  	}
  	
  	public static Plot createPlot(String title, double [] yValues){
  		return createPlot(yValues, title, "Column", "Power");
  	}
  
  	public static Plot createPlot(double [] yValues, String title, String xLabel, String yLabel){
  		double [] xValues = new double [yValues.length];
  		double min = Double.MAX_VALUE;
  		double max = -Double.MAX_VALUE;
  		for (int i = 0; i < xValues.length; i ++){
  			min = (yValues[i] < min) ? yValues[i] : min;
  			max = (yValues[i] > max) ? yValues[i] : max;
  			xValues[i] = i + 1;
  		}
  		if (min == max){
  			max++;
  		}
  		Plot plot = new Plot(title, xLabel, yLabel, xValues, yValues, Plot.DEFAULT_FLAGS);
  		plot.setLimits(1, xValues.length, min, max);
  		return plot;
  	}
  	
  	public Measure_Slanted_Edge(){
  
  	}
  
  	@Override
  	public void run(String arg) {
  		image = IJ.getImage();
  		roi = image.getRoi();
  		if (roi != null){
  			evaluate();
  		}
  
  	}
  
  }
  
  /*
   * Copyright (C) 2010-2014 - Andreas Maier 
   * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
  */
  