import java.awt.Frame;
import java.util.ArrayList;
import ij.IJ;
import ij.ImageJ;
import ij.ImagePlus;
import ij.gui.ImageWindow;
import ij.gui.Plot;
import ij.gui.PointRoi;
import ij.gui.Roi;
import ij.plugin.PlugIn;
import ij.process.ImageProcessor;
import javax.swing.JOptionPane;
import edu.emory.mathcs.jtransforms.fft.DoubleFFT_1D;
public class Measure_MTF_Wire implements PlugIn {
public static ArrayList<ImagePlus> getAvailableImagePlus(){
ArrayList<ImagePlus> list = new ArrayList<ImagePlus>();
Frame [] frames = ImageJ.getFrames();
for (Frame frame: frames){
if (frame instanceof ImageWindow){
ImageWindow window = (ImageWindow)frame;
if (! window.isClosed()){
list.add(window.getImagePlus());
}
}
}
return list;
}
public static ImagePlus [] getAvailableImagePlusAsArray(){
ArrayList<ImagePlus> list = getAvailableImagePlus();
ImagePlus [] array = new ImagePlus[list.size()];
for(int i=0; i< list.size(); i++){
array[i] = list.get(i);
}
return array;
}
@Override
public String toString() {
return "Measure Droege MTF";
}
@Override
public void run(String arg) {
try {
configure();
evaluate();
} catch (Exception e) {
// TODO Auto-generated catch block
e.printStackTrace();
}
}
double r = 16, start=-15, end=15;
int sizeOfBead = 10;
boolean configured = true;
ImagePlus image;
Roi roi;
/**
* 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;
}
private boolean init = false;
private int vW;
private int vH;
private int vD;
/**
* 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 );
}
/**
* Returns the minimal and the maximal value in a given array
* @param array the array
* @return an array with minimal and maximal value
*/
public static double [] minAndMaxOfArray(double [] array){
double [] revan = new double [2];
revan[0] = Double.MAX_VALUE;
revan[1] = -Double.MAX_VALUE;
for (int i = 0; i < array.length; i++){
if (array[i] < revan[0]) {
revan[0] = array[i];
}
if (array[i] > revan[1]) {
revan[1] = array[i];
}
}
return revan;
}
/**
* Adds a constant to the first array
* @param sum the first array
* @param toAdd the constant to add
*/
public static double [] add(double[] sum, double toAdd) {
for (int i =0; i < sum.length; i++){
sum[i] += toAdd;
}
return sum;
}
/**
* Performs a 1-D convolution of the input array with the kernel array.<BR>
* New array will be only of size <br>
* <pre>
* output.lenght = input.length - (2 * (kernel.length/2));
* </pre>
* (Note that integer arithmetic is used here)<br>
* @param input the array to be convolved
* @param kernel the kernel
* @return the output array.
*/
public static double [] convolve(double [] input, double [] kernel){
int offset = ((kernel.length) / 2);
double [] revan = new double [input.length - (2* offset)];
double weightSum = 0;
for (int j = 0; j < kernel.length; j++){
weightSum += kernel[j];
}
if (weightSum == 0) weightSum = 1;
for (int i = offset; i < input.length-offset;i++){
double sum = 0;
for (int j = -offset; j <= offset; j++){
sum += kernel[offset+j] * input[i+j];
}
sum /= weightSum;
revan [i-offset] = sum;
}
return revan;
}
/**
* Returns the next power of 2 given a certain integer value
*
* Code was partially taken from ij.plugin.FFT.java::pad().
* Thanks for the inspiration!
*
* @param value the input number.
* @return the next power of two.
*/
public static int getNextPowerOfTwo(int value){
int i = 2;
while (i < value) {
i *= 2;
}
return i;
}
public static double [] fft(double [] array, int padding){
DoubleFFT_1D fft = new DoubleFFT_1D(getNextPowerOfTwo(array.length + padding));
double [] test = new double [getNextPowerOfTwo(array.length+padding) * 2];
for (int i = 0; i < array.length; i++){
test[i*2] = array[i];
}
fft.complexForward(test);
return test;
}
public double [] computeMTF(double [] pixels, int padding){
// remove minimum for frequency analysis:
double [] minmax = minAndMaxOfArray(pixels);
add(pixels, -minmax[0]);
double [] kernel = {-1, 0, 1};
double [] edge = convolve(pixels, kernel);
// FFT
double [] fft = fft(edge, padding);
return fft;
}
/**
* Adds one array to the first array
* @param sum the first array
* @param toAdd the array to add
*/
public static void add(double[] sum, double[] toAdd) {
for (int i =0; i < sum.length; i++){
sum[i] += toAdd[i];
}
}
/**
* Divides all entries of array by divident.
* @param array the array
* @param divident the number used for division.
*/
public static double [] divide(double[] array, double divident) {
for (int i =0; i < array.length; i++){
array[i] /= divident;
}
return array;
}
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;
}
/**
* 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 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 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 Object evaluate() {
if (configured) {
if (roi instanceof PointRoi){
PointRoi point = (PointRoi) roi;
double [] sum = null;
double min = 0;
double max = 0;
int px = point.getBounds().x;
int py = point.getBounds().y;
int pz = image.getCurrentSlice()-1;
int smallstep = 1;
for (int i = (int) start; i < end; i++){
double beta = ((i*smallstep) / 180.0) * Math.PI;
// direction vector
double x = r * Math.cos(beta);
double y = r * Math.sin(beta);
double z = 0;
// Interpolate along line
double [] pixels = getPixels(image, px, px+x, py,py+y,pz,pz+z, (int) (2*r));
double [] minmax = minAndMaxOfArray(pixels);
min += minmax[0];
max += minmax[1];
double [] mtf = computeMTF(pixels, pixels.length * 16);
if (sum == null){
sum = mtf;
} else {
add(sum, mtf);
}
}
double steps = Math.abs(end - start);
divide(sum, steps);
min /= steps;
max /= steps;
double range = max - min;
System.out.println("MTF:");
for (int i=0; i<(sum.length/4); i++) {
System.out.println(abs(i, sum));
}
double pixelSize = image.getCalibration().pixelWidth;
double [] measuredFrequencies = computeComplexFrequencies(sum, pixelSize);
double cutOffFrequency = 1.0 / (2*sizeOfBead*pixelSize);
int cutOffIndex = 0;
for (int i = 0; i < measuredFrequencies.length; i++){
if (measuredFrequencies[i]> cutOffFrequency) {
cutOffIndex = i;
break;
}
}
double mtf [] = new double[cutOffIndex];
double xValues [] = new double[cutOffIndex];
double modelmtf [] = computeModelMTF(sum, range,sizeOfBead);
for(int i=0; i< cutOffIndex;i++){
mtf[i] = abs(i, sum)/ abs(i, modelmtf);
xValues[i] = measuredFrequencies[i];
}
Plot plot = createPlot(xValues, mtf, "2D Bead MTF of " + image.getTitle(), "Frequency", "Power");
try{
plot.show();
} catch (Exception e){
}
}
}
return null;
}
/**
* Queries the User for an Integer value using Swing.
* @param message
* @param initialValue
* @return the chosen int
* @throws Exception
*/
public static int queryInt(String message, int initialValue) throws Exception{
String input = JOptionPane.showInputDialog(message, "" + initialValue);
if (input == null) throw new Exception("Selection aborted");
return Integer.parseInt(input);
}
/**
* Queries the User for a Double values using Swing.
* @param message
* @param initialValue
* @return the chosen double
* @throws Exception
*/
public static double queryDouble(String message, double initialValue) throws Exception{
String input = JOptionPane.showInputDialog(message, "" + initialValue);
if (input == null) throw new Exception("Selection aborted");
return Double.parseDouble(input);
}
public void configure() throws Exception {
image = IJ.getImage();
roi = image.getRoi();
if (roi != null){
r = queryDouble("Radius in pixels: ", r);
sizeOfBead = queryInt("Size of Bead in [px]: ", sizeOfBead);
start = queryDouble("start angle (deg): ", start);
end = queryDouble("end angle (deg): ", end);
configured = true;
}
}
}
/*
* Copyright (C) 2010-2014 - Andreas Maier
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/