Dave Mitchell's DigitalMicrograph™ Scripting Website

// Function to find the centre of gravity of the passed in image. Cogx and cogy are passed in by reference and

// report the coordinates of the centre

// D. R. G. Mitchell, adminnospam@dmscripting.com (remove the no spam to make this address work)

// version:150807, v1.0, August 2015, www.dmscripting.com

// Function to find the centre of gravity of the passed in image (central roi). The coordinates of the centre are returned

// via the cogs and cogy variables, which are passed in by reference

void findcentreofgravity(image centralroi, number &cogx, number &cogy)

{

// For information on this calculation see John Russ - Image Processing Handbook, 2nd Ed. p489

number xpos, ypos, maxval, imgsum, xsize, ysize, i

image tempimg

maxval=max(centralroi, xpos, ypos)

imgsum=sum(centralroi)

getsize(centralroi, xsize, ysize)

// Traps for a blank image

if(imgsum==0) // the image is blank so set the CoGs to the centre of the image and return

{

cogx=(xsize-1)/2 //minus one since the centre of a 2 x 2 image is 0.5,0.5 - 0,0 is a position

cogy=(ysize-1)/2 //minus one since the centre of a 2 x 2 image is 0.5,0.5 - 0,0 is a position

return

}

// Collapse the image down onto the x axis

image xproj=realimage("",4,xsize,1)

xproj[icol,0]+=centralroi

// Rotate the passed in image through 90 degs so that rows become columns

// Then collapse that image down onto the x axis (was the y axis)

tempimg=realimage("", 4, ysize, xsize)

tempimg=centralroi[irow,icol]

image yproj=realimage("",4,ysize, 1)

yproj[icol,0]+=tempimg

yproj=yproj*(icol+1) // NB the +1 ensures that for the left column and top row, where

xproj=xproj*(icol+1) // icol=0 are included in the weighting. 1 must be subtracted from

// the final position to compensate for this shift

cogx=sum(xproj)

cogy=sum(yproj)

cogx=(cogx/imgsum)-1 // compensation for the above +1 to deal with icol=0

cogy=(cogy/imgsum)-1

}

// Main script - have an image front-most which contains some kind of blob or feature, which you

// wish to locate the centre of gravity of

image temp:=getfrontimage()

number cogx, cogy

findcentreofgravity(temp, cogx,cogy)

result("\nCog x="+cogx+" cogy="+cogy)