{
                       BOX COUNTING MACRO

        The Box Count Macro is a macro that runs under NIHImage, v. 1.58
and later.  It performs the measurements that allows one to calculate the
so-called "capacity" fractal dimension.  The algorithm is called, in
fractal parlance, the "box counting" method.  In the simplest terms, the
routine counts the number of boxes of a given size that are needed to cover
a one pixel wide, binary (black on white) border.  The procedure is
repeated from boxes that are 2 to 64 pixels wide.  Two files are created.
One is called "Box Counts" and the other is called "Summary".  Either can
be saved as a text file for analysis with any application capable of curve
fitting and simple spread sheet manipulations, like KaleidaGraph.  The
Summary file is used to determine the fractal dimension (D).  It consists
of two columns labeled "size" and "#boxes".  One plots the log of size on
the x-axis and the log of #boxes on the y-axis and fits the data with a
power function.  The slope (S) of the line is the negative of the fractal
dimension, i.e. D= -S.

        A full description of the technique can be found in T. G. Smith,
Jr., G. D. Lange and W. B. Marks, Fractal Methods and Results in Cellular Morphology,
which appeared in J. Neurosci. Methods, 69:1123-126, 1996.

If you have any questions or problems, you may contact me, Tom Smith, at

        Bldg. 36 Room 2C02
        N.I.H
        Bethesda, MD 20892
        USA

or via e-mail at tgs@codon.nih.gov

}

var  {Global variables}
   size, MaxCount, ImagePid:integer;
   top, left, bottom, right:integer;
   MaxBoxSize :integer;
   tab,cr:string;


procedure GetHisto(size:integer);
var
   tTop, tLeft: integer;
   count, i: integer;
   done:boolean;
begin
   if ImagePid<>0 then ChoosePic(ImagePid);
   tLeft:=left;
   tTop:=top;
   ResetCounter;
   MaxCount:=sqr(size);
   SetUser1Label('Count');
   for i:=1 to MaxCount do rUser1[i]:=0;
   done:=false;
   repeat
      GetHistogram(tLeft, tTop, size,size);
      count:=histogram[255];
      if count>0 then
         rUser1[count]:=rUser1[count]+1;
      tLeft:=tLeft+size;
      if tLeft>=right then begin
         tLeft:=left;
         tTop:=tTop+size;
         done:=tTop>=bottom
      end;
   until done;
   SetCounter(MaxCount);
 end;

procedure WriteReal(r:real);
var
  e: integer;
begin
  e := 0;
  while r >= 10 do begin
    r := r / 10;
    e := e + 1;
  end;
  while r < 1 do begin
    r := r * 10;
    e := e - 1;
  end;
  write(r:8:6);
  if e >= 0 then
				write('e+')
  else begin
    e := -e;
				write('e-')
  end;
  if e <10 then
     write('0');
  write(e:1);
end;

procedure SaveCounts(size:integer);
var
  i, count, nBoxes, BoxSum:integer;
begin
   BoxSum := 0;
   for i := 1 to 10 do rUser2[i] := 0;
   SelectWindow('Box Counts');
   for count := 1 to MaxCount do begin
      nBoxes := rUser1[count];
      if nBoxes <> 0 then begin
         BoxSum := BoxSum + nBoxes;
         {
         rUser2[1] := rUser2[1] + count * nBoxes;
         rUser2[2] := rUser2[2] + count * count * nBoxes;
         rUser2[3] := rUser2[3] + count * count * count * nBoxes;
         rUser2[4] := rUser2[4] + count * count * count * count * nBoxes;
         rUser2[5] := rUser2[5] + count * count * count * count * count * nBoxes;
         rUser2[6] := rUser2[6] + (1/(count)) * nBoxes;
         rUser2[7] := rUser2[7] + (1/(count * count)) * nBoxes;
         rUser2[8] := rUser2[8] + (1/(count * count * count)) * nBoxes;
         rUser2[9] := rUser2[9] + (1/(count * count * count * count)) * nBoxes;
         rUser2[10]:=rUser2[10] + (1/(count * count * count * count * count)) * nBoxes;
         }
         write(size:2);
         write(tab);
         write(count:6);
         write(tab);
         write(nBoxes:6);
         write(cr);
      end;
   end;
   SelectWindow('Summary');
   write(size:2);
   write(tab);
   write(BoxSum:6);
   write(tab);
   {
   for i := 1 to 10 do begin
      WriteReal(rUser2[i]);
      if i <> 10 then write(tab);
   end;
   }
   write(cr);
end;

procedure DoBoxCounts(size:integer);
begin
   GetHisto(size);
   SaveCounts(size);
end;

procedure FindMargins;
var
   iWidth, iHeight:integer;
begin
   GetPicSize(iWidth, iHeight);
   ResetCounter;
   left:=-1;
   repeat  {Find left edge}
      left:=left+1;
      GetHistogram(left, 0, 1, iHeight);
   until histogram[255]<>0;
   ResetCounter;
   top:=-1;
   repeat  {Find upper edge}
      top:=top+1;
      GetHistogram(left, top, iWidth-left, 1);
   until histogram[255]<>0;
   ResetCounter;
   right:=iWidth+1;
   repeat  {Find right edge}
      right:=right-1;
      getHistogram(right-1, top, 1, iHeight-top);
   until histogram[255]<>0;
   ResetCounter;
   bottom:=iHeight+1;
   repeat  {Find bottom edge}
      bottom:=bottom-1;
      GetHistogram(left, bottom-1, right-left, 1);
   until histogram[255]<>0;
   if (iWidth-right<MaxBoxSize) or (iHeight-bottom<MaxBoxSize) then begin
      PutMessage('Right or bottom margin is less than ', MaxBoxSize:1, ' pixels.');
      exit;
   end;
   {MakeRoi(left, top, right-left, bottom-top);}
end;

procedure Setup;
begin
  SetFont('Monaco');
  SetFontSize(9);
  MaxBoxSize := 64;
  if Get('MaxMeasurements') < sqr(MaxBoxSize) then begin
     beep;
     PutMessage('Max Measurements (Options dialog box) must be set to at least ', sqr(MaxBoxSize):1, '.');
     exit;
  end;
  ImagePid:=PidNumber;
  FindMargins;
  NewTextWindow('Box Counts', 150, 400);
  MoveWindow(50, 50);
  tab:=chr(9);
  cr:=chr(13);
  write('size', tab, '#pixels', tab, '#boxes', cr);
  NewTextWindow('Summary', 150, 400);
  MoveWindow(210, 50);
  writeln('size', tab, '#boxes');
 {
  write('size', tab, '#boxes', tab, '  count^1   ', tab, '  count^2   ', tab, '  count^3   ', tab, '  count^4   ', tab, '  count^5   ', tab, ' count^-1  ', tab, '  count^-2  ', tab, '  count^-3  ', tab, '  count^-4  ', tab, '  count^-5  ', cr);
  }
end;


macro 'Do Box Counts /8';
begin
  SaveState;
  Setup;
  DoBoxCounts(2);
  DoBoxCounts(3);
  DoBoxCounts(4);
  DoBoxCounts(6);
  DoBoxCounts(8);
  DoBoxCounts(12);
  DoBoxCounts(16);
  DoBoxCounts(32);
  DoBoxCounts(64);
  RestoreState;
 end;