// This macro set, which requires ImageJ 1.43m or later, 
// can be used to measure the root to crown ratio in dental
// radiographs. It requires that the user draw and record
// three lines, called 'a', 'b' and 'c':
//
// 1. Draw line 'a' at the top of the tooth parallel to the cusps
//     and press 1 ("Record 'a'");
// 2. Drag line 'a' so it is positioned just underneath the tip of
//     the root and press 2  ("Record 'c'");
// 3. Draw line 'b' across the tooth at the junction of enamel
//     to cementum and press 3  ("Record 'b'").
//
// The macro displays the the root size to crown size ratio on 
// the radiograph and in the Results window. It calculates the
// ratio by finding the root size (distance between the midpoint
// of line 'b' and line 'c'), finding crown size (distance between
// the midpoint of line 'b' and line 'a'), and dividing.
//
// The lines are displayed in an non-destructive overlay. Press 4 
// to delete the lines drawn on the last tooth and the last Result
// table row. Press shift-f (Image>Overlay>Flatten) to create and
// RGB version of the radiograph with the lines rendered onto the
// image. Use the Image>Overlay>Hide Overlay command to
// remove the overlay.
// 
// There is an example digitised panoramic dental radiograph,
// courtesy of Minnie Bannister, available at
//   http://rsb.info.nih.gov/ij/images/DentalRadiograph.png
// 
// The function that calculates the distance between a point and
// a line segment is based on a Python example at stackoverflow.com.
// http://stackoverflow.com/questions/849211/shortest-distance-between-a-point-and-a-line-segment

   var ax1=-1, ay1, ax2, ay2;
   var bx1=-1, by1, bx2, by2;
   var cx1=-1, cy1, cx2, cy2;
   var startingOverlaySize = 0;
   var toothNumber = -1;

   macro "Record 'a' [1]" {
      requires("1.43m");
      if (selectionType!=5)
         exit("Stright line selection required");
      getLine(x1, y1, x2, y2, lineWidth);
      ax1=x1; ay1=y1; ax2=x2; ay2=y2;
      startingOverlaySize = Overlay.size;
      run("Add Selection...", "stroke=red width=1");
      toothNumber = -1;
  }

   macro "Record 'c' [2]" {
      getLine(x1, y1, x2, y2, lineWidth);
      cx1=x1; cy1=y1; cx2=x2; cy2=y2;
      slopeA = (ay2-ay1)/(ax2-ax1);
      slopeC = (cy2-cy1)/(cx2-cx1);
      if (abs(slopeA-slopeC)>0.01) {
          undo();
          exit("Line 'a' and line 'c' must have the same slope");
      }
      run("Add Selection...", "stroke=red width=1");
      run("Select None");
   }

   macro "Record 'b' and Calculate [3]" {
      Dialog.create("Tooth Number");
      Dialog.addNumber("Tooth number:", 0);
      Dialog.show();
      toothNumber = Dialog.getNumber();
      getLine(x1, y1, x2, y2, lineWidth);
      bx1=x1; by1=y1; bx2=x2; by2=y2;
      run("Add Selection...", "stroke=red width=1");
      run("Select None");
      calculate();
   }

   macro "Undo [4]" {
      undo();
   }

   function undo() {
      reset();
      currentOverlaySize = Overlay.size;
      for (i= currentOverlaySize-1; i> startingOverlaySize-1&& i>=0; i--)
          Overlay.removeSelection(i);
      if (toothNumber!=-1) {
          lastRow = nResults - 1;
          IJ.deleteRows(lastRow, lastRow);
      }
      toothNumber = -1;
   }

   function calculate() {
      if (ax1==-1 || bx1==-1 || cx1==-1) {
         reset();
         showMessage("Lines a, b and c not defined");
         exit;
      }
      px = bx1+(bx2-bx1)/2;
      py = by1+(by2-by1)/2;
      makeOval(px-3, py-3, 6, 6);
      run("Add Selection...", "fill=red");
      crown = distance(ax1, ay1, ax2, ay2, px, py);
      root = distance(cx1, cy1, cx2, cy2, px, py);
      ratio = root/crown;
      ax = ax1+(ax2-ax1)/2;
      ay = ay1+(ay2-ay1)/2;
      x = px+(ax-px)/2;
      y = py+(ay-py)/2;
      setFont("SansSerif", 18, " antialiased");
      makeText(d2s(ratio,2), x-15, y-15);
      run("Add Selection...", "stroke=red");
      run("Select None");
      row = nResults;
      setResult("Tooth no.", row, toothNumber);
      setResult("Root", row, root);
      setResult("Crown", row, crown);
      setResult("Ratio", row, ratio);
      updateResults;
      reset();
  }

   function reset() {
        ax1=-1; bx1=-1; cx1=-1;
  }

   function distance(x1, y1, x2, y2, px, py) {
      dx=x2-x1; dy=y2-y1;
      dd = sqrt(dx*dx+dy*dy);
      tx=dx/dd; ty=dy/dd;
      nx=-ty; ny=tx;
      acx = px-x1; acy=py-y1;
      return abs(acx*nx+acy*ny);
   }

// Python example from stackoverflow.com
//def pdis(a, b, c):
//    t = b[0]-a[0], b[1]-a[1]           # Vector ab
//    dd = sqrt(t[0]**2+t[1]**2)         # Length of ab
//    t = t[0]/dd, t[1]/dd               # unit vector of ab
//    n = -t[1], t[0]                    # normal unit vector to ab
//    ac = c[0]-a[0], c[1]-a[1]          # vector ac
//    return fabs(ac[0]*n[0]+ac[1]*n[1]) # Projection of ac to n (the minimum distance)