Index: /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/include/fiximids.h =================================================================== --- /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/include/fiximids.h (revision 35591) +++ /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/include/fiximids.h (revision 35591) @@ -0,0 +1,76 @@ +# include +# include +# include + +typedef struct { + unsigned int start; + unsigned int stop; + Coords coords; +} Mosaic; + +// to determine the image ID, we need time (tzero + trate) and photcode +typedef struct { + unsigned int imageID; + unsigned short NX; // image width + unsigned short NY; // image height + e_time tzero; // readout time (row 0) + unsigned char trate; // scan rate (100 usec/pixel) + short photcode; +} ImageSubset; + +/* global variables set in parameter file */ +# define DVO_MAX_PATH 1024 +char ImageCat[DVO_MAX_PATH]; +char *CATDIR; +int HOST_ID; +char *HOSTDIR; +char *IMAGES; +char *SINGLE_CPT; +int VERBOSE; +int RESET; +int UPDATE; +int PARALLEL; +int PARALLEL_MANUAL; +int PARALLEL_SERIAL; + +SkyRegion UserPatch; + +/***** prototypes ****/ +int main PROTO((int argc, char **argv)); +void ConfigInit PROTO((int *argc, char **argv)); + +void usage_fiximids PROTO(()); +void initialize_fiximids PROTO((int argc, char **argv)); +int args_fiximids PROTO((int argc, char **argv)); + +void usage_fiximids_client PROTO(()); +void initialize_fiximids_client PROTO((int argc, char **argv)); +int args_fiximids_client PROTO((int argc, char **argv)); + +// int update_dvo_fiximids_client PROTO((ImageSubset *image, off_t Nimage, FlatCorrectionTable *flatcorr)); + +ImageSubset *ImageSubsetLoad PROTO((char *filename, off_t *nimage)); +int ImageSubsetSave PROTO((char *filename, ImageSubset *image, off_t Nimage)); +Image *ImagesFromSubset PROTO((ImageSubset *subset, off_t N)); +ImageSubset *ImagesToSubset PROTO((Image *image, off_t N)); + +void lock_image_db PROTO((FITS_DB *db, char *filename)); +void unlock_image_db PROTO((FITS_DB *db)); +void create_image_db PROTO((FITS_DB *db)); +void set_db PROTO((FITS_DB *in)); +int Shutdown PROTO((char *format, ...) OHANA_FORMAT(printf, 1, 2) ); +void TrapSignal PROTO((int sig)); +void SetProtect PROTO((int mode)); +int SetSignals PROTO((void)); + +Image *load_images_fiximids PROTO((FITS_DB *db, off_t *Nimage)); +int update_dvo_fiximids PROTO((ImageSubset *image, off_t Nimage)); +int update_dvo_fiximids_parallel PROTO((SkyTable *sky, ImageSubset *image, off_t Nimage)); + +void update_catalog_fiximids PROTO((Catalog *catalog, ImageSubset *image, off_t *index, off_t Nimage)); +int fiximids_local_astrometry PROTO((float *posAngle, float *pltScale, double x, double y, Coords *mosaic, Coords *coords)); + +off_t getMosaicByTimes PROTO((unsigned int start, unsigned int stop, unsigned int *startMos, unsigned int *stopMos, off_t *indexMos)); +void sort_mosaic_times PROTO((unsigned int *S, unsigned int *E, off_t *I, off_t N)); +void initMosaics PROTO((ImageSubset *image, off_t Nimage)); +Mosaic *getMosaicForImage PROTO((off_t im)); Index: /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/ImageSubsetFixImageIDs.c =================================================================== --- /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/ImageSubsetFixImageIDs.c (revision 35591) +++ /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/ImageSubsetFixImageIDs.c (revision 35591) @@ -0,0 +1,213 @@ +# include "fiximids.h" + +# define GET_COLUMN(OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +ImageSubset *ImageSubsetLoad(char *filename, off_t *nimage) { + + int Ncol; + off_t i, j, k; + off_t Nrow; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + *nimage = 0; + ImageSubset *image = NULL; + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename); + return NULL; + } + + /* load in PHU segment (ignore) */ + if (!gfits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset header\n"); + fclose (f); + return NULL; + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n"); + gfits_free_header (&header); + fclose (f); + return NULL; + } + + ftable.header = &theader; + + // load data for this header + if (!gfits_load_header (f, &theader)) { + fclose (f); + return NULL; + } + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) { + fclose (f); + return NULL; + } + fclose (f); + + char type[16]; + + GET_COLUMN (imageID, "IMAGE_ID", int); + GET_COLUMN (NX, "NX", int); + GET_COLUMN (NY, "NY", int); + GET_COLUMN (tzero, "TZERO", float); + GET_COLUMN (trate, "TRATE", float); + GET_COLUMN (photcode, "PHOTCODE" , short); + + ALLOCATE (image, ImageSubset, Nrow); + for (i = 0; i < Nrow; i++) { + image[i].imageID = imageID[i]; + image[i].NX = NX[i]; + image[i].NY = NY[i]; + image[i].tzero = tzero[i]; + image[i].trate = trate[i]; + image[i].photcode = photcode[i]; + } + fprintf (stderr, "loaded data for %lld images\n", (long long) Nrow); + + free (imageID); + free (NX); + free (NY); + free (tzero); + free (trate); + free (photcode); + + *nimage = Nrow; + return image; +} + +int ImageSubsetSave(char *filename, ImageSubset *image, off_t Nimage) { + + off_t i, j, k; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + gfits_create_table_header (&theader, "BINTABLE", "IMAGE_SUBSET"); + + // an unsigned int needs to have bzero of 0x8000 + gfits_define_bintable_column (&theader, "J", "IMAGE_ID", "image ID", "tmp", 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "NX", "tmp", "tmp", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "NY", "tmp", "tmp", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "TZERO", "tmp", "tmp", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "TRATE", "tmp", "tmp", 1.0, 0.0); + gfits_define_bintable_column (&theader, "I", "PHOTCODE", "tmp", "tmp", 1.0, 0.0); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + unsigned int *imageID; + int *NX, *NY; + float *tzero, *trate; + short *photcode; + + // create intermediate storage arrays + ALLOCATE (imageID, unsigned int, Nimage); + ALLOCATE (NX, int, Nimage); + ALLOCATE (NY, int, Nimage); + ALLOCATE (tzero, float, Nimage); + ALLOCATE (trate, float, Nimage); + ALLOCATE (photcode, int, Nimage); + + // assign the storage arrays + for (i = 0; i < Nimage; i++) { + imageID[i] = image[i].imageID ; + NX[i] = image[i].NX ; + NY[i] = image[i].NY ; + tzero[i] = image[i].tzero ; + trate[i] = image[i].trate ; + photcode[i] = image[i].photcode; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "IMAGE_ID", imageID, Nimage); + gfits_set_bintable_column (&theader, &ftable, "NX", NX, Nimage); + gfits_set_bintable_column (&theader, &ftable, "NY", NY, Nimage); + gfits_set_bintable_column (&theader, &ftable, "TZERO", tzero, Nimage); + gfits_set_bintable_column (&theader, &ftable, "TRATE", trate, Nimage); + gfits_set_bintable_column (&theader, &ftable, "PHOTCODE", photcode, Nimage); + + free (imageID ); + free (NX ); + free (NY ); + free (tzero ); + free (trate ); + free (photcode); + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file for output %s\n", filename); + return FALSE; + } + + gfits_fwrite_header (f, &header); + gfits_fwrite_matrix (f, &matrix); + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + fclose (f); + fflush (f); + + return TRUE; +} + +// XXX drop these two? + +# if (0) +Image *ImagesFromSubset (ImageSubset *subset, off_t N) { + + off_t i; + + // we have been given an ImageSubset array, containing a reduced set of image fields + // create full a Image array and save the needed values + Image *image = NULL; + ALLOCATE (image, Image, N); + for (i = 0; i < N; i++) { + image[i].imageID = subset[i].imageID; + image[i].NX = subset[i].NX ; + image[i].NY = subset[i].NY ; + image[i].tzero = subset[i].tzero ; + image[i].trate = subset[i].trate ; + image[i].coords = subset[i].coords ; + } + return image; +} + +ImageSubset *ImagesToSubset (Image *image, off_t N) { + + off_t i; + + // we have been given an ImageSubset array, containing a reduced set of image fields + // create full a Image array and save the needed values + ImageSubset *subset = NULL; + ALLOCATE (subset, ImageSubset, N); + + for (i = 0; i < N; i++) { + subset[i].imageID = image[i].imageID; + subset[i].NX = image[i].NX ; + subset[i].NY = image[i].NY ; + subset[i].tzero = image[i].tzero ; + subset[i].trate = image[i].trate ; + subset[i].coords = image[i].coords ; + } + return subset; +} + +# endif (0) Index: /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/fiximids.c =================================================================== --- /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/fiximids.c (revision 35591) +++ /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/fiximids.c (revision 35591) @@ -0,0 +1,36 @@ +# include "fiximids.h" + +int main (int argc, char **argv) { + + off_t Nimage; + int status; + FITS_DB db; + + /* get configuration info, args, lockfile */ + initialize_fiximids (argc, argv); + + set_db (&db); + status = dvo_image_lock (&db, ImageCat, 60.0, (UPDATE ? LCK_XCLD : LCK_SOFT)); + if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename); + if (db.dbstate == LCK_EMPTY) Shutdown ("ERROR: No images in catalog %s (1)", db.filename); + + // load images (images are not modified, but we need the astrometry transformations + Image *image = load_images_fiximids (&db, &Nimage); + dvo_image_unlock (&db); + + // the imageSubset here is a reduced set of fields, not a reduced set of images + ImageSubset *subset = ImagesToSubset (image, Nimage); + + status = update_dvo_fiximids (subset, Nimage); + + if (!status) exit (1); + exit (0); +} + + +/* fiximids : set the posangle & pltscale for measurements in the db + ** load images (save subset with imageID, time, photcode, seq) + ** load catalogs + ** update detection + ** count detections / image ID? + */ Index: /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/fiximids_client.c =================================================================== --- /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/fiximids_client.c (revision 35591) +++ /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/fiximids_client.c (revision 35591) @@ -0,0 +1,24 @@ +# include "setposangle.h" + +// setposangle_client is run on a remote host and is responsible for updating the catalogs +// owned by that host. + +// setposangle_client -hostname (hostname) -catdir (catdir) -hostdir (hostdir) -images (images) +// load the SkyTable and HostTable from (catdir) [contains the host responsibilities] +// load the image data from (images) [contains the zp infomation to apply] +// loop over tables from SkyTable with my host ID +// the calling program tells the client which host they are (or host ID?); we do not actually have to run this on the real host + +int main (int argc, char **argv) { + + // get configuration info, args, lockfile (set CATDIR, HOST_ID, HOSTDIR, IMAGES) + initialize_setposangle_client (argc, argv); + + // load the image subset table from the specified location + off_t Nimage; + ImageSubset *subset = ImageSubsetLoad (IMAGES, &Nimage); + + update_dvo_setposangle (subset, Nimage); + + exit (0); +} Index: /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/initialize_fiximids.c =================================================================== --- /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/initialize_fiximids.c (revision 35591) +++ /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/initialize_fiximids.c (revision 35591) @@ -0,0 +1,267 @@ +# include "fiximids.h" + +void usage_fiximids () { + fprintf (stderr, "USAGE: fiximids [options]\n"); + fprintf (stderr, " options:\n"); + fprintf (stderr, " -v : verbose mode\n"); + fprintf (stderr, " -update : actually write results to detections tables\n"); + fprintf (stderr, " -parallel : run in parallel mode\n"); + fprintf (stderr, " -h : this help list\n"); + fprintf (stderr, " -help : this help list\n"); + fprintf (stderr, " --h : this help list\n"); + fprintf (stderr, " --help : this help list\n"); + fprintf (stderr, " Note that the dvo db can be specified by -D CATDIR (directory)\n"); + exit (2); +} + +void initialize_fiximids (int argc, char **argv) { + + int N; + struct stat statbuffer; + char CatdirPhotcodeFile[256]; + + if ((N = get_argument (argc, argv, "-h"))) usage_fiximids(); + if ((N = get_argument (argc, argv, "-help"))) usage_fiximids(); + if ((N = get_argument (argc, argv, "--h"))) usage_fiximids(); + if ((N = get_argument (argc, argv, "--help"))) usage_fiximids(); + + /* are these set correctly? */ + ConfigInit (&argc, argv); + args_fiximids (argc, argv); + + if (stat (CATDIR, &statbuffer)) { + fprintf (stderr, "error accessing dvo database directory '%s'\n", CATDIR); + exit (1); + } + + // load the photcode table + sprintf (CatdirPhotcodeFile, "%s/Photcodes.dat", CATDIR); + if (!LoadPhotcodes (CatdirPhotcodeFile, NULL, FALSE)) { + fprintf (stderr, "error loading photcode table %s\n", CatdirPhotcodeFile); + exit (1); + } +} + +int args_fiximids (int argc, char **argv) { + + int N; + + VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + // XXX I want to avoid running this on all data all the time + // some options: gpc1 only (hardwired), require photcodes (fairly tedious), + // exclude some photcodes? + // XXX for now, hardwire gpc1 photcodes + + SINGLE_CPT = NULL; + if ((N = get_argument (argc, argv, "-cpt"))) { + remove_argument (N, &argc, argv); + SINGLE_CPT = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + // region of interest + UserPatch.Rmin = 0; + UserPatch.Rmax = 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + if ((N = get_argument (argc, argv, "-region"))) { + remove_argument (N, &argc, argv); + UserPatch.Rmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Rmax = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmax = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + PARALLEL = FALSE; + if ((N = get_argument (argc, argv, "-parallel"))) { + PARALLEL = TRUE; + remove_argument (N, &argc, argv); + } + // this is a test mode : rather than launching the remote jobs and waiting for completion, + // fiximids will simply list the remote command and wait for the user to signal completion + PARALLEL_MANUAL = FALSE; + if ((N = get_argument (argc, argv, "-parallel-manual"))) { + PARALLEL = TRUE; // -parallel-manual implies -parallel + PARALLEL_MANUAL = TRUE; + remove_argument (N, &argc, argv); + } + // this is a test mode : rather than launching the fiximids_client jobs remotely, they are + // run in serial via 'system' + PARALLEL_SERIAL = FALSE; + if ((N = get_argument (argc, argv, "-parallel-serial"))) { + if (PARALLEL_MANUAL) { + fprintf (stderr, "ERROR: cannot mix -parallel-manual and -parallel-serial\n"); + exit (1); + } + PARALLEL = TRUE; // -parallel-serial implies -parallel + PARALLEL_SERIAL = TRUE; + remove_argument (N, &argc, argv); + } + + UPDATE = FALSE; + if ((N = get_argument (argc, argv, "-update"))) { + remove_argument (N, &argc, argv); + UPDATE = TRUE; + } + + if (argc != 1) usage_fiximids(); + + return (TRUE); +} + +void usage_fiximids_client () { + fprintf (stderr, "USAGE: fiximids_client -hostID (hostID) -catdir (catdir) -hostdir (hostdir) -images (images) [options]\n"); + fprintf (stderr, " options:\n"); + fprintf (stderr, " -v : verbose mode\n"); + fprintf (stderr, " -update : actually write results to detections tables\n"); + exit (2); +} + +void initialize_fiximids_client (int argc, char **argv) { + + int N; + struct stat statbuffer; + char CatdirPhotcodeFile[256]; + + if ((N = get_argument (argc, argv, "-h"))) usage_fiximids_client(); + if ((N = get_argument (argc, argv, "-help"))) usage_fiximids_client(); + if ((N = get_argument (argc, argv, "--h"))) usage_fiximids_client(); + if ((N = get_argument (argc, argv, "--help"))) usage_fiximids_client(); + + SetZeroPoint (25.0); // XXX is this needed? + args_fiximids_client (argc, argv); + + if (stat (CATDIR, &statbuffer)) { + fprintf (stderr, "error accessing dvo database directory '%s'\n", CATDIR); + exit (1); + } + + // load the photcode table : XXX needed? + sprintf (CatdirPhotcodeFile, "%s/Photcodes.dat", CATDIR); + if (!LoadPhotcodes (CatdirPhotcodeFile, NULL, FALSE)) { + fprintf (stderr, "error loading photcode table %s\n", CatdirPhotcodeFile); + exit (1); + } +} + +int args_fiximids_client (int argc, char **argv) { + + int N; + + VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + UPDATE = FALSE; + if ((N = get_argument (argc, argv, "-update"))) { + remove_argument (N, &argc, argv); + UPDATE = TRUE; + } + + // region of interest + UserPatch.Rmin = 0; + UserPatch.Rmax = 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + if ((N = get_argument (argc, argv, "-region"))) { + remove_argument (N, &argc, argv); + UserPatch.Rmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Rmax = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmax = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + HOST_ID = 0; + if ((N = get_argument (argc, argv, "-hostID"))) { + remove_argument (N, &argc, argv); + HOST_ID = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + if (!HOST_ID) usage_fiximids_client(); + + HOSTDIR = NULL; + if ((N = get_argument (argc, argv, "-hostdir"))) { + remove_argument (N, &argc, argv); + HOSTDIR = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if (!HOSTDIR) usage_fiximids_client(); + + if ((N = get_argument (argc, argv, "-catdir"))) { + remove_argument (N, &argc, argv); + CATDIR = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if (!CATDIR) usage_fiximids_client(); + + IMAGES = NULL; + if ((N = get_argument (argc, argv, "-images"))) { + remove_argument (N, &argc, argv); + IMAGES = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if (!IMAGES) usage_fiximids_client(); + + if (argc != 1) usage_fiximids_client(); + return (TRUE); +} + +void GetConfig (char *config, char *field, char *format, int N, void *ptr); + +void ConfigInit (int *argc, char **argv) { + + double ZERO_POINT; + char *config, *file; + + /*** load configuration info ***/ + file = SelectConfigFile (argc, argv, "ptolemy"); + config = LoadConfigFile (file); + if (config == (char *) NULL) { + fprintf (stderr, "ERROR: can't find configuration file %s\n", file); + if (file != (char *) NULL) free (file); + exit (0); + } + if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file); + + // force CATDIR to be absolute (so parallel mode will work) + char *tmpcatdir = NULL; + ALLOCATE (tmpcatdir, char, DVO_MAX_PATH); + GetConfig (config, "CATDIR", "%s", 0, tmpcatdir); + CATDIR = abspath (tmpcatdir, DVO_MAX_PATH); + free (tmpcatdir); + + sprintf (ImageCat, "%s/Images.dat", CATDIR); + + ScanConfig (config, "ZERO_PT", "%lf", 0, &ZERO_POINT); + SetZeroPoint (ZERO_POINT); + + free (config); + free (file); +} + +void GetConfig (char *config, char *field, char *format, int N, void *ptr) { + + char *status; + + status = ScanConfig (config, field, format, N, ptr); + if (status == NULL) { + fprintf (stderr, "error in config, cannot find %s\n", field); + exit (1); + } + return; +} Index: /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/load_images_fiximids.c =================================================================== --- /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/load_images_fiximids.c (revision 35591) +++ /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/load_images_fiximids.c (revision 35591) @@ -0,0 +1,174 @@ +# include "fiximids.h" + +// array of mosaic definition structures +static off_t Nmosaic; +static Mosaic *mosaic; + +// list of mosaic associated with each image +static off_t Nmosaic_for_images; // number of images (for off_ternal checks) +static off_t *mosaic_for_images; // array of: image -> mosaic + +Image *load_images_fiximids (FITS_DB *db, off_t *Nimage) { + + Image *image; + + if (VERBOSE) fprintf (stderr, "finding images\n"); + + /* read entire db table */ + if (!dvo_image_load (db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db[0].filename); + + /* use a vtable to keep the images to be calibrated */ + image = gfits_table_get_Image (&db[0].ftable, Nimage, &db[0].swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } + + fprintf (stderr, "loaded "OFF_T_FMT" images\n", *Nimage); + + return (image); +} + +// I need an index to go from image to associated mosaic (if existant) +// This code is copied from relastro/src/MosaicOps.c + +off_t getMosaicByTimes (unsigned int start, unsigned int stop, unsigned int *startMos, unsigned int *stopMos, off_t *indexMos) { + + // use bisection to find the overlapping mosaic + + off_t Nlo, Nhi, N; + + // find the last mosaic before start + Nlo = 0; Nhi = Nmosaic; + while (Nhi - Nlo > 10) { + N = 0.5*(Nlo + Nhi); + if (startMos[N] < start) { + Nlo = MAX(N, 0); + } else { + Nhi = MIN(N + 1, Nmosaic); + } + } + + // check for the matched mosaic starting from Nlo + // we may have to go much beyond Nlo since stop is not sorted + // can we use a sorted version of stop to check when we are beyond the valid range?? + for (N = Nlo; N < Nmosaic; N++) { + if (stop < stopMos[N]) continue; + if (start > startMos[N]) continue; + if (stop < startMos[N]) return (-1); + return (indexMos[N]); + } + + return (-1); +} + +// sort two times vectors and an index by first time vector +void sort_mosaic_times (unsigned int *S, unsigned int *E, off_t *I, off_t N) { + +# define SWAPFUNC(A,B){ unsigned int tmp_t; off_t tmp_i; \ + tmp_t = S[A]; S[A] = S[B]; S[B] = tmp_t; \ + tmp_t = E[A]; E[A] = E[B]; E[B] = tmp_t; \ + tmp_i = I[A]; I[A] = I[B]; I[B] = tmp_i; \ +} +# define COMPARE(A,B)(S[A] < S[B]) + + OHANA_SORT (N, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + +} + +// use the time to make the match, but use bracketing to make it faster: + +// find mosaic frames (unique time periods & photcode name matches mosaic) +void initMosaics (ImageSubset *image, off_t Nimage) { + + off_t i, Nmos, NMOSAIC; + unsigned int start, stop; + + unsigned int *startMos, *stopMos; + off_t *indexMos; + + Nmosaic = 0; + NMOSAIC = 10; + ALLOCATE (mosaic, Mosaic, NMOSAIC); + + ALLOCATE (startMos, unsigned int, NMOSAIC); + ALLOCATE (stopMos, unsigned int, NMOSAIC); + ALLOCATE (indexMos, off_t, NMOSAIC); + + /* find the mosaic images (coords.ctype = DIS); generate list of unique mosaics */ + for (i = 0; i < Nimage; i++) { + if (strcmp(&image[i].coords.ctype[4], "-DIS")) continue; + + /* set image time range */ + start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1); + stop = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1); + + /* a new mosaic, define ranges */ + mosaic[Nmosaic].start = start; + mosaic[Nmosaic].stop = stop; + mosaic[Nmosaic].coords = image[i].coords; + + startMos[Nmosaic] = start; + stopMos[Nmosaic] = stop; + indexMos[Nmosaic] = Nmosaic; + + Nmosaic ++; + if (Nmosaic == NMOSAIC) { + NMOSAIC += 10; + REALLOCATE (mosaic, Mosaic, NMOSAIC); + REALLOCATE (startMos, unsigned int, NMOSAIC); + REALLOCATE (stopMos, unsigned int, NMOSAIC); + REALLOCATE (indexMos, off_t, NMOSAIC); + } + } + + // sort the index, start, and stop by the start times: + sort_mosaic_times (startMos, stopMos, indexMos, Nmosaic); + + // array to store image->mosaic index + Nmosaic_for_images = Nimage; + ALLOCATE (mosaic_for_images, off_t, Nmosaic_for_images); + + /* now assign the WRP images to these mosaics */ + for (i = 0; i < Nimage; i++) { + mosaic_for_images[i] = -1; // default value for no mosaic found + + if (strcmp(&image[i].coords.ctype[4], "-WRP")) continue; + + /* set image time range */ + start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1); + stop = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1); + + Nmos = getMosaicByTimes (start, stop, startMos, stopMos, indexMos); + if (Nmos == -1) { + fprintf (stderr, "cannot match mosaic for %d (ID = %d)\n", (int) i, (int) image[i].imageID); + continue; + } + + // mosaic corresponding to this image + mosaic_for_images[i] = Nmos; + } + + free (startMos); + free (stopMos); + free (indexMos); + return; +} + +Mosaic *getMosaicForImage (off_t im) { + + off_t mos; + + if (im < 0) abort(); + if (im >= Nmosaic_for_images) abort(); + + // search for the mosaic that + mos = mosaic_for_images[im]; + if (mos < 0) return NULL; + + return &mosaic[mos]; +} + Index: /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/update_catalog_fiximids.c =================================================================== --- /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/update_catalog_fiximids.c (revision 35591) +++ /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/update_catalog_fiximids.c (revision 35591) @@ -0,0 +1,111 @@ +# include "fiximids.h" + +void update_catalog_fiximids (Catalog *catalog, ImageSubset *image, off_t *index, off_t Nimage) { + + float posAngle, pltScale; + off_t i, j, found; + + found = 0; + for (i = 0; i < catalog[0].Naverage; i++) { + + off_t m = catalog[0].average[i].measureOffset; + for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) { + off_t ID = catalog[0].measure[m].imageID; + if (ID <= 0) continue; // detections with imageID == 0 do not have a valid image (eg, ref photcode) + + // XXX only do GPC1 data for now + catalog[0].measure[m].pltscale = pltScale; + if (catalog[0].measure[m].photcode < 10000) continue; + if (catalog[0].measure[m].photcode > 10600) continue; + + // index[ID] = seqN + + off_t seq = index[ID]; + + Mosaic *mosaic = getMosaicForImage(seq); + Coords *coords = &image[seq].coords; + Coords *mosaicCoords = mosaic ? &mosaic->coords : NULL; + + double Xccd = catalog[0].measure[m].Xccd; + double Yccd = catalog[0].measure[m].Yccd; + + fiximids_local_astrometry (&posAngle, &pltScale, Xccd, Yccd, mosaicCoords, coords); + + catalog[0].measure[m].posangle = ToShortDegrees(posAngle); + catalog[0].measure[m].pltscale = pltScale; + myAssert(isfinite(catalog[0].measure[m].posangle), "oops: fiximids made a nan"); + myAssert(isfinite(catalog[0].measure[m].pltscale), "oops: fiximids made a nan"); + found ++; + } + } + + if (found) { + fprintf (stderr, "found "OFF_T_FMT" matches\n", found); + } +} + +// this is basically a re-write / adaptation of pmSourceLocalAstrometry +// posangle in degrees, plate scale in arcseconds/pixel +int fiximids_local_astrometry (float *posAngle, float *pltScale, double x, double y, Coords *mosaic, Coords *coords) { + + double Lx, Mx, Po, Qo, Px, Qx, Py, Qy; + + // calculate the astrometry for the coordinate of interest + XY_to_LM (&Lx, &Mx, x, y, coords); + if (mosaic) { + XY_to_LM (&Po, &Qo, Lx, Mx, mosaic); + } else { + Po = Lx; + Qo = Mx; + } + + XY_to_LM (&Lx, &Mx, x + 1.0, y, coords); + if (mosaic) { + XY_to_LM (&Px, &Qx, Lx, Mx, mosaic); + } else { + Px = Lx; + Qx = Mx; + } + + XY_to_LM (&Lx, &Mx, x, y + 1.0, coords); + if (mosaic) { + XY_to_LM (&Py, &Qy, Lx, Mx, mosaic); + } else { + Py = Lx; + Qy = Mx; + } + + // XXX units for the resulting Tangent Plane coordinates?? + + double dPdX = Px - Po; + double dPdY = Py - Po; + + double dQdX = Qx - Qo; + double dQdY = Qy - Qo; + + double pltScale_x = hypot(dPdX, dQdX); + double pltScale_y = hypot(dPdY, dQdY); + *pltScale = 3600.0*0.5*(pltScale_x + pltScale_y); + + double posAngle_x, posAngle_y; + double crossProduct = dPdX * dQdY - dPdY * dQdX; + if (crossProduct > 0.) { + *pltScale *= -1.0; + posAngle_x = atan2 (dQdX, dPdX); + posAngle_y = atan2 (dQdY, dPdY) - M_PI_2; + } else { + posAngle_x = atan2 (dQdX, -dPdX); + posAngle_y = atan2 (dQdY, dPdY) - M_PI_2; + } + + // with errors, these may end up on opposite sides of the M_PI boundary. + if (posAngle_x - posAngle_y > M_PI) { + posAngle_y += 2.0 * M_PI; + } + if (posAngle_y - posAngle_x > M_PI) { + posAngle_x += 2.0 * M_PI; + } + *posAngle = 0.5*(posAngle_x + posAngle_y)*DEG_RAD; + + return TRUE; +} Index: /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/update_dvo_fiximids.c =================================================================== --- /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/update_dvo_fiximids.c (revision 35591) +++ /branches/eam_branches/ipp-20130509/Ohana/src/uniphot/src/update_dvo_fiximids.c (revision 35591) @@ -0,0 +1,157 @@ +# include "fiximids.h" + +int update_dvo_fiximids (ImageSubset *image, off_t Nimage) { + + SkyTable *sky = NULL; + SkyList *skylist = NULL; + Catalog catalog; + off_t i, maxID, *index; + + // load the current sky table (layout of all SkyRegions) + sky = SkyTableLoadOptimal (CATDIR, NULL, NULL, FALSE, -1, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + + if (PARALLEL && !HOST_ID) { + update_dvo_fiximids_parallel (sky, image, Nimage); + return TRUE; + } + + // create an index for the image IDs + maxID = 0; + for (i = 0; i < Nimage; i++) { + maxID = MAX(maxID, image[i].imageID); + } + ALLOCATE (index, off_t, maxID + 1); + for (i = 0; i < maxID + 1; i++) { + index[i] = -1; + } + for (i = 0; i < Nimage; i++) { + if (!image[i].imageID) continue; // images with ID == 0 are virtual + index[image[i].imageID] = i; + } + + // set up the links between images and mosaics + initMosaics (image, Nimage); + + // determine the populated SkyRegions overlapping the requested area (default depth) + if (SINGLE_CPT) { + skylist = SkyRegionByCPT (sky, SINGLE_CPT); + } else { + skylist = SkyListByPatch (sky, -1, &UserPatch); + } + myAssert (skylist, "ooops!"); + + // update measurements for each populated catalog + for (i = 0; i < skylist[0].Nregions; i++) { + + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist[0].regions[i], HOST_ID)) continue; + + char localFilename[1024]; + snprintf (localFilename, 1024, "%s/%s.cpt", HOSTDIR, skylist->regions[i]->name); + + // set up the basic catalog info + catalog.filename = HOST_ID ? localFilename : skylist[0].filename[i]; + catalog.catformat = DVO_FORMAT_UNDEF; // not needed since we skip empty catalogs + catalog.catmode = DVO_MODE_UNDEF; // not needed since we skip empty catalogs + catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; + catalog.Nsecfilt = GetPhotcodeNsecfilt (); + + if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) { + fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename); + exit (1); + } + if (!catalog.Naves_disk) { + if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", catalog.filename); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + continue; + } + + update_catalog_fiximids (&catalog, image, index, Nimage); + + if (!UPDATE) { + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + continue; + } + + if (VERBOSE) fprintf (stderr, "saving catalog %s\n", catalog.filename); + + dvo_catalog_save (&catalog, VERBOSE); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + } + return (TRUE); +} + +# define DEBUG 1 + +int update_dvo_fiximids_parallel (SkyTable *sky, ImageSubset *image, off_t Nimage) { + + // write out the subset table of image information + char imageFile[512]; + snprintf (imageFile, 512, "%s/Images.subset.dat", CATDIR); + + if (!ImageSubsetSave (imageFile, image, Nimage)) { + fprintf (stderr, "failed to write image subset\n"); + exit (1); + } + + // now launch the fiximids_client jobs to the parallel hosts + + // load the list of hosts + HostTable *table = HostTableLoad (CATDIR, sky->hosts); + + int i; + for (i = 0; i < table->Nhosts; i++) { + + // ensure that the paths are absolute path names + char *tmppath = abspath (table->hosts[i].pathname, DVO_MAX_PATH); + free (table->hosts[i].pathname); + table->hosts[i].pathname = tmppath; + + char command[1024]; + snprintf (command, 1024, "fiximids_client -hostID %d -catdir %s -hostdir %s -images %s -region %f %f %f %f", + table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, imageFile, + UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); + + fprintf (stderr, "command: %s\n", command); + + char tmpline[1024]; + if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } + if (UPDATE) { snprintf (tmpline, 1024, "%s -update", command); strcpy (command, tmpline); } + + if (PARALLEL_MANUAL) continue; + + if (PARALLEL_SERIAL) { + int status = system (command); + if (status) { + fprintf (stderr, "ERROR running fiximids_client\n"); + exit (2); + } + } else { + // launch the job on the remote machine (no handshake) + int errorInfo = 0; + int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE); + if (!pid) { + if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo); + exit (1); + } + table->hosts[i].pid = pid; // save for future reference + } + } + + // wait for the remote jobs to be completed + if (PARALLEL_MANUAL) { + fprintf (stderr, "run the fiximids_client commands above. when these are done, hit return\n"); + getchar(); + } + if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) { + int status = HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + if (!status) return FALSE; + } + + return (TRUE); +} +