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).
*/