ImageJ macro to calculate MTF (only with USAF target)

macro_MTF.ijm

requires("1.41j");

function dprint(content){
    if (DEBUG){
        print(content);
    }
}

function average(array){
    x = 0;
    res = 0;
    if (array.length == 0){
        return (0);
    }
    while (x < array.length){
        res += array[x];
        x += 1;
    }
    res /= array.length;
    return (res);
}

function isInArray(array, variable){
    var x = 0;
    while (x < array.length){
        if (variable == array[x]){
            return (true);
        }
        x += 1;
    }
    return (false);
}

function makePlot(pixels, indexLowPeaks, indexHighPeaks, text){
    var x = 0;
    var dots = newArray();
    while (x < pixels.length){
        if (isInArray(indexLowPeaks, x) ||
            isInArray(indexHighPeaks, x)){
            dots = Array.concat(dots, pixels[x]);
        }
        else{
            dots = Array.concat(dots, 255);
        }
        x += 1;
    }
    Plot.create("MTF result", "X", "Y");
    Plot.setLimits(0, pixels.length, 0, 255);
    Plot.setLineWidth(2);
    Plot.setColor("lightGray");
    Plot.add("line", pixels);
    Plot.setColor("red");
    Plot.add("circles", dots);
    Plot.addText(text, 0.8, 0.5);
    setJustification("center");
    Plot.show();
}

function getMtf(startX, startY, width, height, calib, tolerance){
    var size = 0;
    if (width > height){
        size = width;
    }
    else{
        size = height;
    }
    var lowPeaks = newArray();
    var highPeaks = newArray();
    var indexLowPeaks = newArray();
    var indexHighPeaks = newArray();
    var pixels = newArray();
    var x = 0;
    var lastLow = 255;
    var lastHigh = 0;
    while (x < size){
        //y selection
        var cur = 0;
        if (width > height){
            cur = getPixel(startX + x, startY);
        }
        else{
            cur = getPixel(startX, startY + x);
        }
        pixels = Array.concat(pixels, cur);
        //peak algorithm
        if (lastLow == 255){
            if (lastHigh - cur > tolerance || lowPeaks.length == 0){
                dprint("LOW PEAK: " + cur);
                lastLow = cur;
                lastHigh = 0;
                lowPeaks = Array.concat(lowPeaks, cur);
                indexLowPeaks = Array.concat(indexLowPeaks, x);
            }
            else if (highPeaks.length > 0){
                if (cur > highPeaks[highPeaks.length - 1]){
                    highPeaks[highPeaks.length - 1] = cur;
                    indexHighPeaks[indexHighPeaks.length - 1] = x;
                    lastHigh = cur;
                }
            }
        }
        if (lastHigh == 0){
            if (cur - lastLow > tolerance){
                dprint("HIGH PEAK: " + cur);
                lastLow = 255;
                lastHigh = cur;
                highPeaks = Array.concat(highPeaks, cur);
                indexHighPeaks = Array.concat(indexHighPeaks, x);
            }
            else if (lowPeaks.length > 0){
                if (cur < lowPeaks[lowPeaks.length - 1]){
                    lowPeaks[lowPeaks.length - 1] = cur;
                    indexLowPeaks[indexLowPeaks.length - 1] = x;
                    lastLow = cur;
                }
            }
        }
        x += 1;
        dprint("n :" + cur);
    }
    if (DEBUG){
        Array.print(lowPeaks);
        Array.print(highPeaks);
    }
    
    //exclude border low peaks
    if (lowPeaks.length == 4){
        var min = (lowPeaks[1] + lowPeaks[2]) / 2;
    }
    else{
        var min = average(lowPeaks);
    }
    
    var max = average(highPeaks);
    var diff = max - min;
    if (calib){
        File.saveString(diff, CALIB_FILE);
        //dialogMessage("Calib diff = " + diff);
        text = "Calib diff = " + diff + "\nmin = " + min + "\nmax = " + max;
        makePlot(pixels, indexLowPeaks, indexHighPeaks, text);
    }
    else{
        if (File.exists(CALIB_FILE)){
            filestr = File.openAsString(CALIB_FILE);
            cal_diff = parseFloat(filestr);
            
            mtf = diff / cal_diff;
            
            if (lowPeaks.length != 4 || highPeaks.length != 3){
                if (tolerance < MIN_TOLERANCE){
                    mtf = 0;
                    text = "diff = " + diff + "\nMTF= " + mtf;
                    makePlot(pixels, indexLowPeaks, indexHighPeaks, text);
                }
                else{
                    getMtf(startX, startY, width, height, calib, tolerance - SUB_TOLERANCE);
                }
            }
            else{
                var text = "diff = " + diff + "\nMTF= " + mtf;
                makePlot(pixels, indexLowPeaks, indexHighPeaks, text);
            }
        }
        else{
            dialogMessage("Not calibrated!");
        }
    }
    //calcMtf(lowPeaks, highPeaks);
}

function askForCalib(){
    Dialog.create("MTF calc");
    if (File.exists(CALIB_FILE)){
        Dialog.addCheckbox("Calibrate", false);
    }
    else{
        Dialog.addCheckbox("Calibrate", true);
    }
    Dialog.show();
    calib = Dialog.getCheckbox();
    return (calib);
}

function dialogMessage(message){
    Dialog.create("MTF calc");
    Dialog.addMessage(message);
    Dialog.show();
}

function main(){
    CALIB_FILE = "./calib.txt";
    DEBUG = false;
    TOLERANCE = 50;
    MIN_TOLERANCE = 5;
    SUB_TOLERANCE = 2;

    getSelectionBounds(rectX, rectY, width, height);

    var calib = askForCalib();

    getMtf(rectX, rectY, width, height, calib, TOLERANCE);
}


main();

macro "MTF calc [m]" {
    main();
}