package edu.stanford.rsl.apps.gui.roi; import ij.IJ; import java.awt.Rectangle; import edu.stanford.rsl.conrad.utils.FFTUtil; import edu.stanford.rsl.conrad.utils.VisualizationUtil; /** * Estimate the pre-sampled (pre-pixelated) MTF of an X-ray detector if a slanted edge is used. * * This code is ported from C. To see the original code, login to roentgen.stanford.edu and look at: * /home/group/data_fahrig/fahrigDirectoryFromLucas/mtf/original/mtfslant2.c * Be warned, that file is poorly commented! I have done my best to infer the intention of the code. * * @author Derek Macklin * */ public class MeasureSlantMTF extends EvaluateROI { /* 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; } @Override public Object evaluate() { if (configured) { 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 VisualizationUtil.createPlot("ERF", sledge).show(); VisualizationUtil.createPlot("LSF", tdata_plot).show(); double MTF1[]=FFTUtil.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] = FFTUtil.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 VisualizationUtil.createPlot(xValues, yValues, "MTF", "Frequency [mm^(-1)]", "").show(); } return null; } public void configure() throws Exception { image = IJ.getImage(); roi = image.getRoi(); if (roi != null){ configured = true; } } @Override public String toString() { return "Measure MTF using a slanted edge"; } } /* * Copyright (C) 2010-2014 - Andreas Maier * CONRAD is developed as an Open Source project under the GNU General Public License (GPL). */