{NIH Image Macro to calculate fractal dimension D of Roi}
{0-255 grey scale}
{Paolo B. DePetrillo, October 27, 1999}
{DePetrillo PB, Yang Q, Rackoff J, SanMiguel A, Karimullah K.}
{Surface fractal computation and its application to immunofluorescent} {histochemical studies of calpain and calpastatin in PC12 cells.}
{J Neurosci Meth 103(2):191-197, 2000.}

var
	i:integer
	left:integer
	top:integer
	w:integer
	h:integer
	numiter:real
	index:integer
	lefto:integer
	topo:integer
	hi:integer
	wi:integer
	wo:integer
	ho:integer
	size:integer
	sindex:integer
	pixCount:integer
	small:integer
	leftb:integer
	topb:integer
	wb:integer
	hb:integer
	mean:real
	meanxx:real
	mcount:real
	sindexr:real
	dsd:real
	n:real
	sumx:real
	sumy:real
	sumy2:real
	sumx2:real
	sumxy:real
	b:real
	b1:real
	b2:real
	b3:real
	b4:real
	b5:real
	b6:real
	b7:real
	b8:real
	b9:real
	b10:real
	sxx:real
	sxy:real
	syy:real
	coeff:real
	a:real
	sigmab:real
	sigmaa:real
	see:real
	r2:real
	
{rUser1 will hold slope values and rUser2 will hold r2 values for the regression}
	
	procedure Fit;   {this is a least mean squares fit, with y=ln(stdev*area) and x=ln(area)}
						{used for multiple slope calculation}
  begin
  	b:=0; sumx:=0; sumy:=0; sumx2:=0;sumxy:=0; sumy2:=0;
  	numiter := rCount;
  	for index:=2 to numiter do begin
      sumx:=sumx+ln(rArea[index]);
      sumy:=sumy+ln(rStdDev[index]*rArea[index]);
      sumx2:=sumx2+(ln(rArea[index])*ln(rArea[index]));
      sumxy:=sumxy+(ln(rArea[index])*ln(rStdDev[index]*rArea[index]));
      sumy2:=sumy2+(ln(rStdDev[index]*rArea[index])*ln(rStdDev[index]*rArea[index]));
   	end;
   	numiter:=numiter-1;
   	sxx:=sumx2-(sumx*sumx/numiter);
   	sxy:=sumxy-(sumx*sumy/numiter);
   	syy:=sumy2-(sumy*sumy/numiter);
   	b:=sxy/sxx;
   	a:=((sumx2*sumy-sumx*sumxy)/numiter)/sxx;
   	coeff:=sxy/sqrt(sxx*sxy);
   	r2:=(sxy*sxy)/(sxx*syy);
   	see:=sqrt((sumy2 - a*sumy - b*sumxy)/(numiter-2));
   	sigmab:=see/sqrt(sxx);
   	sigmaa:=sqrt(sumx2/numiter);
   end;
   
   procedure Fits;   {this is a least mean squares fit, with y=ln(stdev*area) and x=ln(area)}
						{used for single slope caclulation}
  begin
  	b:=0; sumx:=0; sumy:=0; sumx2:=0;sumxy:=0; sumy2:=0;
  	numiter := rCount;
  	for index:=1 to numiter do begin
      sumx:=sumx+ln(rArea[index]);
      sumy:=sumy+ln(rStdDev[index]*rArea[index]);
      sumx2:=sumx2+(ln(rArea[index])*ln(rArea[index]));
      sumxy:=sumxy+(ln(rArea[index])*ln(rStdDev[index]*rArea[index]));
      sumy2:=sumy2+(ln(rStdDev[index]*rArea[index])*ln(rStdDev[index]*rArea[index]));
   	end;
   	sxx:=sumx2-(sumx*sumx/numiter);
   	sxy:=sumxy-(sumx*sumy/numiter);
   	syy:=sumy2-(sumy*sumy/numiter);
   	b:=sxy/sxx;
   	a:=((sumx2*sumy-sumx*sumxy)/numiter)/sxx;
   	coeff:=sxy/sqrt(sxx*sxy);
   	r2:=(sxy*sxy)/(sxx*syy);
   	see:=sqrt((sumy2 - a*sumy - b*sumxy)/(numiter-1));
   	sigmab:=see/sqrt(sxx);
   	sigmaa:=sqrt(sumx2/numiter);
   end;
   
    macro 'MakeRoi'
   
	{get user info and make an Roi}
begin
KillRoi;
hi:=GetNumber('Enter height of Roi',25,0);
wi:=GetNumber('Enter width of Roi',25,0);
MakeRoi(left,top,wi,hi);
end;
	
	macro 'Fractal, multiple slopes'

begin
sindex:=1;
{this resets counter for shrinking rectangle}
	ResetCounter;
	GetRoi(left,top,w,h);
{but we want to remember the original Roi}
	lefto:=left;topo:=top;wo:=w;ho:=h;small:=0;
{this shrinks the Roi and measures}
		while (w-2>0) and (h-2>0) do begin
			Measure;
			small:=small+1;
			InsetRoi(1);
			GetRoi(left,top,w,h);
		end;
		GetRoi(left,top,w,h);
		leftb:=left;topb:=top;wb:=w;hb:=h;
		Fit;
	{we have the first estimate of D}
	rUser1[sindex]:=b; rUser2[sindex]:=r2;
	sindex:=sindex+1;
   	
   	
{let's do it again, this time we shrink from the bottom left to the top right}
	KillRoi;
   	MakeRoi(lefto,topo,wo,ho);
   	ResetCounter;
   	i:=1;
	GetRoi(left,top,w,h);
		while (w-1 >1) and (h-1>1) do begin
			GetRoi(left,top,w,h);
			Measure;
			MakeRoi(left+1,top,w-1,h-1);
			GetRoi(left,top,w,h);
		end;
		Fit;
{we have the second estimate of D}
   		rUser1[sindex]:=b; rUser2[sindex]:=r2;
		sindex:=sindex+1;
   	
   	
{let's do it again, this time we shrink the rectangle from the top left to the bottom right}
    KillRoi;
    MakeRoi(lefto,topo,wo,ho);
   	ResetCounter;
   	i:=1;
	GetRoi(left,top,w,h);
		while (w-1 >1) and (h-1>1) do begin
			GetRoi(left,top,w,h);
			Measure;
			MakeRoi(left+1,top+1,w-1,h-1);
			GetRoi(left,top,w,h);
		end;
		Fit;
{we have our third estimate of D}
   		rUser1[sindex]:=b; rUser2[sindex]:=r2;
		sindex:=sindex+1;
   	
   	
{let's do it again, this time we shrink the rectangle from the top right to the bottom left}
	KillRoi;
    MakeRoi(lefto,topo,wo,ho);
   	ResetCounter;
   	i:=1;
	GetRoi(left,top,w,h);
		while (w-1 >1) and (h-1>1) do begin
			GetRoi(left,top,w,h);
			Measure;
			MakeRoi(left,top+1,w-1,h-1);
			GetRoi(left,top,w,h);
		end;
		Fit;
{we have our fourth estimate of D}
   		rUser1[sindex]:=b; rUser2[sindex]:=r2;
		sindex:=sindex+1;
   		
   		
{let's do it again, this time we shrink the rectangle from the bottom right to the top left}
   	KillRoi;
    MakeRoi(lefto,topo,wo,ho);
   	ResetCounter;
   	i:=1;
	GetRoi(left,top,w,h);
		while (w-1 >1) and (h-1>1) do begin
			GetRoi(left,top,w,h);
			Measure;
			MakeRoi(left,top,w-1,h-1);
			GetRoi(left,top,w,h);
		end;
		Fit;
{we have our fifth estimate of D}
   		rUser1[sindex]:=b; rUser2[sindex]:=r2;
		sindex:=sindex+1;
   		
{let's do it again, this time we shrink the rectangle from top to bottom}
   	KillRoi;
    MakeRoi(lefto,topo,wo,ho);
   	ResetCounter;
   	i:=1;
	GetRoi(left,top,w,h);
		while (w-2 >1) and (h-2>1) do begin
			GetRoi(left,top,w,h);
			Measure;
			MakeRoi(left+1,top+1,w-1,h-1);
			GetRoi(left,top,w,h);
		end;
		Fit;
{we have our sixth estimate of D}
   		rUser1[sindex]:=b; rUser2[sindex]:=r2;
		sindex:=sindex+1;
   		
{let's do it again, this time we shrink the rectangle from bottom to top}
   	KillRoi;
    MakeRoi(lefto,topo,wo,ho);
   	ResetCounter;
   	i:=1;
	GetRoi(left,top,w,h);
		while (w-2 >1) and (h-2>1) do begin
			GetRoi(left,top,w,h);
			Measure;
			MakeRoi(left+1,top,w-1,h-1);
			GetRoi(left,top,w,h);
			i:=i+1;
		end;
		Fit;
{we have our seventh estimate of D}
   		rUser1[sindex]:=b; rUser2[sindex]:=r2;
		sindex:=sindex+1;
   		
{let's do it again, this time we shrink the rectangle left to right}
   	KillRoi;
    MakeRoi(lefto,topo,wo,ho);
   	ResetCounter;
   	i:=1;
	GetRoi(left,top,w,h);
		while (w-2 >1) and (h-2>1) do begin
			GetRoi(left,top,w,h);
			Measure;
			MakeRoi(left+1,top+1,w-1,h-1);
			GetRoi(left,top,w,h);
		end;
		Fit;
{we have our eighth estimate of D}
   		rUser1[sindex]:=b; rUser2[sindex]:=r2;
		sindex:=sindex+1;
   		
{let's do it again, this time we shrink the rectangle right to left}
   	KillRoi;
    MakeRoi(lefto,topo,wo,ho);
   	ResetCounter;
   	i:=1;
	GetRoi(left,top,w,h);
		while (w-2 >1) and (h-2>1) do begin
			GetRoi(left,top,w,h);
			Measure;
			MakeRoi(left,top-1,w-1,h-1);
			GetRoi(left,top,w,h);
		end;
		Fit;
{we have our ninth estimate of D}
   		rUser1[sindex]:=b; rUser2[sindex]:=r2;

   	
   	{and the average is:}
	mean:=0;
	meanxx:=0;
	mcount:=0;
	i:=1;
	for i:=1 to sindex do begin
		if rUser1[i]<.9 then begin
				mcount:=mcount+1;
		end;
		if (rUser1[i]>0.9) or (rUser1[i]=0.9) then begin
				mean:=mean+rUser1[i];
				meanxx:=meanxx+sqr(rUser1[i]);
		end;
	end;
	sindexr:=sindex;
	n:=sindexr-mcount;
	dsd:=sqrt((n*meanxx-(mean*mean))/(n*(n-1)));
	mean:= mean/(sindexr-mcount);
	NewTextWindow('Data Measurements', 500, 600);
    SetFont('Monaco');
    SetFontSize(12);
    for i:=1 to sindex do begin
    Writeln(rUser1[i],'  ',rUser2[i])
    end;
    Writeln(mean,' ','±',dsd/sqrt(n),'  mean fractal dimension ± sem');
    KillRoi;
    MakeRoi(lefto,topo,wo,ho);
    
{this figures out the median pixel value}
	KillRoi;
    MakeRoi(lefto,topo,wo,ho);
    ResetCounter;
	pixCount :=0;
	Measure;
	{rUser1 array is now used for pixel grey values from 0 to 255 and rUser2 array for histogram values}
		for i := 0 to 255 do begin
			rUser1[i+1] := i;
			pixCount := pixCount + Histogram[i];
			rUser2[i+1] := pixCount
		end;
		i:=1;
		while (rUser2[i]<wo*ho/2) and (i+1<256) do begin
			i:=i+1
		end;
		Writeln(rMean[1],'  mean density','  ',i,'  'rUser2[i],'  ',wi,hi);
		Writeln(rUser1[i-1],'   median pixel density','  ', rUser2[i-1],'  occurred at this pixelcount');
end;

	
	macro 'Fractal, single slope'

begin
sindex:=1;
{this resets counter for shrinking rectangle}
ResetCounter;
	GetRoi(left,top,w,h);
{but we want to remember the original Roi}
	lefto:=left;topo:=top;wo:=w;ho:=h;small:=0;
{this shrinks the Roi and measures}
		while (w-2>0) and (h-2>0) do begin
			Measure;
			small:=small+1;
			InsetRoi(1);
			GetRoi(left,top,w,h);
		end;

   	
   	
{let's do it again, this time we shrink from the bottom left to the top right}
	KillRoi;
   	MakeRoi(lefto,topo,wo,ho);
   	i:=1;
	GetRoi(left,top,w,h);
		while (w-1 >1) and (h-1>1) do begin
			GetRoi(left,top,w,h);
			Measure;
			MakeRoi(left+1,top,w-1,h-1);
			GetRoi(left,top,w,h);
		end;

   	
   	
{let's do it again, this time we shrink the rectangle from the top left to the bottom right}
    KillRoi;
    MakeRoi(lefto,topo,wo,ho);
   	i:=1;
	GetRoi(left,top,w,h);
		while (w-1 >1) and (h-1>1) do begin
			GetRoi(left,top,w,h);
			Measure;
			MakeRoi(left+1,top+1,w-1,h-1);
			GetRoi(left,top,w,h);
		end;

   	
   	
{let's do it again, this time we shrink the rectangle from the top right to the bottom left}
	KillRoi;
    MakeRoi(lefto,topo,wo,ho);
   	i:=1;
	GetRoi(left,top,w,h);
		while (w-1 >1) and (h-1>1) do begin
			GetRoi(left,top,w,h);
			Measure;
			MakeRoi(left,top+1,w-1,h-1);
			GetRoi(left,top,w,h);
		end;

   		
   		
{let's do it again, this time we shrink the rectangle from the bottom right to the top left}
   	KillRoi;
    MakeRoi(lefto,topo,wo,ho);
   	i:=1;
	GetRoi(left,top,w,h);
		while (w-1 >1) and (h-1>1) do begin
			GetRoi(left,top,w,h);
			Measure;
			MakeRoi(left,top,w-1,h-1);
			GetRoi(left,top,w,h);
		end;

   		
{let's do it again, this time we shrink the rectangle from top to bottom}
   	KillRoi;
    MakeRoi(lefto,topo,wo,ho);
   	i:=1;
	GetRoi(left,top,w,h);
		while (w-2 >1) and (h-2>1) do begin
			GetRoi(left,top,w,h);
			Measure;
			MakeRoi(left+1,top+1,w-1,h-1);
			GetRoi(left,top,w,h);
		end;

   		
{let's do it again, this time we shrink the rectangle from bottom to top}
   	KillRoi;
    MakeRoi(lefto,topo,wo,ho);
   	i:=1;
	GetRoi(left,top,w,h);
		while (w-2 >1) and (h-2>1) do begin
			GetRoi(left,top,w,h);
			Measure;
			MakeRoi(left+1,top,w-1,h-1);
			GetRoi(left,top,w,h);
			i:=i+1;
		end;

   		
{let's do it again, this time we shrink the rectangle left to right}
   	KillRoi;
    MakeRoi(lefto,topo,wo,ho);
   	i:=1;
	GetRoi(left,top,w,h);
		while (w-2 >1) and (h-2>1) do begin
			GetRoi(left,top,w,h);
			Measure;
			MakeRoi(left+1,top+1,w-1,h-1);
			GetRoi(left,top,w,h);
		end;

   		
{let's do it again, this time we shrink the rectangle right to left}
   	KillRoi;
    MakeRoi(lefto,topo,wo,ho);
   	i:=1;
	GetRoi(left,top,w,h);
		while (w-2 >1) and (h-2>1) do begin
			GetRoi(left,top,w,h);
			Measure;
			MakeRoi(left,top-1,w-1,h-1);
			GetRoi(left,top,w,h);
		end;
		{fit the whole bunch of slopes}
		Fits;
   	
   	
	NewTextWindow('Data Measurements', 500, 600);
    SetFont('Monaco');
    SetFontSize(12);
    Writeln(b,' ','±',sigmab,'  mean fractal dimension ± sem','  r2 is  ',r2);
    KillRoi;
    MakeRoi(lefto,topo,wo,ho);
    
{this figures out the median pixel value}
	KillRoi;
    MakeRoi(lefto,topo,wo,ho);
    ResetCounter;
	pixCount :=0;
	Measure;
	{rUser1 array is now used for pixel grey values from 0 to 255 and rUser2 array for histogram values}
		for i := 0 to 255 do begin
			rUser1[i+1] := i;
			pixCount := pixCount + Histogram[i];
			rUser2[i+1] := pixCount
		end;
		i:=1;
		while (rUser2[i]<wo*ho/2) and (i+1<256) do begin
			i:=i+1
		end;
		Writeln(rMean[1],'  mean density','  ',i,'  'rUser2[i],'  ',wi,hi);
		Writeln(rUser1[i-1],'   median pixel density','  ', rUser2[i-1],'  occurred at this pixelcount');
end;