package edu.stanford.rsl.apps.gui.roi;
import ij.IJ;
import ij.ImagePlus;
import ij.gui.Plot;
import ij.gui.PointRoi;
import ij.process.ImageProcessor;
import edu.stanford.rsl.conrad.utils.DoubleArrayUtil;
import edu.stanford.rsl.conrad.utils.FFTUtil;
import edu.stanford.rsl.conrad.utils.UserUtil;
import edu.stanford.rsl.conrad.utils.VisualizationUtil;
public class Measure3DBeadMTF extends EvaluateROI {
private boolean init = false;
private int vW;
private int vH;
private int vD;
private double r = 16;
public double [] computeComplexFrequencies(double [] fft, double voxelsize){
double nyquistFrequency = 1 / (2* voxelsize);
double stepsize = nyquistFrequency / (fft.length/4.0);
double [] reval = new double [fft.length/4];
for(int i = 0; i<fft.length/4;i++){
reval[i]= i * stepsize;
}
return reval;
}
public double [] computeModelMTF(double [] fft, double range, int pixelsize){
double [] reval = new double [fft.length/2];
for (int i = 0; i < pixelsize/2; i++)
reval[i] = range;
return computeMTF(reval, 0);
}
public double [] computeMTF(double [] pixels, int padding){
// remove minimum for frequency analysis:
double [] minmax = DoubleArrayUtil.minAndMaxOfArray(pixels);
DoubleArrayUtil.add(pixels, -minmax[0]);
double [] kernel = {-1, 0, 1};
double [] edge = DoubleArrayUtil.convolve(pixels, kernel);
// FFT
double [] fft = FFTUtil.fft(edge, padding);
return fft;
}
/**
* Method to perform trilinear interpolation along a line through an ImagePlus.
* @param image the ImagePlus
* @param x1 start x
* @param x2 end x
* @param y1 start y
* @param y2 end y
* @param z1 start z
* @param z2 end z
* @param numberOfQuantizationSteps
* @return the array with the interpolated values
*/
public double [] getPixels(ImagePlus image, double x1, double x2, double y1, double y2, double z1, double z2, int numberOfQuantizationSteps){
double [] revan = new double[numberOfQuantizationSteps];
// direction
double x = (x2 - x1) / (numberOfQuantizationSteps-1);
double y = (y2 - y1) / (numberOfQuantizationSteps-1);
double z = (z2 - z1) / (numberOfQuantizationSteps-1);
for (int i = 0; i< numberOfQuantizationSteps; i++) {
revan[i] = trilinear(image, x1 + (i*x), y1+ (i*y), z1+ (i*z));
}
return revan;
}
@Override
public Object evaluate() {
if (configured) {
if (roi instanceof PointRoi){
PointRoi point = (PointRoi) roi;
double [] sum = null;
int px = point.getBounds().x;
int py = point.getBounds().y;
int pz = image.getCurrentSlice()-1;
int bigstep = 360; // 60;
int smallstep = 1; //3.0;
for (int i = 0; i < bigstep; i++){
double alpha = ((i*smallstep) / 180.0) * Math.PI;
for (int j = 0; j < bigstep; j++){
double beta = ((i*smallstep) / 180.0) * Math.PI;
// direction vector
double x = r * Math.sin(alpha) * Math.cos(beta);
double y = r * Math.sin(alpha) * Math.sin(beta);
double z = r * Math.cos(alpha);
// Interpolate along line
double [] pixels = getPixels(image, px, px+x, py,py+y,pz,pz+z, (int) (2*r));
// remove mean for frequency analysis:
double [] minmax = DoubleArrayUtil.minAndMaxOfArray(pixels);
DoubleArrayUtil.add(pixels, -minmax[0]);
double [] kernel = {-1, 0, 1};
double [] edge = DoubleArrayUtil.convolve(pixels, kernel);
// FFT
double [] fft = FFTUtil.fft(edge);
if (sum == null){
sum = fft;
} else {
DoubleArrayUtil.add(sum, fft);
}
}
}
DoubleArrayUtil.divide(sum, Math.pow(bigstep, 2));
System.out.println("MTF:");
for (int i=0; i<(sum.length/4); i++) {
System.out.println(FFTUtil.abs(i, sum));
}
Plot plot = VisualizationUtil.createHalfComplexPowerPlot(sum, "Edge MTF of " + image.getTitle());
try{
plot.show();
} catch (Exception e){
}
} else {
throw new RuntimeException("A PointRoi is required to measure the 3D MTF.");
}
}
return null;
}
public void configure() throws Exception {
image = IJ.getImage();
roi = image.getRoi();
if (roi != null){
r = UserUtil.queryDouble("Radius in pixels: ", r);
configured = true;
}
}
@Override
public String toString() {
return "Measure 3-D MTF of a bead";
}
/**
* Method to initialize the trilinear interpolation.
* @param data3D
*/
private void init(ImagePlus data3D){
if (!init){
vW = data3D.getWidth()-1;
vH = data3D.getWidth()-1;
vD = data3D.getStackSize()-1;
init = true;
}
}
/**
* Trilinear Interpolation in an ImagePlus.<BR>
* Adopted from Volume Viewer by Kai Uwe Barthel: barthel (at) fhtw-berlin.de
*
* This method is initialized in the first call with the volume dimensions to save computation time.<br>
* If this interpolation method is used from somewhere else, please use this method only on volumes of the same dimension. Instantiate one Measure3DBeadMTF Object per distinct volume dimension.
*
* @param data3D the ImagePlus
* @param x the x coordinate
* @param y the y coordinate
* @param z the z coordinate
* @return the interpolated value
*/
public double trilinear(ImagePlus data3D, double x, double y, double z) {
init(data3D);
int tx = (int)x;
double dx = x - tx;
int tx1 = (tx < vW) ? tx+1 : tx;
int ty = (int)y;
double dy = y - ty;
int ty1 = (ty < vH) ? ty+1 : ty;
int tz = (int)z;
double dz = z - tz;
int tz1 = (tz < vD) ? tz+1 : tz;
ImageProcessor ptz = data3D.getStack().getProcessor(tz+1);
ImageProcessor ptz1 = data3D.getStack().getProcessor(tz1+1);
float v000 = ptz.getPixelValue(tx, ty);
float v001 = ptz1.getPixelValue(tx, ty);
float v010 = ptz.getPixelValue(tx, ty1);
float v011 = ptz1.getPixelValue(tx, ty1);
float v100 = ptz.getPixelValue(tx1, ty);
float v101 = ptz1.getPixelValue(tx1, ty);
float v110 = ptz.getPixelValue(tx1, ty1);
float v111 = ptz1.getPixelValue(tx1, ty1);
return (
(v100 - v000)*dx +
(v010 - v000)*dy +
(v001 - v000)*dz +
(v110 - v100 - v010 + v000)*dx*dy +
(v011 - v010 - v001 + v000)*dy*dz +
(v101 - v100 - v001 + v000)*dx*dz +
(v111 + v100 + v010 + v001 - v110 - v101 - v011 - v000)*dx*dy*dz + v000 );
}
}
/*
* Copyright (C) 2010-2014 - Andreas Maier
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/