Homec4science
Diffusion PTS Analyzer (master) Browse

ijp-pts/src/ch/epfl/biop/pts/PTSAnalyzer.java
c1442a86c48fmaster

Last Change
Raw File

PTSAnalyzer.java
View Options

package ch.epfl.biop.pts;
import java.awt.geom.Point2D;
import java.util.ArrayList;
import java.util.Arrays;
import ij.IJ;
import ij.ImagePlus;
import ij.gui.PolygonRoi;
import ij.gui.Roi;
import ij.measure.Calibration;
import ij.plugin.RoiScaler;
import ij.plugin.Straightener;
import ij.plugin.filter.MaximumFinder;
import ij.process.FloatPolygon;
import ij.process.ImageProcessor;
public class PTSAnalyzer {
private static double scale = 2;
/**
* Main method for analyzing BIOP Pizza Target Slide Images.
* @param ip ImageProcessor of the PTS image
* @param settings PTSSettings object with all relevant settings
* @param criteria PTSCriteria object with convergence conditions
* @return a ResultTable with a single Row, use {@link ij.measure.ResultsTable#addResults()} to append to the Results ResultsTable
*/
public static PTSResult run(ImagePlus imp, PTSSettings settings, PTSCriteria criteria) {
Roi roi = null;
if(imp.getRoi() != null) {
roi = imp.getRoi();
imp.killRoi();
}
ImagePlus imp2 = imp.duplicate();
// If needs scaling do it here
if (settings.isOversample()) {
//IJ.log("Calibration before: "+imp.getCalibration().pixelWidth);
imp2.setProcessor(imp2.getProcessor().resize((int) (imp2.getWidth()*scale)));
Calibration cal = imp2.getCalibration().copy();
cal.pixelHeight /= scale;
cal.pixelWidth /= scale;
imp2.setCalibration(cal);
//IJ.log("Calibration after: "+imp2.getCalibration().pixelWidth);
}
// Find center
Point2D center = new Point2D.Double();
if(roi != null) {
center = new Point2D.Double(roi.getBounds().getCenterX(),roi.getBounds().getCenterY());
//center = settings.getCenterDetector().optimizeCenter(imp.getProcessor(), center);
if(settings.isOversample()) {
center.setLocation(center.getX()*scale, center.getY()*scale);
}
} else {
center = settings.getCenterDetector().computeCenter(imp2.getProcessor());
}
double rad;
boolean is_done = false;
// Do some smart assumption about where to start and stop looking for values;
double min_distance = getMinLookupDist();
// Go to the edge of the image, as defined by the center
double max_distance = getMaxLookupDist(imp2.getProcessor(), center);
rad = (max_distance+min_distance)/2;
PTSResult res = new PTSResult();
double new_rad = Double.MAX_VALUE;
int n_cycles=0;
while(!is_done) {
n_cycles++;
if(n_cycles >100) break;
// Iterate to find the thing
// Define minimum and maximum search diameters
// Make a circle
Roi circ = cleanCircle(center, rad);
imp2.setRoi(circ);
//IJ.log(""+max_distance+" "+min_distance);
// Straighten, until here nothing is calibrated
ImageProcessor band = new Straightener().straighten(imp2, circ, (int) settings.getLineWidth());
res = analyzePeaks(band, settings, rad, imp2.getCalibration());
// Check criteria and continue
if( criteria.checkCriteria(res) ) {
max_distance = rad;
new_rad = ( rad + min_distance ) / 2;
if (Math.abs(rad - new_rad) < criteria.getTolerance()) {
is_done=true;
} else {
rad = new_rad;
}
} else {
min_distance = rad;
rad = ( rad + max_distance ) / 2;
}
}
Roi circ = cleanCircle(center, rad);
if( settings.isOversample()) {
circ = RoiScaler.scale(circ, 1/scale, 1/scale, false);
center.setLocation(center.getX()/scale, center.getY()/scale);
}
res.setCenter(center);
res.setConverged(is_done);
imp.setRoi(circ );
imp2.close();
return res;
}
private static PTSResult analyzePeaks(ImageProcessor band, PTSSettings settings, double current_radius, Calibration cal) {
double[] values = divideArray(sumPixelRows(band), settings.getLineWidth());
// Find Maxima and Minima: These are the POSITIONS
int[] maxima = MaximumFinder.findMaxima(values, 2, true);
// These are SORTED by amplitude, we must sort them again
Arrays.sort(maxima);
int[] minima = MaximumFinder.findMinima(values, 2, true);
Arrays.sort(minima);
// Apply criteria for good and bad detections
double goodDist = current_radius*Math.sin(2*Math.PI/360)*cal.pixelWidth;
double minDist = goodDist*(1-settings.getP2pVarPct()/100);
double maxDist = goodDist*(1+settings.getP2pVarPct()/100);
int min_good_cnt = 0;
int max_good_cnt = 0;
double min_good_avg= 0;
double max_good_avg= 0;
ArrayList<Double> min_d = new ArrayList<Double>();
ArrayList<Double> max_d = new ArrayList<Double>();
// Check peak to Peak Distance Maxima
for(int i=1; i<maxima.length;i++) {
double d = (maxima[i]-maxima[i-1])*cal.pixelWidth;
if( d >= minDist && d <= maxDist) {
max_good_cnt++;
max_good_avg+=values[maxima[i-1]];
max_d.add(d);
}
}
// Check peak to Peak Distance Minima
for(int i=1; i<minima.length;i++) {
double d = (minima[i]-minima[i-1])*cal.pixelWidth;
if( d >= minDist && d <= maxDist) {
min_good_cnt++;
min_good_avg+=values[minima[i-1]];
min_d.add(d);
}
}
max_good_avg /= maxima.length;
min_good_avg /= minima.length;
// Check Michelson Contrast
double michelson = (max_good_avg-min_good_avg)/(max_good_avg+min_good_avg);
// Peak to Peak Statistics
// Need mean and standard deviation
max_d.addAll(min_d); // Use all distances from minima and maxima
double p2p_mean = getMean(toDoubleArray(max_d));
double p2p_stdev = getStdDev(toDoubleArray(max_d));
// Expected Peak 2 Peak
double ep_p2p = current_radius*Math.sin(2*Math.PI/360)*cal.pixelWidth;
// Expected FWHM based on distance
double ep_fwhm = 2*ep_p2p/3;
// Calculated FWHM
double fwhm = 2*p2p_mean/3;
// Total Peaks
double tot_peaks = (double) (maxima.length+minima.length)/2.0;
// Good Peaks
double good_peaks = (double) (min_good_cnt+max_good_cnt)/2.0;
PTSResult res = new PTSResult(ep_p2p, p2p_mean, p2p_stdev, ep_fwhm, fwhm, tot_peaks, good_peaks, good_peaks/tot_peaks * 100, michelson, cal.pixelWidth, current_radius*cal.pixelWidth, new Point2D.Double());
return res;
}
private static double[] divideArray(double[] values, double factor) {
for(int i=0; i<values.length; i++)
values[i] /= factor;
return values;
}
private static double[] sumPixelRows(ImageProcessor band) {
double[] sum = new double[band.getWidth()];
for(int i=0; i<band.getWidth(); i++) {
for(int j=0; j<band.getHeight(); j++) {
sum[i]+= band.getf(i, j);
}
}
return sum;
}
private static Roi cleanCircle(Point2D center, double radius) {
int n_steps = 360*2;
double increment = 2*Math.PI/n_steps;
float[] px = new float[n_steps+1];
float[] py = new float[n_steps+1];
for(int i=0; i<=n_steps; i++) {
px[i] = (float) (radius*Math.cos(i*increment)+center.getX());
py[i] = (float) (radius*Math.sin(i*increment)+center.getY());
}
return new PolygonRoi(new FloatPolygon( px, py), Roi.POLYLINE);
}
private static double getMinLookupDist() {
return 2 / Math.sin(2*Math.PI/360) * 0.80;
}
private static double getMaxLookupDist(ImageProcessor ip, Point2D center) {
// check distances of center to each side of the image
double vertical = Math.min(center.distance(center.getY(), (double)ip.getHeight()),center.distance(center.getY(), 0));
double horizontal = Math.min(center.distance(center.getX(), (double) ip.getWidth()),center.distance(center.getX(), 0));
return Math.min(vertical, horizontal);
}
private static double getMean(double[] nums) {
double total = 0;
for (double value : nums) {
total += value;
}
return total / nums.length;
}
// sample version
private static double getStdDev(double[] nums) {
if (nums.length <= 1) {
return 0;
}
double sum = 0;
double mean = getMean(nums);
for (double value : nums) {
sum += Math.pow(value - mean, 2);
}
return Math.sqrt(sum / (nums.length - 1)); // notice the -1
}
private static double[] toDoubleArray(ArrayList<Double> al) {
double[] target = new double[al.size()];
for (int i = 0; i < target.length; i++) {
target[i] = al.get(i); // java 1.5+ style (outboxing)
}
return target;
}
}
c4science · Help