#ifndef lint
static char SccsId[] = "%W% %G%";
#endif
/* Module: prntcent.c (Print Intensity-Weighted Center of Region)
* Purpose: Compute a weighted center around a pixel
* Subroutine: centroid (xcent,ycent) returns: void
* Copyright: 1999 Smithsonian Astrophysical Observatory
* You may do anything you like with this file except remove
* this copyright. The Smithsonian Astrophysical Observatory
* makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without
* express or implied warranty.
* Note: Based on Mike VanHilst's print_table
* Modified: {0} Doug Mink initial version 25 October 1994
* {1} Doug Mink changed wcs args 7 July 1995
* {2} Doug Mink fixed centroid, wcs check 23 February 1996
* {3} Doug Mink declared set_center_param 24 January 1996
* {4} Paul Sydney drop unused xarray and yarray 9 July 1998
* {5} Doug Mink bump WCS string to 48 from32 6 May 1999
* {6} Doug Mink retyrn x and y 20 Aug 1999
* {n} <who> -- <does what> -- <when>
*/
#include <stdio.h> /* get stderr */
#include <math.h>
#include <X11/Xlib.h> /* X window stuff */
#include <X11/Xutil.h> /* X window manager stuff */
#include "hfiles/struct.h" /* declare structure types */
#include "hfiles/extern.h" /* extern main parameter structures */
/*
* Subroutine: print_center
* Purpose: Finds weighted file coordinate center around given pixel
* Currently uses 9x9 array
* Returns 1 if successful, else 0
* Note: Uses event coords in control struct (control.event.xkey)
*/
static int set_center_param();
static int get_key_buf_coord();
static void comp_center ();
void
print_center ( xc, yc )
double *xc, *yc; /* Returned centroided X and Y */
{
double dval;
double trow[25];
double tcol[25];
double xcent, ycent, xoff, yoff;
double xcbuf,ycbuf,xerr, yerr;
int x, y;
int nbox;
int rot;
int xbuf, ybuf;
int xbuf0, ybuf0; /* buffer coord of upper left in table */
int ival;
int i, j;
int clip;
int nval;
int get_pixel_val();
int xfile, yfile;
char string[64];
int lstr = 48;
void d_trans();
int iswcs();
*xc = 0.0;
*yc = 0.0;
/* determine the buffer coordinates of the event */
if (get_key_buf_coord (&control.event.xkey, &xbuf, &ybuf) == 0 ) {
(void)printf("SORRY: portion of image not in buffer.\n");
return;
}
/* note if there is something special about the values */
if (buffer.shortbuf_summing > 1 ) {
/* (void)printf("Image buffer has summed values:\n");
(void)printf("Each pixel summed from %d by %d square of pixels in file\n",
coord.fb.block, coord.fb.block);
(void)printf("Square starts at file coordinates shown in table\n"); */
return;
}
/* decide which pixels and get the relevant file coordinates */
nbox = 9;
rot = set_center_param (xbuf, ybuf, nbox, nbox,
&xbuf0, &ybuf0, &xfile, &yfile);
/* identify the pixel and note any coordinate oddities */
/* (void)printf("\nFile Column: %d, Row: %d\n", xfile, yfile); */
if (rot > 0 )
return;
/* Accumulate the marginal distributions */
nval = 0;
for (i=0; i<nbox; i++ )
trow[i] = 0.0;
for (j=0; j<nbox; j++ )
tcol[j] = 0.0;
y = xbuf0;
for (i=0; i<nbox; i++ ) {
x = xbuf0;
for (j=0; j<nbox; j++ ) {
if (get_pixel_val(x++, y, &ival, &dval, &clip) ) {
dval = (double) ival; /* convert if value is an integer */
}
trow[j] = trow[j] + dval;
tcol[i] = tcol[i] + dval;
nval = nval + 1;
}
y++;
}
/* Compute the center and estimate its error */
if (nval > 0) {
comp_center (tcol,nbox,&xoff,&xerr);
xcbuf = xbuf0 + xoff - 1.0;
comp_center (trow,nbox,&yoff,&yerr);
ycbuf = ybuf0 + yoff - 1.0;
d_trans (&coord.buftofile, xcbuf, ycbuf, &xcent, &ycent);
if (iswcs (wcs)) {
(void)pix2wcst (wcs,xcent,ycent,string,lstr);
(void)printf ("%s center from x,y= %.3f,%.3f (%d)\n",
string, xcent, ycent, nval);
}
else
(void)printf ("%.3f,%.3f is object center (%d)\n",xcent,ycent,nval);
*xc = xcent;
*yc = ycent;
}
else
(void)printf ("Cannot find center around %d %d\n",xfile,yfile);
return;
}
/*
* Subroutine: get_key_buf_coord
* Purpose: Determine the buffer coordinates of the pixel identified
* with a key event.
* Returns: 1 if the pixel is in the image data buffer, else 0.
*/
static int
get_key_buf_coord ( xkey, bufX, bufY )
XKeyEvent *xkey;
int *bufX, *bufY;
{
float x, y;
void i_transform(), d_transform();
/* translate event to buffer coordinates */
if( xkey->window == dispbox.ID ) {
i_transform(&coord.disptobuf, xkey->x, xkey->y, &x, &y);
} else if( xkey->window == panbox.ID ) {
i_transform(&coord.pantoimg, xkey->x, xkey->y, &x, &y);
d_transform(&coord.imgtobuf,(double)x, (double)y, &x, &y);
} else
return (0);
*bufX = (int)x;
*bufY = (int)y;
/* check if within buffer */
if ((*bufX < 0) || (*bufX >= coord.buf.width) ||
(*bufY < 0) || (*bufY >= coord.buf.height) )
return (0);
return (1);
}
/*
* Subroutine: set_center_param
* Purpose: Put pixval table parameters in the x and y arrays
* 0 to (dim-1) has file coords,
* _array[dim]: main file coord, _array[dim+1]: starting buf coord
*/
static int
set_center_param ( bufx, bufy, xdim, ydim,
xbuf, ybuf, xfile, yfile)
int bufx, bufy;
int xdim, ydim;
int *xbuf, *ybuf;
int *xfile, *yfile;
{
float x0, x1, x2, y0, y1, y2;
/* int xinc, yinc; */
int fx0, fx1, fx2, fy0, fy1, fy2;
int bufx0, bufy0, bufoff;
int rot;
void i_transform();
/* determine starting buffer position */
bufoff = (xdim - 1) / 2;
bufx0 = bufx - bufoff;
if( bufx0 < 0 )
bufx0 = 0;
else if( (bufoff = bufx0 + xdim - coord.buf.width) > 0 )
bufx0 -= bufoff;
bufoff = (ydim - 1) / 2;
bufy0 = bufy - bufoff;
if( bufy0 < 0 )
bufy0 = 0;
else if( (bufoff = bufy0 + ydim - coord.buf.height) > 0 )
bufy0 -= bufoff;
/* store the coordinates of the first buffer element */
*xbuf = bufx0;
*ybuf = bufy0;
/* set buf0 to offset from focus pixel */
bufx0 = bufx0 - bufx;
bufy0 = bufy0 - bufy;
/* determine file coordinates of focus and two other pixels */
i_transform (&coord.buftofile, bufx, bufy, &x0, &y0);
i_transform (&coord.buftofile, bufx+1, bufy, &x1, &y1);
i_transform (&coord.buftofile, bufx, bufy+1, &x2, &y2);
if( coord.file.ioff > 0.1 ) {
fx0 = (int)x0;
fy0 = (int)y0;
fx1 = (int)x1;
fy1 = (int)y1;
fx2 = (int)x2;
fy2 = (int)y2;
}
else {
fx0 = (int)(x0 + 0.5);
fy0 = (int)(y0 + 0.5);
fx1 = (int)(x1 + 0.5);
fy1 = (int)(y1 + 0.5);
fx2 = (int)(x2 + 0.5);
fy2 = (int)(y2 + 0.5);
}
/* store the focus coordinates */
*xfile = fx0;
*yfile = fy0;
/* determine the file increments along both axes and rotation */
if( (fx1 != fx0) && (fy1 == fy0) ) {
/* xinc = fx1 - fx0;
yinc = fy2 - fy0; */
rot = 0;
}
else if( (fx2 != fx0) && (fy2 = fy0) ) {
/* xinc = fx2 - fx0;
yinc = fy1 - fy0; */
rot = 1;
}
else {
rot = 2;
}
return( rot );
}
/* Compute the center and estimate its error given the marginal distribution
and the number of points. */
static void
comp_center (a,npts,coff,err)
double *a; /* array */
int npts; /* number of points */
double *coff; /* center value */
double *err; /* center error */
{
int i, npos;
double val, sumi, sumix, sumix2, di, derr;
/* Initialize */
npos = 0;
sumi = (double) 0;
sumix = (double) 0;
sumix2 = (double) 0;
/* Accumulate the sums */
for (i = 0; i < npts; i++) {
di = (double) (i + 1);
val = (double) a[i];
if (val > (double) 0) {
npos = npos + 1;
sumi = sumi + val;
sumix = sumix + (val * di);
sumix2 = sumix2 + (val * (di * di));
}
}
/* Compute the position and the error */
if (npos <= 0) {
*coff = (1.0 + (double) npts) / (double) 2;
*err = 0.0;
}
else {
*coff = sumix / sumi;
*err = (sumix2 / sumi - (*coff * *coff));
if (*err <= 0.0)
*err = 0.0;
else {
derr = *err / sumi;
*err = sqrt (derr);
if (*err > (double) npts)
*err = (double) 0;
}
}
}
syntax highlighted by Code2HTML, v. 0.9.1