On this page, you can find the NIH Image macro for the algorithms described in the article "Computerized Methods for X-ray Based Small Bone Densitometry", published in Computerized Methods and Programs in Biomedicine, 73(1): 35-42 (2004)
Right-click on this link and select "save as" to download the macro
To download the free image processing program NIH Image (Macintosh patform only), click http://rsb.info.nih.gov/nih-image/index.html. For the PC/Windows version, visit the download page of Scion Corp.
{
Macro to calibrate bone density measurements with an aluminum wedge.
The first macro, Wedge, takes a ROI from the wedge and determines
two major parameters: the background density, D0, and the Al
X-ray attenuation, mu-Al.
The second macro, CalcDens, averages the image values in a given
ROI (bone area) and calibrates them according to the above values.
Rules of operation:
1. The aluminum wedge must be oriented vertically, the dividing
lines between the steps are horizontal.
2. If there are spots / speckles in the aluminum wedge area,
remove them (set them to the image value of the corresponding step)
3. When selecting the wedge ROI, stay away from the wedge borders
but provide enough area (make it at least 12 pixels wide)
4. Start the ROI way in the black area above the wedge. The black
"overhead" is required for the algorithm!
5. Make the ROI high enough to include at least 10 steps.
6. After the algorithm finished, check the "Wedge" window if the
values found are distribured correctly. Visually verify the
results!
}
{*****************************************************************}
VAR MuAl, D0: Real; {Final result variables: Al Xray attenuation, Background}
WedgePID, WOpPid, BonePid: Integer;
X,SX,SXX,SXXX,SXXXX,SY,SXY,SXXY: Real; {For the LSQ analysis}
Dmax: Real; {Darkest/lightest part of image}
Thick: Real; {Thickness of aluminum step}
Guessofs: Real; {Offset from last step to first D0 guess}
X1,Y1,X2,Y2,YSt,YEn: Integer;
S,SPrev,Thresh,Delta,Avg,Sqfac: Real;
Maxsteps,N,K,Stepno,Status,Dist: Integer;
Hicnt, Locnt, Marg: Integer;
mCount, mMean, mMode, mMin, mMax: Real;
Abort, CalibDone: Boolean;
{*****************************************************************}
Procedure LSQ (Dzero: Real);
VAR I: Integer;
Ni, Buf: Real;
BEGIN
SX:=0; SXX:=0; SXXX:=0; SXXXX:=0; SY:=0; SXY:=0; SXXY:=0;
Ni := Stepno; {a dummy to make it more transparent}
FOR I:=1 TO Stepno DO
BEGIN
Buf := LN (rUser2[I]-D0);
Sx := SX+Buf;
SXX := SXX + Buf*Buf;
SXXX := SXXX + Buf*Buf*Buf;
SXXXX := SXXXX + Buf*Buf*Buf*Buf;
SY := SY + rUser1[I];
SXY := SXY + Buf*rUser1[I];
SXXY := SXXY + Buf*Buf*rUser1[I];
END;
{Calculation of c from (a,b,c) = inv(A)*Y, see below}
Buf := Ni*SXX*SXXY + SX*SXY*SXX + SY*SX*SXXX - SXX*SXX*SY - SX*SX*SXXY - Ni*SXXX*SXY;
Buf := Buf / (Ni*SXX*SXXXX + 2*SX*SXX*SXXX - SXX*SXX*SXX - SX*SX*SXXXX - Ni*SXXX*SXXX);
Sqfac := Buf;
END;
{*****************************************************************}
Procedure CheckROI;
VAR Xa,Xb,Ya,Yb: Integer;
BEGIN
Abort:=False;
IF NOT Get ('ROIType') THEN
BEGIN
PutMessage ('Please specify the ROI first');
Abort:=True;
END
ELSE
BEGIN
GetRoi (Xa, Ya, Xb, Yb);
IF (Xb < 12) THEN
BEGIN
PutMessage ('Insufficient ROI width');
Abort:=True;
END
ELSE IF (Yb < 140) THEN
BEGIN
PutMessage ('Insufficient ROI height');
Abort:=True;
END;
END;
END;
Procedure CheckROI2;
BEGIN
Abort:=False;
IF NOT CalibDone THEN
BEGIN
PutMessage ('Perform the wedge calibration first!');
Abort:=True;
END
ELSE
BEGIN
IF NOT Get ('ROIType') THEN
BEGIN
PutMessage ('Please specify the ROI first');
Abort:=True;
END
END;
END;
{*****************************************************************}
Macro 'Wedge [W]';
BEGIN
{Specify some parameters, may be interactive later}
Thresh:=40; {Line sum threshold for recognition of step}
Marg:=3; {Minimum of abv/blw threshold to be a valid stel/flat}
Maxsteps:=6; {Number of wedge steps to be taken into account}
Thick := 0.5; {Wedge step thickness in mm}
Guessofs := 10; {Percent off the lowest step value to be used as first guess}
CalibDone:=False;
CheckROI;
IF NOT Abort THEN
BEGIN
{This is step 1: Create two child images from the
ROI, one of them for the edge finding algorithm.
Use a couple of filters to extract the edges and
clear the borders & do some debris cleanup}
COPY; {Temporary save ROI contents}
GetRoi (X1, Y1, X2, Y2);
SetNewSize (X2, Y2);
MakeNewWindow ('Wedge');
WedgePid := PidNumber;
Paste; KillROI;
SelectPic (WedgePid);
GetRow (1,1,X2-1);
FOR K:=1 TO X2 DO Avg:=Avg+LineBuffer[K];
GetRow (1,Y2-1,X2-1);
FOR K:=1 TO X2 DO S:=S+LineBuffer[K];
IF SThresh) AND (Locnt>Marg) THEN
BEGIN
Status:=1;
SelectPic (WedgePid);
MakeROI (3,Y2+1-Y1, X2-6,Y1-4);
Measure;
GetResults (mCount, mMean, mMode, mMin, mMax);
Dmax := mMean;
SetForegroundColor (0);
MoveTo (3,Y2+3-Y1); Write (mMean:3:0);
SetForegroundCOlor (255);
SelectPic (WOpPid);
END;
END
ELSE IF Status=1 THEN
BEGIN
{Skip while the threshold is exceeded, let YSt follow}
YSt := Y1;
Hicnt:=Hicnt+1; Locnt:=0;
FOR X1:=1 TO X2 DO LineBuffer[X1]:=255;
PutRow (3,Y2+1-Y1,X2-6);
IF (SMarg) THEN Status:=Status+1;
END
ELSE IF Status=2 THEN
BEGIN
Dist:=Dist+1; Locnt:=Locnt+1; Hicnt:=0;
IF (S>Thresh) AND (Locnt>Marg) THEN Status:=Status+1;
IF Status=3 THEN
BEGIN
YEn:=YSt+Dist-3;
SelectPic (WedgePid);
MakeROI (3,Y2+3-YEn, X2-6,Dist-4);
Measure; Stepno:=Stepno+1;
GetResults (mCount, mMean, mMode, mMin, mMax);
rUser1[Stepno] := Thick*Stepno;
rUser2[Stepno] := mMean;
IF mMean>100 THEN SetForegroundColor (0);
MoveTo (3,Y2+3-YEn); Write (mMean:3:0);
SetForegroundCOlor (255);
SelectPic (WOpPid);
Status:=1; Dist:=0;
END;
END;
Y1:=Y1+1; SPrev:=S;
END;
{Now we got the step wedge values (density over thickness)
in the rUser arrays [1..Stepno]. We guess an initial D0
10 % below the last found density value}
D0:= rUser2[Stepno] * (1-0.01*Guessofs);
{Step 3 is the iterative process to minimize the quadratic
factor of a parabolic least-squares fit by variation of D0.
So far, we've got x(i), the aluminum step thicknesses, and
y(i) = ln (D(i) - D0). Now we try to obtain a,b,c for the
parabolic equation yp = a + bx + cxx so that the sum of squares
(yp(x(i)) - y(I))**2 becomes minimal. If we solve this,
basically, we solve the following linear equation system:
| N SX SXX | | a | | SY |
| SX SXX SXXX | * | b | = | SXY |
| SXX SXXX SXXXX | | c | | SXXY |
Or, in shorter form, A * s = Y with s=(a,b,c),
A is a geometry-dependent 3x3 matrix, and Y depends both
on geometry and the guessed density. The goal is to
minimize c. The above equation may be rewritten with inv(A)
as the inverse of A as: s = inv(A) * Y. Cramer's rule
leads to the following solution of c:
| N SX SY | | N SX SXX |
c = det | SX SXX SXY | / det | SX SXX SXXX |
| SXX SXXX SXXY | | SXX SXXX SXXXX |
Only the first determinant changes with D0.}
LSQ (D0); S:=Sqfac;
REPEAT
D0 := D0-0.1; LSQ (D0);
Delta:=1;
IF ABS(Sqfac)<0.001 THEN Delta:=-1
ELSE IF ((s>0) AND (Sqfac<0)) OR ((S<0) AND (Sqfac>0)) THEN Delta:=-1
{ELSE Delta := ABS(S) - ABS(Sqfac)};
ShowMessage (D0:5:2,' ',Sqfac:7:3,' ',S:7:3,' ',Delta:7:3);
{FOR K:=1 TO 20000 DO ;}
S:=Sqfac;
UNTIL (Delta<0) OR (D0 <= 0);
{Done with D0? Then let's calculate MuAl out of the above}
MuAl := 0;
FOR K:=1 TO Stepno DO
BEGIN
MuAl := MuAl + (LN(Dmax-D0)-LN(rUser2[K]-D0))/rUser1[K];
END;
MuAl := MuAl / StepNo;
ShowMessage (Stepno:1,' steps; D0=',D0:3:0,' Dmax=',
Dmax:3:0,' \MuAl=',MuAl:5:2,' Sq=',Sqfac:5:2);
CalibDone := True;
END;
END;
{*****************************************************************}
Macro 'DiscardWnd [D]';
BEGIN
IF WedgePID < 0 THEN BEGIN SelectPIC (WedgePID); Dispose; WedgePid:=0; END;
IF WOpPID < 0 THEN BEGIN SelectPIC (WOpPID ); Dispose; WOpPid := 0; END;
IF BonePID < 0 THEN BEGIN SelectPIC (BonePID ); Dispose; BonePid:= 0; END;
END;
{*****************************************************************}
{Calculate the plain density of an area. No tissue correction,
just the aluminum wedge equivalents corrected for area}
Macro 'PlainDens [P]';
VAR X,Y: Integer;
Dtissue: Real;
Sum, BMA, Buf: Real;
Cnt: Integer;
BEGIN
Dtissue:=Dmax; {No tissue correction in this simple macro}
CheckROI2;
IF NOT Abort THEN
BEGIN
Cnt:=0; Sum:=0;
GetROI (X1,Y1,X2,Y2);
FOR Y:=Y1 TO Y1+Y2-1 DO
FOR X:=X1 TO X1+X2-1 DO
BEGIN
Buf := GetPixel (X,Y);
IF Buf>D0 THEN
BEGIN
Sum := Sum + LN (Dtissue-D0) - LN (Getpixel(X,Y)-D0);
Cnt := Cnt+1;
END;
END;
IF (Cnt=0) THEN
PutMessage ('Error: Area too small or too white')
ELSE
BEGIN
BMA := Sum / (Cnt*MuAl);
PutMessage ('Bone Area Density is ',BMA:7:2);
END;
END;
END;
{*****************************************************************}
{Analyze an area of bone with tissue correction. The bone is
supposed to be oriented vertically, and the rectangular ROI
must include large areas of tissue and background to the left
and right.}
Macro 'CalcDens [C]';
VAR I,X,Y,X1,Y1,X2,Y2,XL,XR,XML,XMR,PrevX: Integer;
Dtissue: Real;
Sum, BMC,BMS,BMT, Buf, Thresh: Real;
AreaC, AreaS, AreaT: Integer;
Cnt,Status: Integer;
R,R1,R2, Av, Max, Max1, Max2, Min, Sum: Real;
ThreshRatio, CortThres: Real;
Differentiate: Boolean;
BEGIN
{As before, we start with several parameters for easy access}
ThreshRatio:=0.4;
CortThres := 0.50;
BMC:=0; BMS:=0; BMT:=0; AreaC:=0; AreaS:=0; AreaT:=0;
CheckROI2;
IF NOT Abort THEN
BEGIN
COPY; {Temporary save ROI contents}
GetRoi (X1, Y1, X2, Y2);
SetNewSize (X2, Y2);
MakeNewWindow ('Bone');
BonePid := PidNumber;
Paste; KillROI;
SelectPic (BonePid);
Filter ('smooth more');
{Scan thru the image and read it line by line}
PrevX:=X2;
FOR Y:=1 TO Y2 DO
BEGIN
{Perform Basic statistics for each line}
GetRow (1,Y,X2-1);
Max:=LineBuffer[3]; Min:=Max; Avg:=0; Cnt:=0;
FOR X:=3 TO X2-4 DO
BEGIN
IF MaxLineBuffer[X] THEN Min:=LineBuffer[X];
Avg := Avg + LineBuffer[X]; Cnt := Cnt+1;
END;
IF (Cnt>0) THEN Avg:=Avg/Cnt ELSE Avg:=0;
{Start running from each of the sides to find
a threshold relative to the Max / Average.
Please keep in mind that bright areas carry small values}
XL:=1; Sum := 0; Cnt:=0;
REPEAT
XL:=XL+1; Buf := LineBuffer[XL];
Sum:=Sum+Buf; Cnt:=Cnt+1;
Thresh := ThreshRatio*Min + (1-ThreshRatio)*Sum/Cnt;
UNTIL (XL>=X2-1) OR ((Buf < Thresh) AND (Cnt>5));
XR:=X2-1; Sum:=0; Cnt:=0;
REPEAT
XR:=XR-1; Buf := LineBuffer[XR];
Sum:=Sum+Buf; Cnt:=Cnt+1;
Thresh := ThreshRatio*Min + (1-ThreshRatio)*Sum/Cnt;
UNTIL (XR<=XL) OR ((Buf < Thresh) AND (Cnt>5));
{Geometry check: Plausible only if XL XL.
Otherwise we skip the line}
IF (XLXL+4) AND (XL>7) AND (XRXML) AND (XMR>3) AND
((LineBuffer[XMR-1]XML) AND (XR>XML) AND (XLD0 THEN
BEGIN
BMT := BMT + LN (Dtissue-D0) - LN (Linebuffer[X]-D0);
AreaT:=AreaT+1;
END;
END;
IF Differentiate THEN
BEGIN
Min := LineBuffer[XL];
FOR X:=XML+1 TO XMR DO IF Min < LineBuffer[X] THEN Min:=LineBuffer[X];
Status:=0;
R1 := Min*(1-CortThres) + Max1*CortThres;
R2 := Min*(1-CortThres) + Max2*CortThres;
ShowMessage (XL:3,' ',XR:3,' ',min:3:0,' ',Avg:3:0,'\',
XML:3,' ',XMR:3,' ',Max1:3:0,' ',Max2:3:0,'\',
R1:3:0,' ',R2:3:0);
X:=XL; Cnt:=0;
REPEAT X:=X+1; UNTIL (LineBuffer[X]XR);
REPEAT
Buf := Linebuffer[X];
IF BufD0 THEN BEGIN BMC:=BMC + LN (Dtissue-D0) - LN (Buf-D0); AreaC:=AreaC+1; END;
LineBuffer[X]:=80;
END;
X:=X+1; Cnt:=Cnt+1;
UNTIL ((Buf>=R1) OR (X>XR) OR (X>PrevX) OR (Buf >= (Min+Max2)/2)) AND (Cnt>3);
PrevX:=X; Cnt:=0;
REPEAT
R := (R1*(X-XL) + R2*(XR-X))/(XR-XL); Buf:=LineBuffer[X];
IF (Buf>R) THEN
BEGIN
IF Buf>D0 THEN BEGIN BMS := BMS + LN (Dtissue-D0) - LN (Buf-D0); AreaS:=AreaS+1; END;
LineBuffer[X]:=180;
END;
X:=X+1; Cnt:=Cnt+1;
UNTIL ((Buf<=R2) AND (Cnt>2)) OR (X>XR);
Cnt:=0;
REPEAT
Buf:=LineBuffer[X];
IF BufD0 THEN BEGIN BMC:=BMC + LN (Dtissue-D0) - LN (Buf-D0); AreaC:=AreaC+1; END;
Cnt:=0; LineBuffer[X]:=100;
END;
X:=X+1; Cnt:=Cnt+1;
UNTIL ((Buf>R2) AND (Cnt>2)) OR (X>XR);
END;
{Finally, a visual confirmation of the segmentation process}
LineBuffer[XL-3]:=0; LineBuffer[XL-4]:=0; LineBuffer[XL-5]:=0;
LineBuffer[XR+3]:=0; LineBuffer[XR+4]:=0; LineBuffer[XR+5]:=0;
IF Differentiate THEN BEGIN LineBuffer[XML]:=255; LineBuffer[XMR]:=255; END;
PutRow (1,Y,X2-1);
END;
END;
IF (AreaT=0) THEN
PutMessage ('Error: Segmentation failed')
ELSE
BEGIN
BMT := BMT / (AreaT*MuAl);
IF AreaC=0 THEN BMC:=0 ELSE BMC:=BMC/(AreaC*MuAl);
IF AreaS=0 THEN BMS:=0 ELSE BMS:=BMS/(AreaS*MuAl);
ShowMessage ('Cortical density: ',BMC:7:3,'\',
'Spongiosa density: ',BMS:7:3,'\',
'Average (whole area): ',BMT:7:3);
PutMessage ('Bone Area Density is ',BMT:7:3,' ',
'Cortical part: ',BMC:7:3,' ',
'Spongy part: ',BMS:7:3);
END;
END;
END;