{Volume of rotation macros. Written by Jeremy  R. Young, April 1996. Regard this as emailware - if you use them I would be grateful for an email saying so, and any other feedback. jy@nhm.ac.uk
N.B. The original purpose is to enable accurate estimates of coccolith volume to be made.If used for that purpose please keep me informed of progess}

var
left, top, width,height,n,th,dx,cx,maxx,minx,minrad,maxrad,scale,LorR:integer;
a1,a2,b1,b2,elliptx,k,e,darea,dvol,area,vol,slicerad,mass,longaxis,shortaxis,density,pi,calcarea:real;
unit,str:string;

Macro 'Explanation [4]';
Begin;
PutMessage('The purpose of this macro set is to calculate the volume of an object with rotational symmetry (e.g. a cylinder, pear, torus) from its cross section. The cross section should be  red and the axis green');
End;


Macro 'Demo [5]'
Begin
SetNewSize(600,400);
MakeNewWindow('demo');
SetPalette('Greyscale',6);
SetScale(28.35,'cm');
SetForeGroundColor(1);
MakeOvalRoi(200,200,100,50);
Fill;
MakeOvalRoi(400,200,100,50);
Fill;
SetForeGroundColor(2);
MakeLineRoi(350,200,350,250);
Fill;
makeRoi (180, 180,400,100);
PutMessage('This has defined the cross section of torus(doughnut). Now run the main macro [1] to measure its volume. Results are given in the info window');
End;


Macro 'Extra features [6]';
Begin;

SetForegroundColor(1);
MakeOvalRoi(0,0,200,300);
DrawBoundary;
MakeOvalRoi(50,50,100,200);
DrawBoundary;
SetForegroundColor(2);
MakeLineRoi(100,0,100,150); Fill;
MoveTo(92,90); Write('a');

MakeLineRoi(100,150,200,150); Fill;
MoveTo(130,155); Write('b');

MoveTo(20,325); Writeln('Axial ratio = a/b');
Write('Note that separation of inner and outer edges is constant');

PutMessage('1. To measure elliptical shapes (with parallel margins) set axial ratio of OUTER margin in the dialog box. 2. By default the right half of  x-section is used. To use left half press option after  setting axial ratio.');
End;


Macro '(-';

Procedure procsetup;
BEGIN
SaveState;
SetForegroundColor(2);
SetOptions('area,mean,x-y center,Min/Max');
SetPrecision(3);
GetRoi(left, top, width,height);
pi:=3.14156;
density:=2.71 {this is value for calcite, in g/cm^3 OR pg/micron^3};
area:=0; vol:=0; 
GetScale(scale,unit);
If e=0 THEN e:=1.0;
ShowMessage('To use left half of profile press while clicking OK','\','\','Ellipticity is ratio of seen diam:diameter perpendicular to view','\','i.e. ellipticity of >1 implies you are viewing a section across the width of the ellipse');   
e:=GetNumber('If Left Profile press Opt \ Ellipticity =',e);
{ellipticity is defined from OUTER edge}
End;

Procedure procgetvalues;
BEGIN
     GetHistogram(dx,top,1,height);            {find thicknes of lith and display graphically} 
     th:=histogram[1];

     a1:=slicerad+0.5;  b1:=a1+k;
     a2:=slicerad-0.5;   b2:=a2+k;
     darea:=(pi *a1*b1 -  pi*a2*b2);         {calculate area & vol of slice} {++}
     dvol:=darea*th;
   area:=area+darea;                                            {calculate cumulative results}
   vol:=vol+dvol;
   PutPixel(dx,top-5-th,2);
   PutPixel(dx,top-5-dvol/500,6);
END;

Procedure procresults;
var
diam,kval:real;

BEGIN
SetCounter(rCount+1);
area:=area/sqr(scale);
vol:=vol/(sqr(scale)*scale);
mass:=vol*density;
calcarea:=calcarea/sqr(scale);
maxrad:=maxrad/scale;
minrad:=minrad/scale;
diam:=2*maxrad;
IF e>=1 then begin 
    longaxis:= diam*e;
    shortaxis:= diam;
 end else begin
    longaxis:=diam;
    shortaxis:=diam*e;
 end;
kval:=vol/(longaxis*longaxis*longaxis);  {kval is a shape constant}
ShowMessage('vol=',vol:10:2,unit,'3','\','mass=',mass:10:2,'pg','\','kval=',kval:10:3,'\','e=',e:10:2,'\','width=',shortaxis:10:2,unit,'\','length=',longaxis:10:2,unit,'\','\','screen diam=',diam:10:2,unit,'\','\','Green curve is thickness','\','Purple curve is volume');
{'area='area:10:2,unit,'2','\','calcarea= ',calcarea:10:2,'\','Outer radius = ',maxrad:10:2,unit,'\','Inner radius=',minrad:10:2,unit,'\',}
RestoreState;
rArea[rCount]:=vol;
rMean[rCount]:=mass;
{rstddev[rCount]:=mass;}
rX[rCount]:=kval;
rY[rCount]:=e;
rMin[rCount]:=shortaxis;
rMax[rCount]:=longaxis;
MakeRoi(left, top, width,height);
setPrecision(3);
END;


Macro 'Main Macro  [1]';
{function: given cross section of a coccolith calculate its volume}
{N.B. Data needed is a half cross section in red (i.e. 1) plus centre mark in green (i.e. 2)}
{Cross-section can be of left half, right half or total lith. To measure left half press option - before closing ellipticity  dialog box. For complete lith one or other half is measured according to which value is pressed}
{Cross-section should be along the short axis and ellipticity entered as long axis/short axis , i.e. values >1. If section parallel to long axis is used then ellpiticity should be entered as short/long axis (i.e. <1), but  errors will be introduced in central area volume calculation}

BEGIN
procsetup;
IF KeyDown('option') THEN LorR:=-1 ELSE LorR:=1;


{REPEAT}   {XXXXXXXXXXXXXXXXX   extra repeat for changing ellpticity loop}
IF LorR=-1 THEN dx:=left ELSE dx:=left+width;                         {++}


REPEAT                                                                          {find outer edge of profile}
    dx:=dx-(1*LorR);  {++}
     GetHistogram(dx,top,1,height);           
UNTIL histogram[1]>0;  maxx:=dx;

IF LorR=-1 THEN dx:=left+width ELSE dx:=left; {++}
REPEAT                                                                        {find marker for centre of profile}
    dx:=dx+(1*LorR);        {++}
     GetHistogram(dx,top,1,height);           
UNTIL histogram[2]>0;  cx:=dx;

REPEAT                                                                           {find inner  edge of profile}
      GetHistogram(dx,top,1,height);     
    dx:=dx+(1*LorR);        {++}    
UNTIL histogram[1]>0; minx:=dx;

maxrad:=(maxx - cx)*LorR;  {++}                                                 {calculate params relative to centre of profile}             
minrad:=(minx-cx)*LorR;        {++}
elliptx:=maxrad-minrad;
calcarea:= (pi * maxrad * maxrad*e) -  (pi * minrad  * (maxrad*e-elliptx));
k:=maxrad*e-elliptx-minrad;
MoveTo(cx,top-5);LineTo(maxx,top-5);

dx:=minx;
REPEAT
  dx:=dx+(1*LorR);
  slicerad:=(dx-cx)*LorR;
  procgetvalues;
UNTIL dx=maxx;

 procresults; 

{e:=e-0.1       }                          
{UNTIL e<.45}

END;

macro 'ScaleBar [2]'
{function - draw scale bar. N.B. It assumes the image has been scaled before, using the Set Scale.. command in NIH-Image }
   var
barlength:integer;
Begin
  GetPicSize (width,height);
  GetScale(scale,unit);
  str:=concat('Bar length (',unit,')');
  barlength:=GetNumber(str,5);
    MakeRoi (10,height-15, barlength*scale,5);
    SetForegroundColor(0);
    SetBackGroundColor(255);
   Fill;
   SetForeGroundColor(255);
   DrawBoundary;
    SetForegroundColor(0);
     MoveTo(barlength*scale+20,height-12);
    Write(barlength,unit);
  end;
end;

Macro 'Rotate ROI [3]'
{function - rotate specimen to vertical, angle needs to have been determined in advance}
var
angle:real;
Begin
angle:=GetNumber('Input Angle',5);
setscaling('bilinear same window');
ScaleandRotate(1,1,angle);
End;


macro 'erode [E]' begin erode; end; 
macro 'dilate [D]' begin dilate; end; 

MACRO 'Bootmenubar [M]'
{if the menubar is blackened this macro will usually restore it to a healthy grey}
var
test:integer;
Begin
test:=244;
{test:=GetNumber('Menu bar LUT value',244);}
{For test :=  244 to 244 do begin}
RedLUT[test]:=200;GreenLUT[test]:=200;BlueLUT[test]:=200;
UpdateLUT;
end;
end;