COCCOLITH BIOMETRICS MACROS - PROGRAM NOTES 

A set of  macros for NIH-Image developed by Jeremy Young (Natural History Museum, London),  with much assistance from Michal Kucera (Geology Dept., Prague, Goteborg Univs.) and Hsiao-Wen Chung (Bioegineering, Pennsylvania Univ.). Purpose measuring Emiliania huxleyi coccoliths from cross-polarised light images. Queries to jy@nhm.ic.ac.uk. Some macros require Image version 1.54.
* After macro name indicates that operation is influenced by cursor position. 
** Indicates operation will be repeated on every tile of mosaic if cursor is in top left corner (red square). 
Revision 1.0b19, May 1994

This file consists of notes explaining the various macros, procedures and variables. It was formed by extracting comments from the macros file, then adding some extra notes. So it is organised in the same sequence as the macro file itself.   


var {Global Variables}
z:integer {main array reference} 
x,y,x1,x2,x3,x4,y1,y2,y3,y4:real;  {x-y pixel locations, used in various contexts}
cx,cy,xpos,ypos:integer;                         {x-y pixel locations, used for ellipse center and reporting proc results}
tcount,fcount,pcount:integer;                   {counters used by auto-measure macro}
lw,pitdist,oldpitdist,pitval:integer;       {variables for proclpit and procpit} 
N,n1,n2:integer;                                 {counters used in various contexts}
xMos,yMos,Nmos,xN,yN:integer;              {mosaic variables - tile x/y centre, ref no., x,y ref nos.} 
qn1,qn2:integer                                   {tile ref from procsequence};
{N.B. If it is not a mosaic image then xMos/yMos is a reference point set by the cursor location}
greyval:integer;               {used to record greyval of overwritten pixel}
edgethreshold:real;          {threshold ratio of min:(max-min) used by edge finding routine} 
radius:real;                     {radius used by proclpit, must be > max likely c.area length - microns} 
pt:integer                       {ca rim depth threshold - greyscale value}
tl:real;                            {length to test for procedge - microns}
maxrimcaratio:real;        {ratio of local ca radius to searched length for procedge}
cut:real;                          {threshold deviation from mean for rejecting edge points}
hscale,vscale,par:real;          {x and y spatial scales, pixel aspect ratio}
edgeL: real;                     {results of procedge - length (µm)}
pi: real;                           {constant pi, set by procinit}
xc,yc:real;                      {coordinate of ellipse center - microns}
xl,yl,xi,yi:real                 {image & local co-ords - pixels}
x8,y8,v:real                    {results from martin edge}
axisa,axisb:real;             {semi-axes of the ellipse - microns}
theta,k:real;                    {inclination of long axis - radians, clockwise from horiz.; k-various uses}
sin2th,sinth,costh:real;   {temp storage for trig fn values}
found,state,init,mosaic:boolean;


Procedure ShowParams
Displays various key customisable parameters in the Info window
Called by Set Parameters

Procedure procinit;
{function - sets a series of default parameter values} 
Called by many macros if not run already (boolean variable "init" records whether it has been).
Var
mag:integer;

{====================Mosaic Macros=============================}

Procedure procn2xy(Nmos);
{given the ref number Nmos finds the x,y loc of the centre of the mosaic tile}

Procedure procxy2n(x,y);
{given an x,y loc finds the ref number Nmos of the mosaic tile, & its centre xMos,yMos}

Procedure procSetSequence(x,y);
{function  - determine which specimen(s) to work on}
determines if it is a mosaic image from image size and then whether to work on only one or all specimens from cursor location.
   If cursor in top left corner then work on all specs}
   If cursor is elsewhere then work on spec in tile}
   If not a mosaic work on spec under cursor} 
var w,h:integer;
Returns xMos, yMos, qn1, qn2.

Macro 'TAKE MOSAIC PHOTO * [3]' 
{function - take "photo" of specimen and save to mosaic} 

Macro 'Make blank mosaic window [4]'
{function - set up blank mosaic window}

Macro 'Set number of next tile * [N]'
{function - set nMos, useful e.g. when compiling mosaics over a period of time} 

Macro 'Retrieve tile from duplicate * [D]' 
{function - retrieve a duplicate tile from another copy of the mosaic}
Useful when checking results. To use (1) close all image windows except the mosaic you are working on (2) Open a copy of the mosaic window. (3) Now whenever the macro is run it will cut and paste a tile from the background image into the active image. The tile pasted is, of course, determined by the mouse location. 

Macro 'Centre mosaic specimen * [C]'
{function - correct position  of off-centre specimens}
Use - Place cursor over centre of specimen, run macro. Specimen is cut & pasted to move the specimen -centre  to the centre of the tile.

Macro 'Show tile centres** [Q]'
Put green pixel(s) in the centre of tile(s) - with context sensitivity from procSetSequence

{======== ========== Biometric Macros ============= ============ =============}

Procedure procEllipseDraw(axisa,axisb,cx,cy,alpha,n2,colour:real,restore:boolean);
var - ang:real;
{draws an elliptical string of n dots, at centre cx,cy, with orientation alpha (degrees), and axial lengths axisa, axisb (pixels) length}
Stores grey vals of pixels in rUser[n2+] and draws in ellipse in colour or erases t accordingto value of boolean. 

procedure procMartinEdge(a,b,x0,y0:real); 
var e,k,q,x3,y3,xs,ys,k2,q2,k3,q3:real;
Given a point x0,y0 calculates its distance to an ellipse of axial lengths a, b. This works in local co-ordinates - i.e. centre of ellipse is at 0,0 and ellipse axes are parallel to the co-ordinate sytem. Hence vertically oriented point x0,y0 has to be calculated first (using procImageTo Local)
Returns x8,y8, position (in local co-ordinates) of nearest point on ellipse (a,b) to x0,y0. & v distance between x0,y0 & x8,y8.

procedure procImageToLocal(xq,yq,alpha,cx,cy);
Translates  & rotates point xq,yq from image co-ordinates into local co-ordinates
Returns xl,yl 

procedure procLocalToImage(xl,yl,alpha,cx,cy);
Rotates & translates  point xl,yl  from local co-ordinates into image co-ordinates.

Procedure procDrawAxes(z1);
displays frm previously caculated and stored (byprocCalcPixelLocs) pixel locs
 ca length, as red line
 ca width, as green line
 rimwidth, as pale blue line}

Procedure procTogglepixelvals;
{for use after running main macro, to show/hide values}
Only works on the most recently measured specimen

procedure procedge(xt1,yt1,xt2,yt2);
{function - find edge of coccolith shield}
{Calculates a threshold value based on the max & min values and looks for series of pixels above it. N.B. Works from outside inwards, i.e. from background toward specimen}
xt1,xt2 is a point outside the coccolith. xt2,yt2 is a point on the ca edge of the coccolith.
var
trip,thresh,xtn1,ytn1:integer


procedure procpit(x1,y1,n:integer);
{finds lowest val in (2n+1)x(2n+1) box, jumps to it, repeats until at local inverse peak, colours it red}
  var
pkx,pky,oldval,xi,yi,val: integer;

Procedure ProcLMinimum (x1,y1,x2,y2); 
{function - find xy loc of minimum value on line}
This is a simple alternative to proclpit. It is used by SEMI-MANUAL MEASURE [6]

Procedure proclpit(x1,y1,x2,y2,pt1); 
var
pitc,peaksep,nn,oldval,val:integer;
{PROCLPIT notes- peaksep is the number of pixels to be checked beyond the local pit, adjust up/down to eliminate pseudo peaks/avoid overlap problems - nn is counter which increases toward peaksep. n2 is reduced value of n used after a likely peak is found. pt1 is threshold value of pit to trip the lower value of n. SOLUTIONS - large specimems with variable values inside centre - raise n & pt1. Small specimens with weak edges - lower pt1. Touching specs - lower n2.};
this is relatively sophsticated pit finding procedure to allow coccolith to be measured when they are close together. A simpler approach would be to just find the minimum value on the line x1,y1 - x2,y2.

Returns xpos,ypos - the location of a point on the ca edge.{N.B. There may be a pixel/rounding error in the precise location of the xpos/ypos pixel};

procedure procEllipseFitHW;
{function - fit ellipse to a ring of points, given the location of the centre}
{N.B. Points come from proclpit & procpit, centre from analysis of pits} 
{pixel co-ordinates of the centre and boundary points are stored in arrays rX and rY, and are converted into micron co-ordinates with the centre as the origin in arrays rUser1 & rUser2}
var
  i,j:integer;                                 {array index general}
  a11,a12,a13:real;                       {for matrix A}
  a21,a22,a23:real;
  a31,a32,a33:real;
  det:real;                                     {determinant of matrix A}
  b1,b2,b3:real;                             {for vector B}
  i11,i12,i13:real;                        {for matrix A^(-1)}
  i21,i22,i23:real;
  i31,i32,i33:real;
  c1,c2,c3:real;                              {for coefficients}

Procedure procCalcPixelLocs(z1)
Begin
  {calculate ellipse co-ordinates and store in array(as pixel values)}
 
Procedure procCalcRimWidth(z);
 {determines a sensible rim width, from a set of local estimates (in microns). 
ave is average, only excluding very large values (>cal) ave2 also excludes values which deviate by > 'cut' from ave}
Var - na,na2,ave,ave2:real;
Output
 yellow pixels - used, & >1 pixel was within cut  within cut, yellow}
 green pixels - used locs, but 0 or 1 pixels were inside cut
 red pixels - values excluded. 
 rim width - rUser1[z]


Procedure procLineError;
Give error message if no line selection

Macro 'AUTO-MEASURE** [5]'
{function - automatically measure central area length and width of coccolith, rim width, and orientation. Imput a single pint inside the central area of the cocolith}
{rX[],rY[],rLength[],tcount -x/y loc,pixel greyvals, counter -  for the pits on the ca edge
 rUser1[],rUser2[],rMax[],fcount -  x/y loc, pixel greyvals, counter - for other points on the ca edge
 before doing ellipse fit the ca edge points are appended to the list of pit points and the counter pcount applies to entire series.
N.B. tcount, fcount & pcount all run from 101 up since they form indices to a work area of the measurement arrays.}
var
ang,Dang:real;                                           {angle and Ĉangle for linepit procedure} 
loc,loc2,nloc,pixelradius,diag,l1,l2,lax,lax1,lax2:integer;
Sequence of operation (repeated for each specimen)
1. Does radial series of transects, of length 'pixelradius' from refernce point inside c.area, using prolpit to find carea rim. Then procpit to migrate from this point to a peak on the carim. Store locs of both points.
2. Finds the two most widely separated peaks and centre of line between them; this is robust to spurious sub-peaks and problems due to saturated pixels, and in trials seemed at least as good as average from all  peaks even on clean specimens
3. Converts all points found on carea rim into a single series and do ellipse fit on these (assuming centre calc. above), using procEllipsefitHW
4. Makes a series of rim width estimates using procedge and estimate rim width from these using procCalcRimWidth
5. Displays results using procDrawAxes, and procTogglePixelVals.

Macro'SEMI-MANUAL MEASURE * [6]'
{N.B. Input is a line region of interest - a line going through the ellipse axis- and a point inside the coccolith ring indicated by the centre of the line selection - i.e. centre of line selection must be inside the c.area. Output results contain central area parameters}
var -  dx,dy,x0,y0,x3,x4,y3,y4,:real;
- This macro is uses the simple routine proclMinimum to find the edge of the carea ellipse, rather than the complex routine, proclpit used by the AUTO-MEASURE macro, this means it should work easily on most clean specimens.
- Combined with "Set Rim Width [9]" it makes a convenient semi-manual technique. 
Sequence of operation
1. Finds minimum values on each half of the selected line, using procLMinimum. Calculates from these, c area length, specimen orientation, and centre position.  
2. Makes new line through centre, perpendicular to original line.
3. Finds minimum values on each half of the new line, using procLMinimum. Calculates from these, c area width.
4. Displays results, using procDrawAxes


Macro 'Correct long axis * [7]';
{function - manually input/correct measurements associated with the long axis} 

Macro 'Correct long axis * [8]';
{function - manually input/correct measurements associated with the long axis} 

Macro 'Correct Rim Width * [9]'
{function - manually input rim width, needs central area measurements to have been made already}
Takes stored data from measurements file plus cursor loc, calculates rim width using procMartinEdge and draws rim as an ellipse "parallel" tothe c.area ellipse. - (using procellipse draw).  Data saved is  rim width in microns, in rUser1 array. N.B. Not yet corrected for Pixel aspect ratio. 

Macro 'Draw Axes** [A]'
uses procdrawaxes and sequencer to draw summary of measurements

Macro 'DrawCaEllipse * [E]'
Takes stored data from measurements file and plots central area ellipse (using procellipse draw) 

Macro 'Toggle pixel vals [T]'
{displays/hides the pixels found by proclpit and procpit on the last measured specimen} Uses proctogglepixelvals

Macro 'Surface plot specimen * [f8]'
{function - produce surface plot of a single specimen -i.e. current ROI, or if no ROI current tile, or if not a mosaic 96x90 box centred on cursor location.}
- The 96x90 box size is size of mosaic tile excluding any influence of edge - this gives good plotting of features with near to background values.  
- Because the background is dark on my images the image is inverted before plotting.
- To save memory the surface plot is discarded after the macro is run

Macro 'Set parameters * [P]'
{function - to allow operator to change key parameters/tune operation for different species or samples}
var 
i: integer;
Note that there is extensive help for this macro in the Info window.

Macro 'Initialise [I]'
Runs procinit - i.e. loads all the constants. Shouldn't be needed during routine work as macros which one might reasonably expect to be run immediately on starting call procinit. It is often handy when developing macros.

{==============Utility macros=====================================}
{These are not all strictly part of the Coccolith biometrics macros set but useful in the work}

macro 'Kill ROI [f2]'

macro 'Magn + ScaleBar [f3]'
{function - set spatial scale, based on microscope lenses used}
Note that the constants in this macro need to be adjusted according to the hardware used.
- hscale:=mag*0.0964;  Mag is the objective lens maginification (x intermediate/optovar lens if present). hscale is the horizontal scale on the image in pixels/micron (find this out by taking an image of a graticule). Constant 0.0964 is the necesssary conversion for my system. 
- The values from SetScale are not actually used by the macros but they set the correct calibrations for the built-in measuring functions of NIH-Image.  
- To draw scale bar without going through the set magnification dialog then enter 2 in answer to dialog 'Draw Scale Bar? (1=Yes, 0=No)'. 

Macro 'Shrink ROI [f5]'

Macro  'Select ROI * [f6]'
Makes a 100x100 Region of Intersts centered on the cursor.

Macro 'Grow ROI [f7]'

macro 'Reset Greymap[f9]'

macro 'Open Grey + 4 ramp LUT[f10]'; 
{this is an LUT of mine which is useful for checking details in dark end of greyscale 

macro 'Add green line to LUT [f11]'
Use to create an LUT with green reserved entry but noothers - so that saturated pixels are NOT red but tile numbers are green. Need to turn off reserved entries first.


EXTRA MACROS - in the file extra macros. These have had to be deleted to save space in the main macro file but may be useful. In particular the test/demo macros may be worth reloading when experimenting with trying out the macro on  different species. 

Macro 'CountSpecs[f4]'
{function - quick, on screen counting of specimens} 

MACRO 'test procpit [P]'
{test/demo macro for procpit. Input-cursor position}
Output - a red pixel on the found pit.

MACRO 'test proclpit [L]'
{test/demo macro for proclpit. Input - lineselection}
Output - red pixels on beginning & end of line used plus green pixel on found minimum.
var
pitc,nn:integer;

Macro 'test procedge [E]'
{demo/test macro for procedure procedge}
Output - grey pixels on beginning & end of line used plus yellow pixel on found edge.