Index: branches/simmosaic_branches/Ohana/src/relastro/Makefile =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/Makefile (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/Makefile (revision 27839) @@ -56,4 +56,6 @@ $(SRC)/save_catalogs.$(ARCH).o \ $(SRC)/write_coords.$(ARCH).o \ +$(SRC)/CoordOps.$(ARCH).o \ +$(SRC)/FixProblemImages.$(ARCH).o \ $(SRC)/relastroVisual.$(ARCH).o Index: branches/simmosaic_branches/Ohana/src/relastro/doc/notes.txt =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/doc/notes.txt (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/doc/notes.txt (revision 27839) @@ -1,2 +1,61 @@ + +2010.03.24 + + I have finished only part of the optimizations listed below: I have + made the detection->image and image->mosaic lookups use bracket + searches. I have NOT converted the image->mosaic lookups to use the + parentID, so these matches are still based on time ranges. I have + also added a bit of positional constraint to the image match + function in 'select_images' (comparing Rmax of the sky regions to + Rmin of the images). These have had a huge improvement. Here is + the execution time for the SAS data (~450 images) with the old vs + the new code on pikake: + + old : 39min + new : 7min + + The bulk of the time for the new code is now disk I/O. + +2010.03.23 + + I've been working on optimization. I've updated ImageOps to make + the measure->image relationship generation use the image IDs along + with a bracket search -- this goes much faster as a result. More + work is needed: + + * select_images needs to apply astrometry to all image corners: + -> cache the image center and radius and use this to narrow down + the searches. + + * MosaicOps uses the time to match mosaics. I've added sorting to + speed this up, but this should be done based on the image and + parent index. + + -> make a tool to create the parent ID for existing DBs. + -> make sure addstar is populating the parent IDs + -> use the parent IDs to make the link. + +load image data: 0.006347 sec + setup sky: 0.091475 sec + convert image table: 0.091498 sec + select images: 2.925978 sec + init images: 2.926957 sec + init mosaics: 2.927319 sec +load images: 2.933762 sec +load catalog data: 9.772259 sec +make image bins: 9.780618 sec +set up image indexes: 10.196519 sec + +load image data: 0.006746 sec + setup sky: 0.091440 sec + convert image table: 0.091467 sec + select images: 0.374193 sec + init images: 0.375355 sec + init mosaics: 0.375753 sec +load images: 0.382606 sec +load catalog data: 7.260094 sec +make image bins: 7.268564 sec +set up image indexes: 7.687485 sec + 2008.03.01 Index: branches/simmosaic_branches/Ohana/src/relastro/include/relastro.h =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/include/relastro.h (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/include/relastro.h (revision 27839) @@ -91,6 +91,8 @@ int VERBOSE; +int VERBOSE2; int RESET; +int NLOOP; int UPDATE; int PLOTSTUFF; @@ -104,4 +106,7 @@ int CHIPORDER; +int UserCatalog; +double UserCatalogRA, UserCatalogDEC; + char *PHOTCODE_KEEP_LIST, *PHOTCODE_SKIP_LIST; int NphotcodesKeep, NphotcodesSkip; @@ -140,5 +145,5 @@ char *GetPhotnamebyCode PROTO((PhotCodeData *photcodes, int code)); void InterpolateGrid PROTO((float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords)); -int *SelectRefMosaic PROTO((Mosaic **refmosaic, int *Nimage)); +off_t *SelectRefMosaic PROTO((Mosaic **refmosaic, off_t *Nimage)); int args PROTO((int argc, char **argv)); int bcatalog PROTO((Catalog *subcatalog, Catalog *catalog)); @@ -151,14 +156,14 @@ void dump_grid PROTO((void)); int edge_check PROTO((double *x1, double *y1, double *x2, double *y2)); -void findImages PROTO((Catalog *catalog, int Ncatalog)); +void findImages PROTO((Catalog *catalog, int Ncatalog, int MATCHCAT)); int findMosaics PROTO((Catalog *catalog, int Ncatalog)); -Image *find_images PROTO((FITS_DB *db, GSCRegion *region, int Nregion, int *Nimage, int **LineNum)); +Image *find_images PROTO((FITS_DB *db, GSCRegion *region, off_t Nregion, off_t *Nimage, off_t **LineNum)); void set_db (FITS_DB *in); -int Shutdown (char *format, ...); +int Shutdown (char *format, ...) OHANA_FORMAT(printf, 1, 2); void TrapSignal (int sig); void SetProtect (int mode); int SetSignals (void); -GSCRegion *find_regions PROTO((Image *image, int Nimage, int *Nregions, GSCRegion *fullregion)); +GSCRegion *find_regions PROTO((Image *image, off_t Nimage, int *Nregions, GSCRegion *fullregion)); void freeGridBins PROTO((int Ncatalog)); void freeImageBins PROTO((int Ncatalog)); @@ -166,21 +171,21 @@ void free_catalogs PROTO((Catalog *catalog, int Ncatalog)); int gcatalog PROTO((Catalog *catalog, int FINAL)); -Coords *getCoords PROTO((int meas, int cat)); -float getMcal PROTO((int meas, int cat)); -float getMgrid PROTO((int meas, int cat)); -float getMmos PROTO((int meas, int cat)); -float getMrel PROTO((Catalog *catalog, int meas, int cat)); -GSCRegion *get_regions PROTO((double minRa, double maxRa, double minDec, double maxDec, int *Nregions)); -void getfullregion PROTO((Image *image, int Nimage, GSCRegion *fullregion)); -Image *getimage PROTO((int N)); -Image *getimages PROTO((int *N)); +// Coords *getCoords PROTO((off_t meas, int cat)); +float getMcal PROTO((off_t meas, int cat)); +float getMgrid PROTO((off_t meas, int cat)); +float getMmos PROTO((off_t meas, int cat)); +float getMrel PROTO((Catalog *catalog, off_t meas, int cat)); +GSCRegion *get_regions PROTO((double minRa, double maxRa, double minDec, double maxDec, off_t *Nregions)); +void getfullregion PROTO((Image *image, off_t Nimage, GSCRegion *fullregion)); +Image *getimage PROTO((off_t N)); +Image *getimages PROTO((off_t *N)); void global_stats PROTO((Catalog *catalog, int Ncatalog)); void initGrid PROTO((int dX, int dY)); void initGridBins PROTO((Catalog *catalog, int Ncatalog)); -void initImageBins PROTO((Catalog *catalog, int Ncatalog)); -void initImages PROTO((Image *input, int N)); +void initImageBins PROTO((Catalog *catalog, int Ncatalog, int FULLINIT)); +void initImages PROTO((Image *input, off_t N)); void initMosaicBins PROTO((Catalog *catalog, int Ncatalog)); -void initMosaicGrid PROTO((Image *image, int Nimage)); -void initMosaics PROTO((Image *image, int Nimage)); +void initMosaicGrid PROTO((Image *image, off_t Nimage)); +void initMosaics PROTO((Image *image, off_t Nimage)); void initMrel PROTO((Catalog *catalog, int Ncatalog)); void initialize PROTO((int argc, char **argv)); @@ -189,5 +194,5 @@ Catalog *load_catalogs PROTO((SkyList *skylist, int *Ncatalog, int subselect)); SkyList *load_images PROTO((FITS_DB *db, SkyRegion *region)); -Image *select_images PROTO((SkyList *skylist, Image *timage, int Ntimage, int **LineNumber, int *Nimage)); +Image *select_images PROTO((SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage)); void check_permissions (char *basefile); @@ -198,8 +203,8 @@ int main PROTO((int argc, char **argv)); -void mark_images PROTO((Image *image, int Nimage, Image *timage, int Ntimage)); -void matchImage PROTO((Catalog *catalog, int meas, int cat)); -void matchMosaics PROTO((Catalog *catalog, int meas, int cat)); -GSCRegion *name_region PROTO((char *name, int *Nregions)); +void mark_images PROTO((Image *image, off_t Nimage, Image *timage, off_t Ntimage)); +void matchImage PROTO((Catalog *catalog, off_t meas, int cat, int MATCHCAT)); +void matchMosaics PROTO((Catalog *catalog, off_t meas, int cat)); +GSCRegion *name_region PROTO((char *name, off_t *Nregions)); double opening_angle PROTO((double x1, double y1, double x2, double y2, double x3, double y3)); void plot_chisq PROTO((Catalog *catalog, int Ncatalog)); @@ -224,5 +229,5 @@ int setMrelOutput PROTO((Catalog *catalog, int Ncatalog, int mark)); void set_ZP PROTO((double ZERO)); -int setrefcode PROTO((Image *image, int Nimage)); +int setrefcode PROTO((Image *image, off_t Nimage)); void skip_measurements PROTO((Catalog *catalog, int pass)); void sortA PROTO((double *X, int N)); @@ -250,5 +255,5 @@ void fit_free (CoordFit *fit); void fit_add (CoordFit *fit, double x1, double y1, double x2, double y2, double wt); -void fit_eval (CoordFit *fit); +int fit_eval (CoordFit *fit); void fit_apply (CoordFit *fit, double *x2, double *y2, double x1, double y1); double **poly2d_dx (double **poly, int Nx, int Ny); @@ -256,8 +261,8 @@ double **poly2d_copy (double **poly, int Nx, int Ny); double poly2d_eval (double **poly, int Nx, int Ny, double x, double y); -CoordFit *fit_apply_coords (CoordFit *fit, Coords *coords); +int fit_apply_coords (CoordFit *fit, Coords *coords); int CoordsGetCenter (CoordFit *fit, double tol, double *xo, double *yo); CoordFit *CoordsSetCenter (CoordFit *input, double Xo, double Yo); -void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords); +int FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords); void FitMosaic (StarData *raw, StarData *ref, int Nmatch, Coords *coords); void FitSimple (StarData *raw, StarData *ref, int Nmatch, Coords *coords); @@ -268,24 +273,24 @@ int UpdateMosaic (Catalog *catalog, int Ncatalog); int UpdateMeasures (Catalog *catalog, int Ncatalog); -void fixImageRaw (Catalog *catalog, int Ncatalog, int im); +void fixImageRaw (Catalog *catalog, int Ncatalog, off_t im); void FlagOutliers(Catalog *catalog); int MeasFilterTest(Measure *measure); int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon); -int ParFactor (double *pR, double *pD, double R, double D, time_t T); +int ParFactor (double *pR, double *pD, double R, double D, double T, double Tmean); int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts); int FitPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *pR, double *pD, int Npts); int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts); -Mosaic *getMosaicForImage (int N); - -StarData *getImageRef (Catalog *catalog, int Ncatalog, int im, int *Nstars, CoordMode mode); -StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, CoordMode mode); - -Mosaic *getmosaics (int *N); -void initMosaics (Image *image, int Nimage); -StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, int mos, int *Nstars); -StarData *getMosaicRef (Catalog *catalog, int Ncatalog, int mos, int *Nstars); -Mosaic *getMosaicForImage (int im); +Mosaic *getMosaicForImage (off_t N); + +StarData *getImageRef (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode); +StarData *getImageRaw (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode); + +Mosaic *getmosaics (off_t *N); +void initMosaics (Image *image, off_t Nimage); +StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars); +StarData *getMosaicRef (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars); +Mosaic *getMosaicForImage (off_t im); double getMeanR (Measure *measure, Average *average, SecFilt *secfilt); @@ -297,2 +302,21 @@ int relastro_objects (void); int UpdateObjectOffsets (SkyList *skylist); + +int relastroVisualPlotRawRef(StarData *raw, StarData *ref, double dRmax, int numObj); +int relastroVisualPlotScatter(double values[], double thresh, int npts); +int relastroVisualPlotOutliers(Catalog *catalog, int offset, int Nmeasure, + StatType statsR, StatType statsD, double thresh); + + + +int FixProblemImages (SkyList *skylist); +int *getCatlist (int *N, off_t im); + +void initCoords (void); +void getOffsets (double *dPos, off_t *nPos, off_t N); +void saveOffsets (double dPos, off_t nPos, off_t N); +void setBadCoords (off_t N); +int badCoords (off_t N); +Coords *getCoords (off_t N); +int saveCoords (Coords *coords, off_t N); +void resetImageRaw (Catalog *catalog, int Ncatalog, off_t im); Index: branches/simmosaic_branches/Ohana/src/relastro/src/ConfigInit.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/ConfigInit.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/ConfigInit.c (revision 27839) @@ -60,5 +60,5 @@ /* update master photcode table if not defined */ sprintf (CatdirPhotcodeFile, "%s/Photcodes.dat", CATDIR); - if (!LoadPhotcodes (CatdirPhotcodeFile, MasterPhotcodeFile)) { + if (!LoadPhotcodes (CatdirPhotcodeFile, MasterPhotcodeFile, TRUE)) { fprintf (stderr, "error loading photcode table %s or master file %s\n", CatdirPhotcodeFile, MasterPhotcodeFile); exit (1); Index: branches/simmosaic_branches/Ohana/src/relastro/src/CoordOps.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/CoordOps.c (revision 27839) +++ branches/simmosaic_branches/Ohana/src/relastro/src/CoordOps.c (revision 27839) @@ -0,0 +1,83 @@ +# include "relastro.h" + +static double *dPosSum = NULL; // sum of dPos^2 for all measures on each image +static off_t *nPosSum = NULL; // sum of measures on each image (used for dPosSum) +static int *isBadCoord = NULL; // keep or reject each image? + +static Coords *oldCoords; // list of available images +static off_t NoldCoords; // number of available images + +void initCoords (void) { + + off_t N; + Image *images; + + images = getimages (&N); + + NoldCoords = N; + ALLOCATE (oldCoords, Coords, NoldCoords); + ALLOCATE (dPosSum, double, NoldCoords); + ALLOCATE (nPosSum, off_t, NoldCoords); + ALLOCATE (isBadCoord, int, NoldCoords); + memset (oldCoords, 0, N*sizeof(Coords)); + memset (dPosSum, 0, N*sizeof(double)); + memset (nPosSum, 0, N*sizeof(off_t)); + memset (isBadCoord, 0, N*sizeof(int)); +} + +int saveCoords (Coords *coords, off_t N) { + + if (N < 0) return FALSE; + if (N >= NoldCoords) return FALSE; + + memcpy (&oldCoords[N], coords, sizeof(Coords)); + return TRUE; +} + +Coords *getCoords (off_t N) { + + if (N < 0) return NULL; + if (N >= NoldCoords) return NULL; + + return (&oldCoords[N]); +} + +int badCoords (off_t N) { + + if (N < 0) return FALSE; + if (N >= NoldCoords) return FALSE; + + return (isBadCoord[N]); +} + +void setBadCoords (off_t N) { + + if (N < 0) return; + if (N >= NoldCoords) return; + + isBadCoord[N] = TRUE; + return; +} + +void saveOffsets (double dPos, off_t nPos, off_t N) { + + if (N < 0) return; + if (N >= NoldCoords) return; + + dPosSum[N] += dPos; + nPosSum[N] += nPos; + + return; +} + +void getOffsets (double *dPos, off_t *nPos, off_t N) { + + if (N < 0) return; + if (N >= NoldCoords) return; + + *dPos = dPosSum[N]; + *nPos = nPosSum[N]; + + return; +} + Index: branches/simmosaic_branches/Ohana/src/relastro/src/FitChip.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/FitChip.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/FitChip.c (revision 27839) @@ -17,5 +17,5 @@ // XXX save measurements of the fit quality (scatter, chisq) in the image table -void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords) { +int FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords) { int i, Nscatter, Niter, skip; @@ -84,5 +84,5 @@ default: fprintf (stderr, "invalid chip order %d\n", coords[0].Npolyterms); - abort (); + skip = TRUE; } if (skip) { @@ -90,12 +90,20 @@ fit_free (fit); free (values); - return; + return FALSE; } - fprintf (stderr, "scatter limit: %f based on %d detections; using %d of %d for fit\n", dRmax, Nscatter, fit[0].Npts, Nmatch); + // fprintf (stderr, "scatter limit: %f based on %d detections; using %d of %d for fit\n", dRmax, Nscatter, fit[0].Npts, Nmatch); // measure the fit, update the coords & object coordinates - fit_eval (fit); - fit_apply_coords (fit, coords); + if (!fit_eval (fit)) { + fprintf (stderr, "failed to fit new model\n"); + return FALSE; + } + + if (!fit_apply_coords (fit, coords)) { + fprintf (stderr, "failed to fit new model\n"); + return FALSE; + } + fit_free (fit); @@ -107,5 +115,5 @@ free (values); - return; + return TRUE; } Index: branches/simmosaic_branches/Ohana/src/relastro/src/FitSimple.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/FitSimple.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/FitSimple.c (revision 27839) @@ -17,5 +17,5 @@ // XXX the value needs to be set in a more intelligent way - if (dR > 0.15) continue; + if (dR > 0.50) continue; fit_add (fit, raw[i].X, raw[i].Y, ref[i].P, ref[i].Q, 1.0); Index: branches/simmosaic_branches/Ohana/src/relastro/src/FixProblemImages.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/FixProblemImages.c (revision 27839) +++ branches/simmosaic_branches/Ohana/src/relastro/src/FixProblemImages.c (revision 27839) @@ -0,0 +1,83 @@ +# include "relastro.h" + +// loop over all images. for any images which have a bad coordinate solution, replace the +// original coordinates and recalculate the positions + +int FixProblemImages (SkyList *skylist) { + + int Nbad; + off_t i, Nimage; + Image *image; + SkyList sublist; + + // allocate so we can reallocate below + ALLOCATE (sublist.regions, SkyRegion *, 1); + ALLOCATE (sublist.filename, char *, 1); + + image = getimages (&Nimage); + + Nbad = 0; + // first check on the dPos reported for each image + for (i = 0; i < Nimage; i++) { + double dPosSum, dPos; + off_t nPos; + + // check if this image should be fixed + if (badCoords(i)) { + Nbad ++; + continue; + } + + getOffsets (&dPosSum, &nPos, i); + dPos = sqrt(dPosSum / nPos); + if (dPos > 4.0) { + setBadCoords (i); + Nbad ++; + } + } + + fprintf (stderr, "fixing %d images\n", Nbad); + + for (i = 0; i < Nimage; i++) { + int j, cat, Ncat, *catlist, Ncatlist; + Catalog *catalog; + + // check if this image should be fixed + if (!badCoords(i)) continue; + + fprintf (stderr, "fixing %s\n", image[i].name); + + // I need a list of the catalogs for this image + catlist = getCatlist(&Ncatlist, i); + + // allocate Ncatlist skylist regions + REALLOCATE (sublist.regions, SkyRegion *, Ncatlist); + REALLOCATE (sublist.filename, char *, Ncatlist); + sublist.Nregions = Ncatlist; + sublist.ownElements = FALSE; // this list is only holding a view to the elements + + // copy the desired catalogs from skylist to skylistSubset + for (j = 0; j < Ncatlist; j++) { + cat = catlist[j]; + sublist.filename[j] = skylist[0].filename[cat]; + sublist.regions[j] = skylist[0].regions[cat]; + } + + catalog = load_catalogs (&sublist, &Ncat, FALSE); + assert (Ncat == Ncatlist); + + // match measurements with images + initImageBins (catalog, Ncat, FALSE); + findImages (catalog, Ncat, FALSE); + + // update the detection coordinates using the new image parameters + resetImageRaw (catalog, Ncat, i); + + freeImageBins (Ncat); + + // write the updated detections to disk + save_catalogs (catalog, Ncat); + } + + return (TRUE); +} Index: branches/simmosaic_branches/Ohana/src/relastro/src/GetAstromError.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/GetAstromError.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/GetAstromError.c (revision 27839) @@ -8,11 +8,11 @@ switch (mode) { case ERROR_MODE_RA: - dPobs = measure[0].dXccd; // need to redefine this as RAerr + dPobs = measure[0].dXccd / 100.0; // need to redefine this as RAerr break; case ERROR_MODE_DEC: - dPobs = measure[0].dYccd; // need to redefine this as RAerr + dPobs = measure[0].dYccd / 100.0; // need to redefine this as RAerr break; case ERROR_MODE_POS: - dPobs = hypot (measure[0].dXccd, measure[0].dYccd); // need to redefine this as RAerr + dPobs = hypot (measure[0].dXccd, measure[0].dYccd) / 100.0; // need to redefine this as RAerr break; default: Index: branches/simmosaic_branches/Ohana/src/relastro/src/ImageOps.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/ImageOps.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/ImageOps.c (revision 27839) @@ -1,35 +1,64 @@ # include "relastro.h" -#define TESTING - + +# define USE_IMAGE_ID 1 + +static off_t **bin; // link from catalog,measure to image +static int **clist; // catalog which supplied measurement on image +static off_t **mlist; // measure reference for measurement on image +static off_t *Nlist; // number of measurements on image +static off_t *NLIST; // allocated number of measurements on image + +static Image *image; // list of available images +static off_t Nimage; // number of available images + +static int *Ncatlist; // catalogs associated with each image +static int *NCATLIST; // catalogs associated with each image +static int **catlist; // catalogs associated with each image + +// if we search by image ID, we sort (imageIDs, imageIdx) by imageIDs to get a sorted +// index + +# if USE_IMAGE_ID +static off_t *imageIDs; // list of all image IDs +static off_t *imageIdx; // list of index for image IDs +# else static unsigned int *start; static unsigned int *stop; -static int **bin; - -static int **clist; -static int **mlist; -static int *Nlist; -static int *NLIST; - -static Image *image; -static int Nimage; - -Image *getimages (int *N) { +# endif + +Image *getimages (off_t *N) { *N = Nimage; return (image); } -Image *getimage (int N) { +Image *getimage (off_t N) { return (&image[N]); } -void initImages (Image *input, int N) { - - int i; +int *getCatlist (int *N, off_t im) { + + *N = Ncatlist[im]; + return (catlist[im]); +} + +void initImages (Image *input, off_t N) { + + off_t i; image = input; Nimage = N; - ALLOCATE (start, unsigned, Nimage); - ALLOCATE (stop, unsigned, Nimage); +# if USE_IMAGE_ID + ALLOCATE (imageIDs, off_t, Nimage); + ALLOCATE (imageIdx, off_t, Nimage); + + for (i = 0; i < Nimage; i++) { + imageIdx[i] = i; + imageIDs[i] = image[i].imageID; + } + llsortpair (imageIDs, imageIdx, Nimage); +# else + ALLOCATE (start, unsigned, Nimage); + ALLOCATE (stop, unsigned, Nimage); for (i = 0; i < Nimage; i++) { @@ -37,26 +66,67 @@ stop[i] = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1); } -} - -void initImageBins (Catalog *catalog, int Ncatalog) { - - int i, j; - - ALLOCATE (bin, int *, Ncatalog); +# endif +} + +off_t getImageByID (off_t ID) { + + // we have a pair of vectors (imageIDs, imageIdx) sorted by imageIDs + // use bisection to find the specified image ID + +# if USE_IMAGE_ID + off_t Nlo, Nhi, N; + + Nlo = 0; Nhi = Nimage; + while (Nhi - Nlo > 10) { + N = 0.5*(Nlo + Nhi); + if (imageIDs[N] < ID) { + Nlo = MAX(N, 0); + } else { + Nhi = MIN(N + 1, Nimage); + } + } + + for (N = Nlo; N < Nhi; N++) { + if (imageIDs[N] == ID) + return (imageIdx[N]); + } +# endif + + return (-1); +} + +// these are really image & catalog indexes +void initImageBins (Catalog *catalog, int Ncatalog, int FULLINIT) { + + off_t i, j; + + ALLOCATE (bin, off_t *, Ncatalog); for (i = 0; i < Ncatalog; i++) { - ALLOCATE (bin[i], int, MAX (catalog[i].Nmeasure, 1)); + ALLOCATE (bin[i], off_t, MAX (catalog[i].Nmeasure, 1)); for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1; } - ALLOCATE (Nlist, int, Nimage); - ALLOCATE (NLIST, int, Nimage); - ALLOCATE (clist, int *, Nimage); - ALLOCATE (mlist, int *, Nimage); + ALLOCATE (Nlist, off_t, Nimage); + ALLOCATE (NLIST, off_t, Nimage); + ALLOCATE (clist, int *, Nimage); + ALLOCATE (mlist, off_t *, Nimage); for (i = 0; i < Nimage; i++) { - Nlist[i] = 0; + Nlist[i] = 0; NLIST[i] = 100; ALLOCATE (clist[i], int, NLIST[i]); - ALLOCATE (mlist[i], int, NLIST[i]); + ALLOCATE (mlist[i], off_t, NLIST[i]); + } + + if (FULLINIT) { + ALLOCATE (Ncatlist, int, Nimage); + ALLOCATE (NCATLIST, int, Nimage); + ALLOCATE (catlist, int *, Nimage); + + for (i = 0; i < Nimage; i++) { + Ncatlist[i] = 0; + NCATLIST[i] = 32; + ALLOCATE (catlist[i], int, NCATLIST[i]); + } } } @@ -64,5 +134,5 @@ void freeImageBins (int Ncatalog) { - int i; + off_t i; for (i = 0; i < Ncatalog; i++) { @@ -79,26 +149,82 @@ /* match measurements to images */ -void findImages (Catalog *catalog, int Ncatalog) { - - int i, j; +void findImages (Catalog *catalog, int Ncatalog, int MATCHCAT) { + + off_t i, j; char *name; for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Nmeasure; j++) { - matchImage (catalog, j, i); - } - } - - for (i = 0; VERBOSE && (i < Nimage); i++) { + // XXX : update to optionally restrict by photcode equiv a la relphot + // ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[j].photcode); + // if (photcode[0].code != ecode) continue; + matchImage (catalog, j, i, MATCHCAT); + } + } + + for (i = 0; VERBOSE2 && (i < Nimage); i++) { name = GetPhotcodeNamebyCode (image[i].photcode); - fprintf (stderr, "image %d has %d measures (%s, %s)\n", i, Nlist[i], + fprintf (stderr, "image %lld has %lld measures (%s, %s)\n", (long long) i, (long long) Nlist[i], ohana_sec_to_date(image[i].tzero), name); } } -/* modify this function to use the measure->imageID field */ -void matchImage (Catalog *catalog, int meas, int cat) { - - int i; +# if USE_IMAGE_ID +// this is the imageID-based match +void matchImage (Catalog *catalog, off_t meas, int cat, int MATCHCAT) { + + off_t idx, ID; + Measure *measure; + int i, found; + + measure = &catalog[cat].measure[meas]; + + ID = measure[0].imageID; + idx = getImageByID (ID); + if (idx == -1) { + if (VERBOSE2) fprintf (stderr, "can't match detection to image?\n"); + return; + } + + // index for (catalog, measure) -> image + bin[cat][meas] = idx; + + // index for image, Nentry -> catalog + clist[idx][Nlist[idx]] = cat; + + // index for image, Nentry -> measure + mlist[idx][Nlist[idx]] = meas; + Nlist[idx] ++; + + if (Nlist[idx] == NLIST[idx]) { + NLIST[idx] += 100; + REALLOCATE (clist[idx], int, NLIST[idx]); + REALLOCATE (mlist[idx], off_t, NLIST[idx]); + } + + if (MATCHCAT) { + // index for image -> catalog list + found = FALSE; + for (i = 0; !found && (i < Ncatlist[idx]); i++) { + if (catlist[idx][i] == cat) found = TRUE; + } + if (!found) { + catlist[idx][Ncatlist[idx]] = cat; + Ncatlist[idx] ++; + if (Ncatlist[idx] == NCATLIST[idx]) { + NCATLIST[idx] += 32; + REALLOCATE (catlist[idx], int, NCATLIST[idx]); + } + } + } + + return; +} + +# else +// this is the time-based match +void matchImage (Catalog *catalog, off_t meas, int cat, int MATCHCAT) { + + off_t i; Measure *measure; @@ -107,4 +233,6 @@ /* find the image that supplied this measurement */ for (i = 0; i < Nimage; i++) { + // let's try the very slow method first before adding a bisection search + // if (image[i].imageID != measure[0].imageID) continue; if (image[0].photcode == -1) continue; if (measure[0].photcode != image[i].photcode) continue; @@ -124,14 +252,18 @@ if (Nlist[i] == NLIST[i]) { NLIST[i] += 100; - REALLOCATE (clist[i], int, NLIST[i]); - REALLOCATE (mlist[i], int, NLIST[i]); + REALLOCATE (clist[i], int, NLIST[i]); + REALLOCATE (mlist[i], off_t, NLIST[i]); } return; } -} - -Coords *getCoords (int meas, int cat) { - - int i; + if (VERBOSE2) fprintf (stderr, "can't match detection to image?\n"); + return; +} +# endif + +/* +Coords *getCoords (off_t meas, int cat) { + + off_t i; i = bin[cat][meas]; @@ -139,8 +271,9 @@ return (&image[i].coords); } +*/ void plot_images () { - int i, bin; + off_t i, bin; double *xlist, *Mlist, *dlist; Graphdata graphdata; @@ -188,18 +321,30 @@ // return StarData values for detections in the specified image, converting coordinates from the // chip positions: X,Y -> L,M -> P,Q -> R,D -void fixImageRaw (Catalog *catalog, int Ncatalog, int im) { - - int i, m, c, n; +void fixImageRaw (Catalog *catalog, int Ncatalog, off_t im) { + + off_t i, m, c, n, nPos; double X, Y, L, M, P, Q, R, D, dR, dD; + double dPos; Mosaic *mosaic; Coords *moscoords, *imcoords; + // check if this image is bad and should be skipped + if (badCoords(im)) return; + + // WRP images need to have an associated mosaic moscoords = NULL; - mosaic = getMosaicForImage (im); - if (mosaic != NULL) { - moscoords = &mosaic[0].coords; + if (!strcmp(&image[im].coords.ctype[4], "-WRP")) { + mosaic = getMosaicForImage (im); + if (mosaic == NULL) return; // if we cannot find the associated image, skip it + moscoords = &mosaic[0].coords; } imcoords = &image[im].coords; + + // accumulate the rms position offsets. if this value, or any specific entry, is too + // large, we will reset the image to the original coords at the end of the analysis + + dPos = 0.0; + nPos = 0; for (i = 0; i < Nlist[im]; i++) { @@ -226,14 +371,30 @@ dD = 3600.0*(catalog[c].average[n].D - D); - if (fabs(catalog[c].measure[m].dR - dR) > 10.0) { - // XXXXX running into this still for last megacam exposure: wrong mosaic? - // ???? inconsistently hitting this???? + // complain if the new location is far from the average location + if (fabs(dR) > 2.0) { fprintf (stderr, "!"); - // abort (); - } - if (fabs(catalog[c].measure[m].dD - dD) > 10.0) { + setBadCoords (im); // report a failure for this image + return; + } + if (fabs(dD) > 2.0) { fprintf (stderr, "*"); - // abort (); - } + setBadCoords (im); // report a failure for this image + return; + } + + // complain if the new location is far from the old location + if (fabs(catalog[c].measure[m].dR - dR) > 2.0) { + fprintf (stderr, "@"); + setBadCoords (im); // report a failure for this image + return; + } + if (fabs(catalog[c].measure[m].dD - dD) > 2.0) { + fprintf (stderr, "#"); + setBadCoords (im); // report a failure for this image + return; + } + + dPos += SQ(catalog[c].measure[m].dR - dR) + SQ(catalog[c].measure[m].dD - dD); + nPos ++; catalog[c].measure[m].dR = dR; @@ -251,4 +412,74 @@ } } + + saveOffsets (dPos, nPos, im); + + return; +} + +// return StarData values for detections in the specified image, converting coordinates from the +// chip positions: X,Y -> L,M -> P,Q -> R,D +void resetImageRaw (Catalog *catalog, int Ncatalog, off_t im) { + + off_t i, m, c, n; + double X, Y, L, M, P, Q, R, D, dR, dD; + + Mosaic *mosaic; + Coords *moscoords, *imcoords, *oldcoords; + + // check if this image is bad and should be skipped + if (!badCoords(im)) { + fprintf (stderr, "ERROR: inconsistent result?"); + abort(); + } + + // replace the current coords with the old coords: + oldcoords = getCoords (im); + memcpy (&image[im].coords, oldcoords, sizeof(Coords)); + + // WRP images need to have an associated mosaic + moscoords = NULL; + if (!strcmp(&image[im].coords.ctype[4], "-WRP")) { + mosaic = getMosaicForImage (im); + if (mosaic == NULL) return; // if we cannot find the associated image, skip it + moscoords = &mosaic[0].coords; + } + imcoords = &image[im].coords; + + for (i = 0; i < Nlist[im]; i++) { + m = mlist[im][i]; + c = clist[im][i]; + + X = catalog[c].measure[m].Xccd; + Y = catalog[c].measure[m].Yccd; + n = catalog[c].measure[m].averef; + + dR = dD = 0.0; + if (moscoords == NULL) { + // this is a Simple image (not a mosaic) + // note that for a Simple image, L,M = P,Q + XY_to_LM (&L, &M, X, Y, imcoords); + LM_to_RD (&R, &D, L, M, imcoords); + } else { + XY_to_LM (&L, &M, X, Y, imcoords); + XY_to_LM (&P, &Q, L, M, moscoords); + LM_to_RD (&R, &D, P, Q, moscoords); + } + + catalog[c].measure[m].dR = dR; + catalog[c].measure[m].dD = dD; + + if (catalog[c].measure[m].dR > +180.0*3600.0) { + // average on high end of boundary, move star up + R += 360.0; + catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R); + } + if (catalog[c].measure[m].dR < -180.0*3600.0) { + // average on low end of boundary, move star down + R -= 360.0; + catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R); + } + } + return; } @@ -257,7 +488,7 @@ // chip positions: X,Y -> L,M -> P,Q -> R,D. This function is used by the image fitting steps, for // which the detections have already been filtered when they were loaded (bcatalog) -StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, CoordMode mode) { - - int i, m, c, n; +StarData *getImageRaw (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode) { + + off_t i, m, c, n; Mosaic *mosaic; @@ -265,16 +496,16 @@ StarData *raw; - ALLOCATE (raw, StarData, Nlist[im]); - mosaic = NULL; moscoords = NULL; if (mode == MODE_MOSAIC) { - mosaic = getMosaicForImage (im); - if (mosaic == NULL) { - fprintf (stderr, "mosaic not found for image %s\n", image[im].name); - exit (1); - } - moscoords = &mosaic[0].coords; - } + mosaic = getMosaicForImage (im); + if (mosaic == NULL) { + fprintf (stderr, "mosaic not found for image %s\n", image[im].name); + return NULL; + } + moscoords = &mosaic[0].coords; + } + + ALLOCATE (raw, StarData, Nlist[im]); for (i = 0; i < Nlist[im]; i++) { @@ -319,7 +550,7 @@ LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, moscoords); break; - default: - fprintf (stderr, "error: invalid mode in getImageRaw"); - abort (); + default: + fprintf (stderr, "error: invalid mode in getImageRaw"); + abort (); } } @@ -332,13 +563,11 @@ // the sky positions: R,D -> P,Q -> L,M -> X,Y -StarData *getImageRef (Catalog *catalog, int Ncatalog, int im, int *Nstars, CoordMode mode) { - - int i, m, c, n; +StarData *getImageRef (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode) { + + off_t i, m, c, n; Mosaic *mosaic; Coords *moscoords; StarData *ref; - - ALLOCATE (ref, StarData, Nlist[im]); mosaic = NULL; @@ -348,8 +577,10 @@ if (mosaic == NULL) { fprintf (stderr, "mosaic not found for image %s\n", image[im].name); - exit (1); + return NULL; } moscoords = &mosaic[0].coords; } + + ALLOCATE (ref, StarData, Nlist[im]); for (i = 0; i < Nlist[im]; i++) { @@ -371,14 +602,14 @@ switch (mode) { case MODE_SIMPLE: - RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, &image[im].coords); - ref[i].L = ref[i].P; - ref[i].M = ref[i].Q; - LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords); - break; + RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, &image[im].coords); + ref[i].L = ref[i].P; + ref[i].M = ref[i].Q; + LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords); + break; case MODE_MOSAIC: - RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, moscoords); - LM_to_XY (&ref[i].L, &ref[i].M, ref[i].P, ref[i].Q, moscoords); - LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords); - break; + RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, moscoords); + LM_to_XY (&ref[i].L, &ref[i].M, ref[i].P, ref[i].Q, moscoords); + LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords); + break; default: fprintf (stderr, "invalid case"); @@ -399,5 +630,6 @@ return; - int i, j, k, m, N, Ndel, Nave, Nmax, TOOFEW, Nsecfilt; + int Ndel, Nave; + off_t i, j, k, m, N, Nmax, TOOFEW, Nsecfilt; double Ns, theta, x, y; double *R, *D, *dR, *dD; @@ -490,5 +722,6 @@ void FlagOutliers2D (Catalog *catalog) { - int i, j, k, m, N, Ndel, Nave, Nmax, TOOFEW, Nsecfilt; + int Ndel, Nave; + off_t i, j, k, m, N, Nmax, TOOFEW, Nsecfilt; double *index; double Ns, theta, x, y; @@ -569,5 +802,5 @@ // recalculate image center, sigma based on closest 50% of points for(k = 0; k < N; k++) { - int ind = (int) index[k]; + off_t ind = (off_t) index[k]; R[k] = catalog[0].measure[ind].dR; D[k] = catalog[0].measure[ind].dD; Index: branches/simmosaic_branches/Ohana/src/relastro/src/MosaicOps.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/MosaicOps.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/MosaicOps.c (revision 27839) @@ -2,26 +2,78 @@ // array of mosaic definition structures -static int Nmosaic; +static off_t Nmosaic; static Mosaic *mosaic; // list of all images associated with a mosaic -static int *Nmosaic_own_images; // number of images for this mosaic -static int *Amosaic_own_images; // size of allocated array -static int **mosaic_own_images; // array of arrays: mosaic -> images +static off_t *Nmosaic_own_images; // number of images for this mosaic +static off_t *Amosaic_own_images; // size of allocated array +static off_t **mosaic_own_images; // array of arrays: mosaic -> images // list of mosaic associated with each image -static int Nmosaic_for_images; // number of images (for internal checks) -static int *mosaic_for_images; // array of: image -> mosaic - -Mosaic *getmosaics (int *N) { +static off_t Nmosaic_for_images; // number of images (for off_ternal checks) +static off_t *mosaic_for_images; // array of: image -> mosaic + +Mosaic *getmosaics (off_t *N) { *N = Nmosaic; return (mosaic); } +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 + +} + +// first, let's continue to 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 (Image *image, int Nimage) { - - int i, j, found, NMOSAIC; +void initMosaics (Image *image, off_t Nimage) { + + off_t i, Nmos, NMOSAIC; unsigned int start, stop; + + unsigned int *startMos, *stopMos; + off_t *indexMos; Nmosaic = 0; @@ -29,7 +81,10 @@ ALLOCATE (mosaic, Mosaic, NMOSAIC); - ALLOCATE (Nmosaic_own_images, int, NMOSAIC); - ALLOCATE (Amosaic_own_images, int, NMOSAIC); - ALLOCATE (mosaic_own_images, int *, NMOSAIC); + ALLOCATE (Nmosaic_own_images, off_t, NMOSAIC); + ALLOCATE (Amosaic_own_images, off_t, NMOSAIC); + ALLOCATE (mosaic_own_images, off_t *, 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 */ @@ -54,5 +109,9 @@ Nmosaic_own_images[Nmosaic] = 0; Amosaic_own_images[Nmosaic] = 10; - ALLOCATE (mosaic_own_images[Nmosaic], int, Amosaic_own_images[Nmosaic]); + ALLOCATE (mosaic_own_images[Nmosaic], off_t, Amosaic_own_images[Nmosaic]); + + startMos[Nmosaic] = start; + stopMos[Nmosaic] = stop; + indexMos[Nmosaic] = Nmosaic; Nmosaic ++; @@ -60,13 +119,19 @@ NMOSAIC += 10; REALLOCATE (mosaic, Mosaic, NMOSAIC); - REALLOCATE (mosaic_own_images, int *, NMOSAIC); - REALLOCATE (Nmosaic_own_images, int, NMOSAIC); - REALLOCATE (Amosaic_own_images, int, NMOSAIC); - } - } - + REALLOCATE (mosaic_own_images, off_t *, NMOSAIC); + REALLOCATE (Nmosaic_own_images, off_t, NMOSAIC); + REALLOCATE (Amosaic_own_images, off_t, 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, int, Nmosaic_for_images); + ALLOCATE (mosaic_for_images, off_t, Nmosaic_for_images); /* now assign the WRP images to these mosaics */ @@ -80,28 +145,26 @@ stop = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1); - /* find existing mosaic with this time range */ - found = FALSE; - for (j = 0; !found && (j < Nmosaic); j++) { - if (stop < mosaic[j].start) continue; - if (start > mosaic[j].stop) continue; - found = TRUE; - break; - } - /* if no matching mosaic exists, skip this image */ - if (!found) continue; + Nmos = getMosaicByTimes (start, stop, startMos, stopMos, indexMos); + if (Nmos == -1) { + fprintf (stderr, "cannot match mosaic for %s\n", image[i].name); + continue; + } // mosaic corresponding to this image - mosaic_for_images[i] = j; + mosaic_for_images[i] = Nmos; // add image to mosaic_own_image list - mosaic_own_images[j][Nmosaic_own_images[j]] = i; - Nmosaic_own_images[j] ++; - if (Nmosaic_own_images[j] == Amosaic_own_images[j]) { - Amosaic_own_images[j] += 10; - REALLOCATE (mosaic_own_images[j], int, Amosaic_own_images[j]); - } - assert (Nmosaic_own_images[j] < Amosaic_own_images[j]); - } - + mosaic_own_images[Nmos][Nmosaic_own_images[Nmos]] = i; + Nmosaic_own_images[Nmos] ++; + if (Nmosaic_own_images[Nmos] == Amosaic_own_images[Nmos]) { + Amosaic_own_images[Nmos] += 10; + REALLOCATE (mosaic_own_images[Nmos], off_t, Amosaic_own_images[Nmos]); + } + assert (Nmosaic_own_images[Nmos] < Amosaic_own_images[Nmos]); + } + + free (startMos); + free (stopMos); + free (indexMos); return; } @@ -109,7 +172,7 @@ // return StarData values for detections in the specified image, converting coordinates from the // chip positions: X,Y -> L,M -> P,Q -> R,D -StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, int mos, int *Nstars) { - - int i, j, im, Nraw, Nnew; +StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars) { + + off_t i, j, im, Nraw, Nnew; StarData *raw, *new; @@ -125,4 +188,8 @@ // this function does the reverse-lookup for the mosaic corresponding to this image new = getImageRaw (catalog, Ncatalog, im, &Nnew, MODE_MOSAIC); + if (!new) { + fprintf (stderr, "inconsistent: missing mosaic for image already associated with a mosaic? (1)\n"); + abort(); + } // merge new and raw @@ -142,7 +209,7 @@ // return StarData values for averages positions in the specified image, converting coordinates from // the sky positions: R,D -> P,Q -> L,M -> X,Y -StarData *getMosaicRef (Catalog *catalog, int Ncatalog, int mos, int *Nstars) { - - int i, j, im, Nref, Nnew; +StarData *getMosaicRef (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars) { + + off_t i, j, im, Nref, Nnew; StarData *ref, *new; @@ -157,4 +224,8 @@ // this function does the reverse-lookup for the mosaic corresponding to this image new = getImageRef (catalog, Ncatalog, im, &Nnew, MODE_MOSAIC); + if (!new) { + fprintf (stderr, "inconsistent: missing mosaic for image already associated with a mosaic? (2)\n"); + abort(); + } // merge new and ref @@ -172,7 +243,7 @@ } -Mosaic *getMosaicForImage (int im) { - - int mos; +Mosaic *getMosaicForImage (off_t im) { + + off_t mos; if (im < 0) abort(); Index: branches/simmosaic_branches/Ohana/src/relastro/src/ParFactor.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/ParFactor.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/ParFactor.c (revision 27839) @@ -42,5 +42,5 @@ /* given RA, DEC, Time, calculate the parallax factor */ -int ParFactor (double *pR, double *pD, double R, double D, time_t T) { +int ParFactor (double *pR, double *pD, double R, double D, double T, double Tmean) { double jd; @@ -49,5 +49,5 @@ /* given a time T in UNIX seconds, determine the solar longitude S */ - jd = ohana_sec_to_jd (T); + jd = ohana_sec_to_jd (365.25*86400.0*(T + Tmean)); sun_ecliptic (jd, &L, &B, &E); Index: branches/simmosaic_branches/Ohana/src/relastro/src/UpdateChips.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/UpdateChips.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/UpdateChips.c (revision 27839) @@ -4,7 +4,8 @@ /* we can measure new image parameters for each non-mosaic chip independently */ - int i, Nimage, Nraw, Nref; + off_t i, Nimage, Nraw, Nref; Image *image; StarData *raw, *ref; + Coords *oldCoords; image = getimages (&Nimage); @@ -17,14 +18,22 @@ /* convert measure coordinates to raw entries */ raw = getImageRaw (catalog, Ncatalog, i, &Nraw, MODE_MOSAIC); + if (!raw) continue; /* convert average coordinates to ref entries */ ref = getImageRef (catalog, Ncatalog, i, &Nref, MODE_MOSAIC); + if (!ref) continue; // note that Nraw & Nref must be equal: if not, we made a programming error in one of these two functions. assert (Nraw == Nref); + saveCoords (&image[i].coords, i); + // FitChip does iterative, clipped fitting - fprintf (stderr, "image %d : Nstars: %d\n", i, Nraw); - FitChip (raw, ref, Nraw, &image[i].coords); + // fprintf (stderr, "image %lld : Nstars: %lld\n", (long long) i, (long long) Nraw); + if (!FitChip (raw, ref, Nraw, &image[i].coords)) { + fprintf (stderr, "reject fit for image %s (%lld) : Nstars: %lld\n", image[i].name, (long long) i, (long long) Nraw); + oldCoords = getCoords (i); + memcpy (&image[i].coords, oldCoords, sizeof(Coords)); + } free (raw); Index: branches/simmosaic_branches/Ohana/src/relastro/src/UpdateMeasures.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/UpdateMeasures.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/UpdateMeasures.c (revision 27839) @@ -3,5 +3,5 @@ int UpdateMeasures (Catalog *catalog, int Ncatalog) { - int i, Nimage; + off_t i, Nimage; Image *image; Index: branches/simmosaic_branches/Ohana/src/relastro/src/UpdateMosaic.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/UpdateMosaic.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/UpdateMosaic.c (revision 27839) @@ -4,5 +4,5 @@ /* we can measure new image parameters for each mosaic independently */ - int i, Nmosaic, Nstars; + off_t i, Nmosaic, Nstars; Mosaic *mosaic; StarData *raw, *ref; Index: branches/simmosaic_branches/Ohana/src/relastro/src/UpdateObjectOffsets.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 27839) @@ -1,3 +1,9 @@ # include "relastro.h" + +// We run through each DVO catalog, updating the measures that come from the modified images +// We need to watch for failures: +// * in UpdateMeasures, in fixImageRaw, we track the cumulative offset for each image +// * after all updates are done, we can check for any bad images and reset them to the +// original coordinates int UpdateObjectOffsets (SkyList *skylist) { @@ -28,6 +34,6 @@ // match measurements with images - initImageBins (&catalog, 1); - findImages (&catalog, 1); + initImageBins (&catalog, 1, FALSE); + findImages (&catalog, 1, FALSE); // update the detection coordinates using the new image parameters Index: branches/simmosaic_branches/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/UpdateObjects.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/UpdateObjects.c (revision 27839) @@ -1,5 +1,5 @@ # include "relastro.h" -static int Nmax; +static off_t Nmax; static double *X, *dX; static double *Y, *dY; @@ -13,5 +13,5 @@ void initObjectData (Catalog *catalog, int Ncatalog) { - int i, j; + off_t i, j; Nmax = 0; @@ -40,15 +40,19 @@ int UpdateObjects (Catalog *catalog, int Ncatalog) { - int i, j, k, m, N, Nsecfilt, found, kp; + off_t j, k, m; + int i, N, Nsecfilt, mode, result, status, XVERB; StatType statsR, statsD; Coords coords; - PMFit fit; - time_t To; - int mode, Nave, Npm, Npar, Nskip; - double Tmin, Tmax; - float mag; - int mask; - PhotCode *code; - + PMFit fitAve, fitPM, fitPAR, fit; + time_t T2000; + off_t Nave, Npm, Npar, Nskip; + off_t NaveSum, NpmSum, NparSum, NskipSum; + double Tmin, Tmax, Tmean, Trange; + + memset (&fit, 0, sizeof(fit)); + memset (&fitAve, 0, sizeof(fitAve)); + memset (&fitPM, 0, sizeof(fitPM)); + memset (&fitPAR, 0, sizeof(fitPAR)); + initObjectData (catalog, Ncatalog); @@ -64,7 +68,8 @@ strcpy (coords.ctype, "RA---SIN"); + XVERB = FALSE; + // use J2000 as a reference time - To = ohana_date_to_sec ("2000/01/01"); - Nave = Npar = Npm = 0; + T2000 = ohana_date_to_sec ("2000/01/01"); // XXX in the future, use catalog[0].Nsecfilt only? allow catalogs to have variable Nsecfilt? @@ -72,9 +77,10 @@ assert (catalog[0].Nsecfilt == Nsecfilt); + NaveSum = NparSum = NpmSum = NskipSum = 0; for (i = 0; i < Ncatalog; i++) { - if (VERBOSE) fprintf (stderr, "astrometrize catalog %d : %d ave, %d meas\n", i, catalog[i].Naverage, catalog[i].Nmeasure); - - Nskip = 0; + if (VERBOSE) fprintf (stderr, "astrometrize catalog %d : %lld ave, %lld meas\n", i, (long long) catalog[i].Naverage, (long long) catalog[i].Nmeasure); + + Nave = Npar = Npm = Nskip = 0; for (j = 0; j < catalog[i].Naverage; j++) { /* calculate the average value of R,D for a single star */ @@ -82,17 +88,19 @@ // skip objects which are known to be problematic // XXX include this code or not? - # if (0) +# if (0) if (catalog[i].average[j].code & STAR_BAD) { Nskip ++; continue; } - # endif +# endif N = 0; m = catalog[i].average[j].measureOffset; - Tmin = Tmax = (catalog[i].measure[m].t - To) / (86400*365.25); + Tmin = Tmax = (catalog[i].measure[m].t - T2000) / (86400*365.25); mode = FIT_MODE; + // find the basic properties of the detections for this object (Tmin, Tmax, Tmean) + Tmean = 0; for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { @@ -111,10 +119,10 @@ // exclude measurements by previous outlier detection // XXX include this code or not? - # if (0) +# if (0) if (catalog[i].measure[m].dbFlags & MEAS_BAD) { catalog[i].measure[m].dbFlags |= ID_MEAS_SKIP_ASTROM; continue; } - # endif +# endif catalog[i].measure[m].dbFlags |= ID_MEAS_USED_OBJ; @@ -122,12 +130,21 @@ R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); D[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); - T[N] = (catalog[i].measure[m].t - To) / (86400*365.25) ; // time relative to J2000 in years + T[N] = (catalog[i].measure[m].t - T2000) / (86400*365.25) ; // time relative to J2000 in years Tmin = MIN(Tmin, T[N]); Tmax = MAX(Tmax, T[N]); - - dR[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_RA); - dD[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_DEC); + Tmean += T[N]; + + // dX, dY : error in arcsec -- + dX[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_RA); + dY[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_DEC); dT[N] = catalog[i].measure[m].dt; + + // XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in + // arcseconds. dR, dD are the errors in those directions in degrees. IF we have + // non-circular errors (different values for X and Y), then dR and dD will be + // incorrect: they would need to be rotated to take out the position angle + dR[k] = dX[k] / 3600.0; + dD[k] = dY[k] / 3600.0; N++; @@ -141,5 +158,6 @@ // XXX add the parallax factor range as a criterion as well - if ((Tmax - Tmin) < PM_DT_MIN) mode = FIT_AVERAGE; + Trange = Tmax - Tmin; + if (Trange < PM_DT_MIN) mode = FIT_AVERAGE; if ((mode == FIT_PM_ONLY) && (N < PM_TOOFEW)) mode = FIT_AVERAGE; @@ -147,5 +165,9 @@ if (N < SRC_MEAS_TOOFEW) { // XXX need to define PHOTOM and ASTROM object flags + // XXX reset the average value fields? catalog[i].average[j].flags |= ID_STAR_FEW; + catalog[i].average[j].ChiSqAve = NAN; + catalog[i].average[j].ChiSqPM = NAN; + catalog[i].average[j].ChiSqPar = NAN; if (N < 2) continue; } @@ -154,103 +176,112 @@ coords.crval1 = R[0]; coords.crval2 = D[0]; + Tmean /= (float) N; - /* project all of the R,D coordinates to a plane centered on this coordinate */ - for (k = 0; k < N; k++) { - RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords); - dX[k] = dR[k]; - dY[k] = dD[k]; - // fprintf (stderr, "%d %f %f %f %f %f\n", k, T[k], R[k], D[k], X[k], Y[k]); + XVERB = FALSE && (catalog[i].measure[m].dM < 0.01) && (N == 6) && (mode == FIT_PM_ONLY); + + // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq + + if ((mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) { + // project all of the R,D coordinates to a plane centered on this coordinate set + // the times to be relative to Tmean (this is required for parallax as well) + for (k = 0; k < N; k++) { + RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords); + T[k] -= Tmean; + if (XVERB) { + fprintf (stderr, "%lld %f %f %f %f %f +/- %f %f\n", (long long) k, T[k], R[k], D[k], X[k], Y[k], dX[k], dY[k]); + } + } + + FitPM (&fitPM, X, dX, Y, dY, T, N); + if (XVERB) fprintf (stderr, "fitted: %f - %f : %f %f : %f %f : %f vs %f\n", Tmin, Tmax, fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.chisq, fitAve.chisq); + + // project Ro, Do back to RA,DEC + XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &coords); + if (XVERB) fprintf (stderr, "project: %f %f : %f %f : %f\n", fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p); + + fitPM.p = fitPM.dp = 0.0; + catalog[i].average[j].flags |= ID_STAR_FIT_PM; + Npm ++; + } + + if (mode == FIT_PM_AND_PAR) { + fprintf (stderr, "parallax fitting is still untested (%s, %d)", __FILE__, __LINE__); + exit (2); + + for (k = 0; k < N; k++) { + ParFactor (&pX[k], &pY[k], R[k], D[k], T[k], Tmean); + } + FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N); + XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords); + catalog[i].average[j].flags |= ID_STAR_FIT_PAR; + Npar ++; } - /* fit the model components as needed */ - switch (mode) { + // fit the average model + if ((mode == FIT_AVERAGE) || (mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) { + liststats (R, dR, N, &statsR); // WARNING: this function modifies R (do not use after here) + liststats (D, dD, N, &statsD); // WARNING: this function modifies D (do not use after here) + + fitAve.Ro = statsR.mean; + fitAve.dRo = 3600.0*statsR.sigma; + + fitAve.Do = statsD.mean; + fitAve.dDo = 3600.0*statsD.sigma; + + fitAve.chisq = 0.5*(statsR.chisq + statsD.chisq); + fitAve.Nfit = N; + + fitAve.uR = fitAve.duR = 0.0; + fitAve.uD = fitAve.duD = 0.0; + fitAve.p = fitAve.dp = 0.0; + catalog[i].average[j].flags |= ID_STAR_FIT_AVE; + Nave ++; + } + + /* choose the result based on the chisq values */ + // XXXX for now, just use the mode as the result: + result = mode; + + switch (result) { case FIT_AVERAGE: - liststats (R, dR, N, &statsR); - liststats (D, dD, N, &statsD); - - fit.Ro = statsR.mean; - fit.dRo = 3600.0*statsR.sigma; - - fit.Do = statsD.mean; - fit.dDo = 3600.0*statsD.sigma; - - fit.chisq = 0.5*(statsR.chisq + statsD.chisq); - fit.Nfit = N; - - fit.uR = fit.duR = 0.0; - fit.uD = fit.duD = 0.0; - fit.p = fit.dp = 0.0; - - Nave ++; + catalog[i].average[j].flags |= ID_STAR_USE_AVE; + fit = fitAve; break; - case FIT_PM_ONLY: - FitPM (&fit, X, dX, Y, dY, T, N); - // fprintf (stderr, "fitted: %f - %f : %f %f : %f %f : %f\n", Tmin, Tmax, fit.Ro, fit.Do, fit.uR, fit.uD, fit.p); - // project Ro, Do back to RA,DEC - XY_to_RD (&fit.Ro, &fit.Do, fit.Ro, fit.Do, &coords); - // fprintf (stderr, "project: %f %f : %f %f : %f\n", fit.Ro, fit.Do, fit.uR, fit.uD, fit.p); - // continue; - - fit.p = fit.dp = 0.0; - - Npm ++; + catalog[i].average[j].flags |= ID_STAR_USE_PM; + fit = fitPM; break; - - case FIT_PAR_ONLY: - fprintf (stderr, "programming error at %s, %d", __FILE__, __LINE__); - exit (2); - - for (k = 0; k < N; k++) { - ParFactor (&pX[k], &pY[k], R[k], D[k], T[k]); - } - FitPar (&fit, X, dX, Y, dY, pX, pY, N); - - // project Ro, Do back to RA,DEC - XY_to_RD (&fit.Ro, &fit.Do, fit.Ro, fit.Do, &coords); - - fit.uR = fit.duR = 0.0; - fit.uD = fit.duD = 0.0; - - Npar ++; + case FIT_PM_AND_PAR: + catalog[i].average[j].flags |= ID_STAR_USE_PAR; + fit = fitPAR; break; - - case FIT_PM_AND_PAR: - fprintf (stderr, "programming error at %s, %d", __FILE__, __LINE__); - exit (2); - - for (k = 0; k < N; k++) { - ParFactor (&pX[k], &pY[k], R[k], D[k], T[k]); - } - FitPMandPar (&fit, X, dX, Y, dY, T, pX, pY, N); - XY_to_RD (&fit.Ro, &fit.Do, fit.Ro, fit.Do, &coords); - Npar ++; - break; - - default: - fprintf (stderr, "programming error at %s, %d", __FILE__, __LINE__); - exit (2); - } - - if (0 && (j < 100)) { - fprintf (stderr, "%f %f -> %f %f (%f,%f)\n", - catalog[i].average[j].R, - catalog[i].average[j].D, - fit.Ro, fit.Do, - 3600*(catalog[i].average[j].R - fit.Ro), - 3600*(catalog[i].average[j].D - fit.Do)); - } + } + + if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f)\n", + catalog[i].average[j].R, + catalog[i].average[j].D, + fit.Ro, fit.Do, + 3600*(catalog[i].average[j].R - fit.Ro), + 3600*(catalog[i].average[j].D - fit.Do)); //make sure that the fit succeeded - assert(finite(fit.Ro) && finite(fit.Do) && - finite(fit.dRo) && finite(fit.dDo) && - finite(fit.uR) && finite(fit.uD) && - finite(fit.duR) && finite(fit.duD) && - finite(fit.p) && finite(fit.dp)); + status = finite(fit.Ro); + status &= finite(fit.Do); + status &= finite(fit.dRo); + status &= finite(fit.dDo); + status &= finite(fit.uR); + status &= finite(fit.uD); + status &= finite(fit.duR); + status &= finite(fit.duD); + status &= finite(fit.p); + status &= finite(fit.dp); + if (!status) { + Nskip ++; + continue; + } // the measure fields must be updated before the average fields m = catalog[i].average[j].measureOffset; for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { - // XXX why was this here?? if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue; setMeanR (fit.Ro, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); setMeanD (fit.Do, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); @@ -270,11 +301,23 @@ catalog[i].average[j].dP = fit.dp; // parallax error in arcsec - catalog[i].average[j].Xp = (fit.Nfit > 1) ? 100.0*log10(fit.chisq) : NAN_S_SHORT; + // Xp is supposed to be the position scatter, not the chisq : fix this: + // catalog[i].average[j].Xp = (fit.Nfit > 1) ? 100.0*log10(fit.chisq) : NAN_S_SHORT; + catalog[i].average[j].ChiSqAve = fitAve.chisq; + catalog[i].average[j].ChiSqPM = fitPM.chisq; + catalog[i].average[j].ChiSqPar = fitPAR.chisq; + catalog[i].average[j].Xp = 0.0; + catalog[i].average[j].Tmean = (Tmean * 86400 * 365.26) + T2000; + catalog[i].average[j].Trange = (Trange * 86400 * 365.26); + catalog[i].average[j].Npos = fit.Nfit; } - if (VERBOSE) fprintf (stderr, "catalog %d : %d ave, %d pm, %d par : Nskip % d\n", i, Nave, Npm, Npar, Nskip); + NaveSum += Nave; + NpmSum += Npm; + NparSum += Npar; + NskipSum += Nskip; + if (VERBOSE) fprintf (stderr, "catalog %lld : %lld ave, %lld pm, %lld par : Nskip %lld\n", (long long) i, (long long) Nave, (long long) Npm, (long long) Npar, (long long) Nskip); } - if (VERBOSE) fprintf (stderr, "fitted %d objects (%d ave, %d pm, %d par)\n", Nave + Npm + Npar, Nave, Npm, Npar); + fprintf (stderr, "fitted %lld objects (%lld ave, %lld pm, %lld par), skipped %lld\n", (long long) (NaveSum + NpmSum + NparSum), (long long) NaveSum, (long long) NpmSum, (long long) NparSum, (long long) NskipSum); return (TRUE); } @@ -282,25 +325,25 @@ /* fitting proper-motion and parallax: -given a source at position r,d, at a time t, we need to calculate a vector (pr,pd) - -let x,y be the coordinate in the linearized frame with y parallel to DEC lines - -L,B are the ecliptic longitude and latitude of the object, -dL and dB are the offsets in the L and B directions - -dL = sin(t - topp) -dB = cos(t - topp)*sin(B) - -these need to be rotated to the R,D frame to yield pR,pD. Then, the equation of motion -for the source in the x,y frame is: - -x = Ro + uR * (t - to) + p * pR -y = Do + uD * (t - to) + p * pD - -the unknowns in these equations are Ro, uR, Do, uD, and p - -XXX think through the concepts for the pole a bit better. all objects near the pole -move the same way with the same phase. choose a projection center and define dL,dB relative -to that center point coordinate system? + given a source at position r,d, at a time t, we need to calculate a vector (pr,pd) + + let x,y be the coordinate in the linearized frame with y parallel to DEC lines + + L,B are the ecliptic longitude and latitude of the object, + dL and dB are the offsets in the L and B directions + + dL = sin(t - topp) + dB = cos(t - topp)*sin(B) + + these need to be rotated to the R,D frame to yield pR,pD. Then, the equation of motion + for the source in the x,y frame is: + + x = Ro + uR * (t - to) + p * pR + y = Do + uD * (t - to) + p * pD + + the unknowns in these equations are Ro, uR, Do, uD, and p + + XXX think through the concepts for the pole a bit better. all objects near the pole + move the same way with the same phase. choose a projection center and define dL,dB relative + to that center point coordinate system? */ Index: branches/simmosaic_branches/Ohana/src/relastro/src/UpdateSimple.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/UpdateSimple.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/UpdateSimple.c (revision 27839) @@ -4,5 +4,5 @@ /* we can measure new image parameters for each non-mosaic chip independently */ - int i, Nimage, Nstars; + off_t i, Nimage, Nstars; Image *image; StarData *raw, *ref; @@ -18,7 +18,9 @@ /* convert measure coordinates to raw entries */ raw = getImageRaw (catalog, Ncatalog, i, &Nstars, MODE_SIMPLE); + if (!raw) continue; /* convert average coordinates to ref entries */ ref = getImageRef (catalog, Ncatalog, i, &Nstars, MODE_SIMPLE); + if (!ref) continue; FitSimple (raw, ref, Nstars, &image[i].coords); Index: branches/simmosaic_branches/Ohana/src/relastro/src/args.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/args.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/args.c (revision 27839) @@ -50,4 +50,14 @@ if (FIT_TARGET == TARGET_NONE) usage(); + + UserCatalog = FALSE; + if ((N = get_argument (argc, argv, "-catalog"))) { + remove_argument (N, &argc, argv); + UserCatalogRA = atof(argv[N]); + remove_argument (N, &argc, argv); + UserCatalogDEC = atof(argv[N]); + remove_argument (N, &argc, argv); + UserCatalog = TRUE; + } /* specify portion of the sky : allow default of all sky? */ @@ -67,5 +77,7 @@ remove_argument (N, &argc, argv); } else { - usage (); + if (!UserCatalog) { + usage (); + } } @@ -111,7 +123,11 @@ } - VERBOSE = FALSE; + VERBOSE = VERBOSE2 = FALSE; if ((N = get_argument (argc, argv, "-v"))) { VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-vv"))) { + VERBOSE = VERBOSE2 = TRUE; remove_argument (N, &argc, argv); } @@ -242,4 +258,11 @@ } + NLOOP = 4; + if ((N = get_argument (argc, argv, "-nloop"))) { + remove_argument (N, &argc, argv); + NLOOP = atof (argv[N]); + remove_argument (N, &argc, argv); + } + if (argc != 1) usage (); return TRUE; @@ -248,4 +271,5 @@ void usage () { fprintf (stderr, "ERROR: USAGE: relastro -region RA RA DEC DEC\n"); + fprintf (stderr, " OR: relastro -catalog (ra) (dec)\n"); fprintf (stderr, " working options: \n"); fprintf (stderr, " -update-objects\n"); @@ -263,4 +287,5 @@ fprintf (stderr, " -statmode (mode)\n"); fprintf (stderr, " -reset"); + fprintf (stderr, " -nloop (N) : number of image-fit iterations"); fprintf (stderr, " -update : apply new fit to database\n"); fprintf (stderr, " -params\n"); Index: branches/simmosaic_branches/Ohana/src/relastro/src/bcatalog.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/bcatalog.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/bcatalog.c (revision 27839) @@ -3,9 +3,7 @@ int bcatalog (Catalog *subcatalog, Catalog *catalog) { - int i, j, k, offset, found; - int NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm, Nsecfilt; - float mag; - int mask; - PhotCode *code; + off_t i, j, offset; + off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm; + int Nsecfilt; // XXX in the future, use catalog[0].Nsecfilt only? allow catalogs to have variable Nsecfilt? @@ -106,6 +104,5 @@ if (VERBOSE) { - fprintf (stderr, "%d: using %d stars (%d measures) for catalog\n", i, - subcatalog[0].Naverage, subcatalog[0].Nmeasure); + fprintf (stderr, "%lld: using %lld stars (%lld measures) for catalog\n", (long long) i, (long long) subcatalog[0].Naverage, (long long) subcatalog[0].Nmeasure); } return (TRUE); Index: branches/simmosaic_branches/Ohana/src/relastro/src/fitpoly.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/fitpoly.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/fitpoly.c (revision 27839) @@ -103,5 +103,5 @@ /* convert the xsum,ysum,sum terms into vector,matrix and solve */ -void fit_eval (CoordFit *fit) { +int fit_eval (CoordFit *fit) { int i, j, ix, iy, jx, jy; @@ -146,6 +146,8 @@ // ix, iy, vector[i][0], ix, iy, vector[i][1]); } - - dgaussjordan (matrix, vector, fit[0].Nelems, 2); + + if (!dgaussjordan (matrix, vector, fit[0].Nelems, 2)) { + return (FALSE); + } for (i = 0; i < fit[0].Nelems; i++) { @@ -166,4 +168,5 @@ array_free (matrix, fit[0].Nelems); array_free (vector, fit[0].Nelems); + return (TRUE); } @@ -271,5 +274,5 @@ /* this should only apply to the polynomial, not the projection terms */ /* compare with psastro supporting code */ -CoordFit *fit_apply_coords (CoordFit *fit, Coords *coords) { +int fit_apply_coords (CoordFit *fit, Coords *coords) { double Xo, Yo, R1, R2; @@ -281,5 +284,8 @@ // L = pc1_1*cd1*(x - cp1) + pc1_2*cd2*(y - cp2) + ... - CoordsGetCenter (fit, 0.001, &Xo, &Yo); + if (!CoordsGetCenter (fit, 0.001, &Xo, &Yo)) { + fprintf (stderr, "failed to modify model\n"); + return (FALSE); + } coords[0].crpix1 = Xo; coords[0].crpix2 = Yo; @@ -329,6 +335,5 @@ /* keep the order and type from initial values */ - // XXX if desired in the future, return modfit (and free above) - return (NULL); -} - + return (TRUE); +} + Index: branches/simmosaic_branches/Ohana/src/relastro/src/initialize.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/initialize.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/initialize.c (revision 27839) @@ -5,4 +5,6 @@ int NPHOTCODES; char *codename, *ptr, *list; + + ptr = NULL; ConfigInit (&argc, argv); Index: branches/simmosaic_branches/Ohana/src/relastro/src/load_images.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/load_images.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/load_images.c (revision 27839) @@ -1,30 +1,49 @@ # include "relastro.h" + +# define MARKTIME(MSG,...) { \ + float dtime; \ + gettimeofday (&stop, (void *) NULL); \ + dtime = DTIME (stop, start); \ + fprintf (stderr, MSG, __VA_ARGS__); } SkyList *load_images (FITS_DB *db, SkyRegion *region) { Image *image, *subset; - int Nimage, Nsubset; - int *LineNumber; + off_t Nimage, Nsubset; + off_t *LineNumber; + struct timeval start, stop; SkyTable *sky = NULL; SkyList *skylist = NULL; + gettimeofday (&start, (void *) NULL); + // load the current sky table (layout of all SkyRegions) - sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, SKY_DEPTH, VERBOSE); + sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); SkyTableSetFilenames (sky, CATDIR, "cpt"); // determine the populated SkyRegions overlapping the requested area - skylist = SkyListByPatch (sky, -1, region); + if (UserCatalog) { + skylist = SkyRegionByPoint (sky, -1, UserCatalogRA, UserCatalogDEC); + } else { + skylist = SkyListByPatch (sky, -1, region); + } + MARKTIME(" setup sky: %f sec\n", dtime); // convert database table to internal structure image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped); + MARKTIME(" convert image table: %f sec\n", dtime); // select the images which overlap the selected sky regions subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset); + MARKTIME(" select images: %f sec\n", dtime); gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, LineNumber, Nsubset); initImages (subset, Nsubset); + MARKTIME(" init images: %f sec\n", dtime); + initMosaics (subset, Nsubset); + MARKTIME(" init mosaics: %f sec\n", dtime); /* unlock, if we can (else, unlocked below) */ Index: branches/simmosaic_branches/Ohana/src/relastro/src/mkpolyterm.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/mkpolyterm.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/mkpolyterm.c (revision 27839) @@ -73,4 +73,5 @@ if (dPos > dPosRef) { fprintf (stderr, "*** warning : non-convergence in model conversion (mkpolyterm.c;73) *** \n"); + return FALSE; } array_free (alpha, 2); Index: branches/simmosaic_branches/Ohana/src/relastro/src/plotstuff.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/plotstuff.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/plotstuff.c (revision 27839) @@ -125,5 +125,5 @@ if (PLOTDELAY > 500000) { fprintf (stdout, "press return\n"); - fscanf (stdin, "%*c"); + if (fscanf (stdin, "%*c") != 0) fprintf (stderr, "\n"); } else { usleep (PLOTDELAY); Index: branches/simmosaic_branches/Ohana/src/relastro/src/relastro.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/relastro.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/relastro.c (revision 27839) @@ -1,11 +1,19 @@ # include "relastro.h" + +# define MARKTIME(MSG,...) { \ + float dtime; \ + gettimeofday (&stop, (void *) NULL); \ + dtime = DTIME (stop, start); \ + fprintf (stderr, MSG, __VA_ARGS__); } int main (int argc, char **argv) { - int status, Ncatalog; + int i, status, Ncatalog; Catalog *catalog; FITS_DB db; + struct timeval start, stop; + SkyList *skylist = NULL; - SkyList *skylist = NULL; + gettimeofday (&start, (void *) NULL); /* get configuration info, args */ @@ -26,14 +34,22 @@ if (db.dbstate == LCK_EMPTY) Shutdown ("ERROR: No images in catalog %s (1)", db.filename); if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename); + MARKTIME("load image data: %f sec\n", dtime); /* load regions and images based on specified sky patch (default depth) */ skylist = load_images (&db, &UserPatch); + MARKTIME("load images: %f sec\n", dtime); + + initCoords(); /* load catalog data from region files : subselect high-quality measurements */ catalog = load_catalogs (skylist, &Ncatalog, TRUE); + MARKTIME("load catalog data: %f sec\n", dtime); /* match measurements with images */ - initImageBins (catalog, Ncatalog); - findImages (catalog, Ncatalog); + initImageBins (catalog, Ncatalog, TRUE); + MARKTIME("make image bins: %f sec\n", dtime); + + findImages (catalog, Ncatalog, TRUE); + MARKTIME("set up image indexes: %f sec\n", dtime); if (PLOTSTUFF) { @@ -45,13 +61,23 @@ switch (FIT_TARGET) { case TARGET_SIMPLE: - UpdateSimple (catalog, Ncatalog); + for (i = 0; i < NLOOP; i++) { + UpdateObjects (catalog, Ncatalog); + UpdateSimple (catalog, Ncatalog); + } break; case TARGET_CHIPS: - UpdateChips (catalog, Ncatalog); + for (i = 0; i < NLOOP; i++) { + UpdateObjects (catalog, Ncatalog); + UpdateChips (catalog, Ncatalog); + MARKTIME("update chips: %f sec\n", dtime); + } break; case TARGET_MOSAICS: - UpdateMosaic (catalog, Ncatalog); + for (i = 0; i < NLOOP; i++) { + UpdateObjects (catalog, Ncatalog); + UpdateMosaic (catalog, Ncatalog); + } break; @@ -69,4 +95,7 @@ UpdateObjectOffsets (skylist); + // iterate over catalogs to make detection coordinates consistant + FixProblemImages (skylist); + // save the updated image parameters dvo_image_update (&db, VERBOSE); Index: branches/simmosaic_branches/Ohana/src/relastro/src/relastro_objects.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/relastro_objects.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/relastro_objects.c (revision 27839) @@ -11,9 +11,13 @@ // load the current sky table (layout of all SkyRegions) - sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, SKY_DEPTH, VERBOSE); + sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); SkyTableSetFilenames (sky, CATDIR, "cpt"); // determine the populated SkyRegions overlapping the requested area (default depth) - skylist = SkyListByPatch (sky, -1, &UserPatch); + if (UserCatalog) { + skylist = SkyRegionByPoint (sky, -1, UserCatalogRA, UserCatalogDEC); + } else { + skylist = SkyListByPatch (sky, -1, &UserPatch); + } // load data from each region file, only use bright stars Index: branches/simmosaic_branches/Ohana/src/relastro/src/select_images.c =================================================================== --- branches/simmosaic_branches/Ohana/src/relastro/src/select_images.c (revision 24860) +++ branches/simmosaic_branches/Ohana/src/relastro/src/select_images.c (revision 27839) @@ -13,14 +13,20 @@ } SkyRegionCoords; -Image *select_images (SkyList *skylist, Image *timage, int Ntimage, int **LineNumber, int *Nimage) { +void dsortindex (double *X, off_t *Y, int N); +off_t getRegionStartByRA (double R, double *Rref, off_t Nregions); + +Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage) { Image *image; - int i, j, k, m, found, nimage, NIMAGE; - int InRange; + off_t i, j, k, m, nStart, iSky, nimage, NIMAGE; + off_t *line_number; + int InRange, found; double Ri[5], Di[5], Xi[5], Yi[5], dx, dy; - int *line_number; Coords tcoords; SkyRegionCoords *skycoords; + double *RmaxSky; + off_t *index; + if (skylist[0].Nregions < 1) { *Nimage = 0; @@ -38,4 +44,7 @@ ALLOCATE (skycoords, SkyRegionCoords, skylist[0].Nregions); + + ALLOCATE (RmaxSky, double, skylist[0].Nregions); + ALLOCATE (index, off_t, skylist[0].Nregions); /* compare with each region file */ @@ -56,4 +65,7 @@ skycoords[i].Yc[4] = skycoords[i].Yc[0]; + RmaxSky[i] = skylist[0].regions[i][0].Rmax; + index[i] = i; + dx = 0.02*(skycoords[i].Xc[2] - skycoords[i].Xc[0]); dy = 0.02*(skycoords[i].Yc[2] - skycoords[i].Yc[0]); @@ -65,4 +77,6 @@ } + dsortindex (RmaxSky, index, skylist[0].Nregions); + if (VERBOSE) fprintf (stderr, "finding images\n"); BuildChipMatch (timage, Ntimage); @@ -71,5 +85,5 @@ NIMAGE = 100; ALLOCATE (image, Image, NIMAGE); - ALLOCATE (line_number, int, NIMAGE); + ALLOCATE (line_number, off_t, NIMAGE); // go through the complete list of images, selecting ones which overlap any region @@ -100,5 +114,8 @@ } - if (!FindMosaicForImage (timage, Ntimage, i)) continue; + if (!FindMosaicForImage (timage, Ntimage, i)) { + fprintf (stderr, "cannot find mosaic for %lld\n", (long long) i); + continue; + } /* define image corners */ @@ -110,11 +127,18 @@ found = FALSE; - /* transform to ra,dec */ + /* transform corners to ra,dec */ + double RminImage = 360.0; for (j = 0; j < 5; j++) { XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords); - } + RminImage = MIN(RminImage, Ri[j]); + } + + // RA(nStart) is guaranteed to be < RminImage: + nStart = getRegionStartByRA (RminImage, RmaxSky, skylist[0].Nregions); /* compare with each region file */ - for (m = 0; (m < skylist[0].Nregions) && !found; m++) { + for (iSky = nStart; (iSky < skylist[0].Nregions) && !found; iSky++) { + + m = index[iSky]; /* we make positional comparisons in the projection of catalog */ @@ -122,5 +146,5 @@ tcoords.crval2 = skycoords[m].Dc; - /* transform to ra,dec */ + /* transform corner coords to X,Y in this catalog system */ InRange = TRUE; for (j = 0; (j < 5) && InRange; j++) { @@ -151,5 +175,5 @@ } if (RESET) { - // XXX do we need this in relastro? + // XXX do we need / want to do this in relastro? assignMcal (&image[nimage], (double *) NULL, -1); image[nimage].dMcal = NAN; @@ -161,13 +185,13 @@ NIMAGE += 100; REALLOCATE (image, Image, NIMAGE); - REALLOCATE (line_number, int, NIMAGE); + REALLOCATE (line_number, off_t, NIMAGE); } } } - if (VERBOSE) fprintf (stderr, "found %d images\n", nimage); + if (VERBOSE) fprintf (stderr, "found %lld images\n", (long long) nimage); REALLOCATE (image, Image, MAX (nimage, 1)); - REALLOCATE (line_number, int, MAX (nimage, 1)); + REALLOCATE (line_number, off_t, MAX (nimage, 1)); free (skycoords); @@ -253,2 +277,35 @@ } +void dsortindex (double *X, off_t *Y, int N) { + +# define SWAPFUNC(A,B){ double tmpf; off_t tmpi; \ + tmpf = X[A]; X[A] = X[B]; X[B] = tmpf; \ + tmpi = Y[A]; Y[A] = Y[B]; Y[B] = tmpi; \ +} +# define COMPARE(A,B)(X[A] < X[B]) + + OHANA_SORT (N, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + +} + +off_t getRegionStartByRA (double R, double *Rref, off_t Nregions) { + + // use bisection to find the overlapping mosaic + + off_t Nlo, Nhi, N; + + // find the last mosaic before start + Nlo = 0; Nhi = Nregions; + while (Nhi - Nlo > 10) { + N = 0.5*(Nlo + Nhi); + if (Rref[N] < R) { + Nlo = MAX(N, 0); + } else { + Nhi = MIN(N, Nregions); + } + } + return (Nlo); +}