// Asteroid Shape Analysis Macro for ImageJ
// This macro was written by Rina Noguchi for Watanabe et al., Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu - a spinning-top-shaped rubble pile, Science (2019).
// This macro was modified after Liu et al., Optimising Shape Analysis to quantify volcanic ash morphology. GeoResJ (2015).

// INPUTS: this macro requires a binary image.

// OUTPUTS: Aspect ratio (minor/major axes diameter) and circularity.
// Mathematical definitions for the circularity given in SM7 of Watanabe et al. (2019).

// Default variables: particle size range = 750-Infinity pixels per particle; 'exclude' edge-touching particles; 
// 'include' holes in the particle area.
{
run("Set Measurements...", "area centroid perimeter bounding fit shape feret's redirect=None decimal=3");
run("Analyze Particles...", "size=1000-Infinity pixel Circularity=0.00-1.00 show=Outlines exclude clear include record stack");


    n = nResults;
    area1 = newArray(n);
    length1 = newArray(n);
    area2 = newArray(n);
    length2 = newArray(n);
    xstart = newArray(n);
    ystart = newArray(n);
    ff1 = newArray(n);
    AR1 = newArray(n);
    round1 = newArray(n);
    BBX1 = newArray(n);
    BBY1 = newArray(n);
    major1 = newArray(n);
    minor1 = newArray(n);
    cCirc = newArray(n);
    feret1 = newArray(n);
    minferet1 = newArray(n);
    solidity1 = newArray(n);
	a = newArray(n);
    for (i=0; i<n; i++) {
      area1[i] = getResult('Area', i);
      length1[i] = getResult('Perim.', i);
      xstart[i] = getResult('XStart', i);
      ystart[i] = getResult('YStart', i);
      ff1[i] = getResult('Circ.', i);
      AR1[i] = getResult('Aspect', i);
      round1[i] = getResult('Round', i);
      BBX1[i] = getResult('Width', i);
      BBY1[i] = getResult('Height', i);
      minor1[i] = getResult('Minor', i);
      major1[i] = getResult('Major', i);
      cCirc[i] = getResult('Circ', i);
      feret1[i] = getResult('Feret', i);
      minferet1[i] = getResult('MinFeret', i);
    }
	  run("RGB Color");
for(i=0; i<nResults; i++) {
x=getResult('X',i);
y=getResult('Y',i);
d=getResult('Major',i);
ah = getResult('Angle',i)*PI/180;
setColor("blue");
drawLine(x+(d/2)*cos(ah),y-(d/2)*sin(ah),x-(d/2)*cos(ah),y+(d/2)*sin(ah));
d=getResult('Minor',i);
ah=ah+PI/2;
setColor("red");
drawLine(x+(d/2)*cos(ah),y-(d/2)*sin(ah),x-(d/2)*cos(ah),y+(d/2)*sin(ah));
}
	
    run("Clear Results");
    for (i=0; i<n; i++) {
      doWand(xstart[i], ystart[i]);
      run("Convex Hull");
      run("Measure");
      area2[i] = getResult('Area', i);
      length2[i] = getResult('Perim.', i);
	  a[i] = getResult('Angle',i);
    }
run("Clear Results")
    run("Select None");
    for (i=0; i<n; i++) {
	setResult("Major axis", i, major1[i]);
	setResult("Minor axis", i, minor1[i]);
	setResult('Major axis angle',i, a[i]);
	setResult('Minor axis angle',i, a[i]-90);
	setResult("Equatorial aspect ratio", i, minor1[i]/major1[i]);
	setResult("Circularity", i, (2*(sqrt(PI*area1[i])))/length1[i]);
    }
     updateResults();
  }