// 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)