// MomentMacroJ_v1_5.txt,a version of the MomentMacro modified for use with ImageJ.
// Runs best with 8-bit images of white cross-sections on black background.
// Requires ImageJ version 1.34g or later.  See http://rsb.info.nih.gov/ij/upgrade/index.html 
// for the latest version of ImageJ.

var threshl = 0;
var threshu = 255;
var scalar = 1;
var units = "pixels";

var Sn, Sx, Sy, Sxx, Syy, Sxy, rot2, density, Centrex, Centrey, Theta, J, Zp, maxRad;

import java.lang.Math

macro "Moment Calculation" {
    // This macro calculates the anisometry and bulkiness of binary particles
    // based on the dynamically equivalent ellipse as defined by Medalia (1970).
    // Currently set to work with 8-bit files only
    requires("1.34g"); // add prompt for earlier versions

    title = getTitle();
    rename(title + "^");
    title = getTitle();
    selectWindow(title);

    snapshot();
    scale = getBoolean("Reset scale?  (" + scalar + units + "/pixel)");
    if (scale == 1) {
        ResetUnits();
    }

    drawaxes = getBoolean("Draw principal axes?");
    name = getString('Enter name of sample, max 14 chars:', "sample");
    threshl = getNumber('Lower Threshold:  ', threshl); //set threshold limits
    threshu = getNumber('Upper Threshold:  ', threshu);
    getSelectionBounds(xmin, ymin, xmax, ymax); //  ymin is measured from top,xmin from left
                                                //  xmax & ymax= diameters along x & y axes

    // Calculate Moments
    run("Copy");
    newImage('Calculating Moments', '8-bit black', xmax + 2, ymax + 2, 1);
    run("Paste");
    CalcSums();
    if (Sn == 0) {
        showMessage('Threshold limits too narrow');
        selectWindow('Calculating TA');
        close();
        selectWindow('Calculating Moments');
        //close();
        exit;
    }

    Cx = Sx / Sn + xmin - 1; // x-coord. of Centroid
    Cy = Sy / Sn + ymin - 1; // y-coord. of Centroid
    // x-coord of Centoid and y-coord of Centroid held for calculation of max radius

    // Set x and y used to draw axis
    Centrex = Cx;
    Centrey = Cy;

    SxSnhold = Sx / Sn;
    SySnhold = Sy / Sn;
    xminhold = xmin;
    yminhold = ymin;
    // following  code calculates y (dist from neutral axis)
    if (((ymax + ymin) - Cy) > (Cy - ymin)) {
        Ymaxrad = ymax + ymin - Cy;
        Yminrad = Cy - ymin;
        BigY = Ymaxrad;
    } else {
        Ymaxrad = Cy - ymin;
        Yminrad = ymax + ymin - Cy;
        BigY = Ymaxrad;
    }
    // following code calculates x (dist from neutral axis)
    if (((xmax + xmin) - Cx) > (Cx - xmin)) {
        Xmaxrad = xmax + xmin - Cx;
        Xminrad = Cx - xmin;
        BigX = Xmaxrad;
    } else {
        Xmaxrad = Cx - xmin;
        Xminrad = xmax + xmin - Cx;
        BigX = Xmaxrad;
    }
    Xmaxrad = Xmaxrad / scalar; //calibrating radii
    Xminrad = Xminrad / scalar;
    Ymaxrad = Ymaxrad / scalar;
    Yminrad = Yminrad / scalar;
    Parea = Sn;
    Myy = Sxx - (Sx * Sx / Sn);
    Mxx = Syy - (Sy * Sy / Sn);
    Mxy = Sxy - (Sx * Sy / Sn);
    if (Mxy == 0) {
        Theta = 0;
    } else {
        Theta = atan(((Mxx - Myy) + sqrt(sqr(Mxx - Myy) + (4 * sqr(Mxy)))) / (2 * Mxy)) * 180 / 3.141592654;
    }
    CA = Parea / sqr(scalar);
    //end CA
    Cx = Cx / scalar; //save x-coord of centroid
    Cy = Cy / scalar; //save y-coord of centroid
    BigX = BigX / scalar; //calibrating x,y,Cx,Cy
    BigY = BigY / scalar;
    Ix = Mxx / (sqr(sqr(scalar))); //save mom about x-axis
    Iy = Myy / (sqr(sqr(scalar)));
    Zx = Ix / BigY; //save section moduli
    Zy = Iy / BigX;

    selectWindow('Calculating Moments');
    rot2 = Theta * 3.141592654 / 180;
    CalcSums();
    M1 = Sxx - (Sx * Sx / Sn);
    M2 = Syy - (Sy * Sy / Sn);
    R1 = sqrt(M1 / Parea);
    R2 = sqrt(M2 / Parea);
    Imax = M1 / (sqr(sqr(scalar)));
    Imin = M2 / (sqr(sqr(scalar)));
    Rmaks = R1 / scalar;
    Rmyn = R2 / scalar;
    rot2 = 0;
    // Theta= -Theta; REMOVED due to no apparent function and gives incorrect principal axes (5/24/2005)

    selectWindow('Calculating Moments');
    close();

    // calculate TA
    newImage("Calculating TA", "8-bit White", xmax, ymax, 1);
    run("Restore Selection");
    setForegroundColor(1, 1, 1);
    run("Fill");
    run("Measure");
    Area = getResult("Area");
    TA = Area / sqr(scalar);
    close();
    //end TA

    selectWindow("Results");
    run("Close");

    // Draw principal axes
    if (drawaxes == 1) {
        selectWindow(title);
        DrawAxis();
    }

    // Calculate the maximum distance of any pixel from the center
    run("Copy");
    newImage('Calculating Max Radius', '8-bit black', xmax + 2, ymax + 2, 1);
    run("Paste");
    MAXRAD();
    close();

    // Calculate the polar moment of inertia
    J = Ix + Iy;
    //Calculate polar modulus
    Zp = J / maxRad;

    // displays the results in the upper left corner of the text window
    setFont('Monaco', 10);
    tab = fromCharCode(0009);
    if (!isOpen("Log")) {
        print("Name" + tab + "Scale" + tab + "TA" + tab + "CA" + tab + "Xbar" + tab + "Ybar" + tab + "Ix" + tab + "Iy" + tab + "Imax" + tab + "Imin" + tab + "Theta" + tab + "Zx" + tab + "Zy" + tab + "MaxXrad" + tab + "MaxYrad" + tab + "J" + tab + "Zp" + tab + "MaxRad");
    }
    print(name + tab + scalar + "pixels/" + units + tab + TA + tab + CA + tab + Cx + tab + Cy + tab + Ix + tab + Iy + tab + Imax + tab + Imin + tab + Theta + tab + Zx + tab + Zy + tab + Xmaxrad + tab + Ymaxrad + tab + J + tab + Zp + tab + maxRad);
    saveLog();

    selectWindow(title);
    reset();
} // End of main macro "Calculating Moments"

// Begin functions
function CalcSums() {
    Sn = 0;
    Sx = 0;
    Sy = 0;
    Sxx = 0;
    Syy = 0;
    Sxy = 0;
    for (y = 1; y <= ymax; y = y + 1) {
        for (x = 1; x <= xmax; x = x + 1) {
            if ((getPixel(x, y) >= threshl) && (getPixel(x, y) <= threshu)) {
                Sn = Sn + 1;
                Sx = Sx + (x * cos(rot2) + y * sin(rot2));
                Sy = Sy + (y * cos(rot2) - x * sin(rot2));
                Sxx = Sxx + sqr((x * cos(rot2) + y * sin(rot2)));
                Syy = Syy + sqr((y * cos(rot2) - x * sin(rot2)));
                Sxy = Sxy + ((y * cos(rot2) - x * sin(rot2)) * (x * cos(rot2) + y * sin(rot2)));
            }
        }
    }
}

// Find pixel furthest away from the center and calculated the distance to that pixel
function MAXRAD() {
    maxRad = 0;
    for (y = 1; y <= ymax; y = y + 1) {
        for (x = 1; x <= xmax; x = x + 1) {
            if ((getPixel(x, y) >= threshl) && (getPixel(x, y) <= threshu)) {
                maxRad1 = (x - SxSnhold) * (x - SxSnhold);


                maxRad2 = (y - SySnhold) * (y - SySnhold);


                maxRad3 = (sqrt(maxRad2 + maxRad1)) / scalar;
                if (maxRad3 > maxRad) {
                    maxRad = maxRad3;
                }
            }
        }
    }
}

function CalcTA() {
    Sn = 0;
    for (y = 1; y <= ymax; y = y + 1) {
        for (x = 1; x <= xmax; x = x + 1) {
            //          if (getPixel(x,y)<255) {
            Sn = Sn + 1;
            // }
        }
    }
}

function saveLog() {
    selectWindow("Log");
    run("Text..."); // File>Save As>Text
}

function sqr(n) {
    return n * n;
}

function ResetUnits() {
    //optional function to set the units to be used in the results for the Moment Calculation Macro.
    units = getString('Enter alternate units as desired:', units);
    scalar = getNumber('Enter number of pixels per ' + units + ':', scalar);
}

function DrawAxis() {
    //optional function to draw major/minor axis of ellipse
    moveTo(Centrex, Centrey);
    setColor(255, 255, 255);
    lineTo(Centrex - (cos((0 - Theta) * 3.141592654 / 180) * 2 * R1),
        Centrey + (sin((0 - Theta) * 3.141592654 / 180) * 2 * R1));
    moveTo(Centrex, Centrey);
    lineTo(Centrex - (cos((Theta + 90) * 3.141592654 / 180) * 2 * R2),
        Centrey - (sin((Theta + 90) * 3.141592654 / 180) * 2 * R2));
    moveTo(Centrex, Centrey);
    lineTo(Centrex + (cos((0 - Theta) * 3.141592654 / 180) * 2 * R1),
        Centrey - (sin((0 - Theta) * 3.141592654 / 180) * 2 * R1));
    moveTo(Centrex, Centrey);
    lineTo(Centrex + (cos((Theta + 90) * 3.141592654 / 180) * 2 * R2),
        Centrey + (sin((Theta + 90) * 3.141592654 / 180) * 2 * R2));
    setThreshold(255, 255);
}

// Written by: Matthew Warfel  - Cornell University - 4/4/97
// Modified by: Stanley Serafin- Johns Hopkins University - 6/30/00
// Modified and adapted for ImageJ by:  Valerie Burke DeLeon - 2/21/05
// Modified and adapted fro ImageJ by: Adam David Sylvester - 9/28/2012
// Modified by: Guy Taylor
//
// Updates:
// 2005/05/24 - v1.2  added "DrawAxis" function modified from MomentMacro (VBD)
// 2006/03/17 - v1.2  renamed "neutral axes" to more correct term "principal axes" (VBD)
// 2006/07/21 - v1.3  replaced "/*" with "//" character to define initial comment lines, following reports of comments read as code (VBD)
// 2012/09/28 - v1.4  added function MAXRAD and calculations of J and Zp (ADS)
// 2013/09/01 - v1.4B corrected function MAXRAD and calculations of J and Zp (ADS)
// 2013/09/01 - v1.4B corrected function MAXRAD and calculations of J and Zp (ADS)
// 2019/12/22 - v1.5  Basic formatting and fix line endings (Windows vs. OSX vs. Unix) (GT)
// 2019/12/22 - v1.6  Fixed axis drawing (GT)