var  {declaration of global variables}
  x,y,picno:integer;
  Sn,Sx,Sy,Sxx,Syy,Sxy,threshu,threshl:integer;
  Cx,Cy,Mxx,Myy,Mxy,Theta,rot2,Parea:integer;
  M1,M2,Theta2,R1,R2,anis,bulk:integer;
  xmin,ymin,xmax,ymax:integer;
  tsumopix:integer;
  meen,moad,myn,maks:integer;
  s1,units:string;
  CA,TA,Ix,Iy,Imax,Imin,Rmyn,Rmaks,scaler:real;
  Answer,blah,name:string;
  Samplenum,BigX,BigY: integer;
  Zx,Zy,Zmax,Zmin,Xmaxrad,Xminrad,Ymaxrad,Yminrad: real;
  centrex,centrey,leny: integer;

procedure CalcSums;   {Calculate sum of pixel values and their squares}
begin
  Sn:=0; Sx:=0; Sy:=0; Sxx:=0; Syy:=0; Sxy:=0; tsumopix:=0;
  For y:=1 to ymax do
    begin
      For x:=1 to xmax do
        begin
          if GetPixel(x,y)=255 then
            begin
              Sn:=Sn+1;
              Sx:=Sx+(x*cos(rot2)+y*sin(rot2));
              Sy:=Sy+(y*cos(rot2)-x*sin(rot2));
              Sxx:=Sxx+sqr((x*cos(rot2)+y*sin(rot2)));
              Syy:=Syy+sqr((y*cos(rot2)-x*sin(rot2)));
              Sxy:=Sxy+((y*cos(rot2)-x*sin(rot2))*(x*cos(rot2)+y*sin(rot2)));
            end;
        end;
    end;
end;

Macro 'Help [F1]';
begin
  picno:=PidNumber;
  SelectWindow('Info');
  s1:='Select object of interest using the wand tool.';
  s1:=concat(s1,'  After selecting object, execute the Moment Calculation Macro');
  s1:=concat(s1,' from the Special Menu or by pressing [F2].  Results are written');
  s1:=concat(s1,' to a text file called "Measurement Results".\--> Click mouse');
  s1:=concat(s1,' button to continue...');
  ShowMessage(s1);
  SelectWindow('Info');
  Repeat until button;
  s1:='After executing the Moment Calculation Macro, you may execute';
  s1:=concat(s1,' the DrawAxis routine [F3] to draw the primary axis of the object.');
  ShowMessage(s1);
  SelectWindow('Info');
  wait(1);
  Repeat until button;
  if picno<>0 then SelectPic(picno);
end;


Macro 'Moment Calculation [F2]';

{This macro calculates the anisometry and bulkiness of binary particles
 based on the dynamically equivalent ellipse as defined by Medalia (1970).}

begin
  if (length(units)=0) or (scaler=0) then
    begin
      scaler:=1;
      units:='pixels';
    end;

  GetRoi(xmin,ymin,xmax,ymax);	{ymin is measured from top,xmin from left}
				{xmax & ymax= diameters along x & y axes}
  Copy;
  SetNewSize(xmax+2,ymax+2);
  MakeNewWindow('Calculating Moments');
  Paste;
  name:= GetString('Enter name of sample, max 10 chars: ');
  leny:=length(name);
  threshl:=GetNumber('Lower Threshold:  ',80,0);   {where actual thresholding}
  threshu:=GetNumber('Upper Threshold:  ',255,0);	{is done}
  if threshu=255 then changevalues(255,255,254);
  SetDensitySlice(threshl,threshu);  {set as desired}
  MakeBinary;
  Rot2:=0;
  CalcSums;
  If Sn=0 then begin PutMessage('Sn=0'); exit; end;
  Cx:=Sx/Sn+xmin-1;  {x-coord. of Centroid}
  Cy:=Sy/Sn+ymin-1;  {y-coord. of Centroid}
  Centrex:=Cx;
  Centrey:=Cy;

  		{following  code calculates y (dist from neutral axis)}

  if (((ymax + ymin) - Cy) > (Cy - ymin))	then
      begin
            Ymaxrad:= ymax + ymin - Cy;
            Yminrad:= Cy - ymin;
            BigY:= Ymaxrad;
      end
   else
	begin
		Ymaxrad:= Cy - ymin;
		Yminrad:= ymax + ymin - Cy;
		BigY:= Ymaxrad;
	end;

		{following code calculates x (dist from neutral axis)}

  if (((xmax + xmin) - Cx) > (Cx - xmin))	then
      begin
	      Xmaxrad:= xmax + xmin - Cx;
	      Xminrad:= Cx - xmin;
            BigX:= Xmaxrad;
	end
    else
      begin
	      Xmaxrad:= Cx - xmin;
	      Xminrad:= xmax + xmin - Cx;
            BigX:= Xmaxrad;
	end;

  Xmaxrad:= Xmaxrad/scaler;		{calibrating radii}
  Xminrad:= Xminrad/scaler;
  Ymaxrad:= Ymaxrad/scaler;
  Yminrad:= Yminrad/scaler;
  Parea:=Sn;
  Myy:=Sxx-(Sx*Sx/Sn);
  Mxx:=Syy-(Sy*Sy/Sn);
  Mxy:=Sxy-(Sx*Sy/Sn);
  If Mxy=0 then Theta:=0 else Theta:=arctan(((Mxx-Myy)+sqrt(sqr(Mxx-Myy)+
         (4*sqr(Mxy))))/(2*Mxy))*180/3.141592654;

  CA:= Parea/sqr(scaler);
  Cx := Cx/scaler;				{save x-coord of centroid}
  Cy := Cy/scaler;				{save y-coord of centroid}
  BigX:= BigX/scaler;				{calibrating x,y,Cx,Cy}
  BigY:= BigY/scaler;

  Ix:= Mxx/(sqr(sqr(scaler))); 		{save moment about x-axis}
  Iy:= Myy/(sqr(sqr(scaler)));
  Zx:= Ix/BigY;					{save section moduli}
  Zy:= Iy/BigX;

  SelectWindow('Calculating Moments');
  rot2:=theta*3.141592654/180;
  CalcSums;
  M1:=Sxx-(Sx*Sx/Sn);
  M2:=Syy-(Sy*Sy/Sn);
  R1:=sqrt(M1/Parea);
  R2:=sqrt(M2/Parea);

  Imax:= M1/(sqr(sqr(scaler)));
  Imin:= M2/(sqr(sqr(scaler)));
  Rmaks:= R1/scaler;
  Rmyn:= R2/scaler;

  fill;
  SetDensitySlice(threshl,threshu);
  MakeBinary;
  Rot2:=0;
  CalcSums;
  Theta:= -Theta;

If Sn=0 then begin PutMessage('Sn=0'); exit; end;

TA:= Sn/sqr(scaler);

SelectWindow('Calculating Moments');
Dispose;
RevertToSaved;


    Measure;

    rArea[1]:= TA;
    rMean[1]:= CA;
    rStdDev[1]:= Cx;		{saving important data into built-in arrays}
    rX[1]:= Cy;
    rY[1]:= Ix;			{names of arrays have no significance}
    rLength[1]:= Iy;
    rMajor[1]:= Imax;
    rMinor[1]:= Imin;
    rAngle[1]:= Theta;
    rMin[1]:= Zx;       	{Image doesn't allow user-defined arrays}
    rMax[1]:= Zy;
    rUser1[1]:= Xmaxrad;
    rUser2[1]:= Ymaxrad;

{displays the results in the upper left corner of the text window}

  SetFont('Monaco');
  SetFontSize(8);
  SetText('Plain; Align Left');
  CopyResults;

  NewTextWindow('Measurement Results',800,50);
  MoveWindow(20,395);
  Write('Name      ');
  Write('   TA     ');
  Write('   CA     ');
  Write('    Xbar     ');
  Write('  Ybar     ');
  Write('  Ix     ');
  Write('   Iy     ');
  Write(' Imax   ');
  Write('   Imin   ');
  Write('   Theta  ');
  Write('   Zx     ');
  Write('   Zy   ');
  Write('   Max x rad');
  Writeln('  Max y rad');
  Write(name);
  for x:=1 to (10 - leny) do
     begin
 	   Write(' ');
     end;
  Paste;
  ResetCounter;

{displays max and min radii about x- and y- axes}

  NewTextWindow('Max and Min Radii',300,40);
  MoveWindow(40,490);
  Write('Sample    ');
  Write(' max x rad ');
  Write(' min x rad ');
  Write(' max y rad ');
  Writeln(' min y rad');
  Write(name:10);
  Write(Xmaxrad:11:3);
  Write(Xminrad:11:3);
  Write(Ymaxrad:11:3);
  Writeln(Yminrad:11:3);

end;

Macro 'DrawAxis [F3]';

{optional macro to draw major/minor axis of ellipse}

begin
  MoveTo(Centrex,Centrey);
  SetForegroundColor(threshl-2);
  LineTo(Centrex-(cos((0-Theta)*3.141592654/180)*2*R1),
         Centrey+(sin((0-Theta)*3.141592654/180)*2*R1));
  MoveTo(Centrex,Centrey);
  LineTo(Centrex-(cos((Theta+90)*3.141592654/180)*2*R2),
         Centrey-(sin((Theta+90)*3.141592654/180)*2*R2));
  MoveTo(Centrex,Centrey);
  LineTo(Centrex+(cos((0-Theta)*3.141592654/180)*2*R1),
         Centrey-(sin((0-Theta)*3.141592654/180)*2*R1));
  MoveTo(Centrex,Centrey);
  LineTo(Centrex+(cos((Theta+90)*3.141592654/180)*2*R2),
         Centrey+(sin((Theta+90)*3.141592654/180)*2*R2));
  MoveTo(Centrex,Centrey);
  SetForegroundColor(255);
  SetDensitySlice(threshl-3,threshl-1);
end;

Macro 'SetUnits [F4]';

{optional macro to set the units to be used in the results for the Moment Calculation Macro.}

begin
  units:=GetString('Enter alternate units as desired:','pixels');
  scaler:=GetNumber(concat('Enter number of pixels per ',units,':'),1,1);
end;


{References:

Russ, John C.  Computer-Assisted Microscopy:  The Meaurement and Analysis
of Images.  Plenum Press:  New York, NY (1990).

Medalia, Avrom I.  Dynamic Shape Factors of Particles.  Powder Technology.
Vol. 4 (1970/71), p. 117-138.

}

{Written by:  Matthew Warfel - Cornell University - 4/4/97}
{Modified by: Stanley Serafin- Johns Hopkins University - 6/30/00}