{

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 S<Avg THEN FlipVertical;
		ShowMessage (Avg,'   ',S);
		
		Duplicate ('Operator');
		WOpPid := PidNumber;

		{Do some stuff with the Operator pic}

		SelectPic (WOpPid);
		Filter ('median'); Filter ('sobel');
		ChangeValues (1,15,0); ChangeValues (16,255,128);
		Filter ('max');
		GetPicSize (X2,Y2);
		GetColumn (1,1,Y2); FOR Y1:=1 TO Y2 DO LineBuffer[Y1]:=0;
		PutColumn (1,1,Y2);
		GetColumn (2,1,Y2); FOR Y1:=1 TO Y2 DO LineBuffer[Y1]:=0;
		PutColumn (2,1,Y2);

		GetColumn (X2,1,Y2); FOR Y1:=1 TO Y2 DO LineBuffer[Y1]:=0;
		PutColumn (X2,1,Y2);
		GetColumn (X2-1,1,Y2); FOR Y1:=1 TO Y2 DO LineBuffer[Y1]:=0;
		PutColumn (X2-1,1,Y2);
		GetColumn (X2-2,1,Y2); FOR Y1:=1 TO Y2 DO LineBuffer[Y1]:=0;
		PutColumn (X2-2,1,Y2);


		{Step 2, the calculation of the average density values
		 for each step of the wedge.
		 Go into the status machine to analyze the wedge steps}

		S := 0; SPrev := 0; Stepno := 0; Status:=-1; Dist:=0;
		Y1:=4; N:=0; Hicnt:=0; Locnt:=0; Dmax:=0;
		WHILE (Y1<=Y2) AND (Stepno<Maxsteps) DO
		BEGIN
			S:=0;
			SelectPic (WOpPid);
			GetRow (1,Y2+1-Y1,X2-1);
			FOR X1:=1 TO X2 DO BEGIN S:=S+LineBuffer[X1]; END;
			S := S/X2; {Average of one line of the Op img}
			Delta:=Abs(S-SPrev);
			
			{Status significance:
				-1: Not initialized - get first value
				0: Wait for Delta S threshold value - skip initial black area
				1: Threshold found, go to next flat area
				2: Do the maths for the flat area, back to 1
			}

			IF Status = -1 THEN BEGIN Status:=0; SPrev:=S; N:=0; END
			ELSE IF Status=0 THEN
			BEGIN
				{Skip as long as we don't reach threshold}
				Locnt:=Locnt+1; Hicnt:=0;
				FOR X1:=1 TO X2 DO LineBuffer[X1]:=80;
				PutRow (3,Y2+1-Y1,X2-6);
				IF (S>Thresh) 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 (S<Thresh) AND (Hicnt>Marg) 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 Max<LineBuffer[X] THEN Max:=LineBuffer[X];
				IF Min>LineBuffer[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<X2 and XR > XL. 
			 Otherwise we skip the line}

			IF (XL<X2) AND (XR>XL+4) AND (XL>7) AND (XR<X2-8) THEN
			BEGIN
				Dtissue := (Linebuffer[XL-3]+Linebuffer[XL-4]+Linebuffer[XL-5]
					   +Linebuffer[XR+3]+Linebuffer[XR+4]+Linebuffer[XR+5])/6;

				{From the current XL, advance until a maximum is found,
				 then do the same with the right side}

				XML := XL;
				WHILE (XML<XR) AND (XML<X2-4) AND 
				      ((LineBuffer[XML+1]<Linebuffer[XML]) OR 
					(LineBuffer[XML+2]<Linebuffer[XML]) OR
					(LineBuffer[XML+3]<Linebuffer[XML])) DO XML := XML+1;
				Max1 := Linebuffer[XML];
				XMR := XR;
				WHILE (XMR>XML) AND (XMR>3) AND 
				      ((LineBuffer[XMR-1]<Linebuffer[XMR]) OR 
					(LineBuffer[XMR-2]<Linebuffer[XMR]) OR
					(LineBuffer[XMR-3]<Linebuffer[XMR])) DO XMR := XMR-1;
				Max2 := Linebuffer[XMR];
				Differentiate := False;
				IF (XMR>XML) AND (XR>XML) AND (XL<X2-4) THEN Differentiate:=True;

				FOR X:=XL TO XR DO
				BEGIN
					IF LineBuffer[X]>D0 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]<R1) OR (X>XR);
					REPEAT
						Buf := Linebuffer[X];
						IF Buf<R1 THEN
						BEGIN
							IF Buf>D0 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 Buf<R2 THEN
						BEGIN
							IF Buf>D0 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;