Index: /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/Makefile =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/Makefile (revision 36744) +++ /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/Makefile (revision 36745) @@ -31,4 +31,5 @@ $(SRC)/args.$(ARCH).o \ $(SRC)/copy_images.$(ARCH).o \ +$(SRC)/select_images.$(ARCH).o \ $(SRC)/Shutdown.$(ARCH).o \ $(SRC)/join_stars.$(ARCH).o \ @@ -41,5 +42,4 @@ $(SRC)/ConfigInit.$(ARCH).o \ $(SRC)/args.$(ARCH).o \ -$(SRC)/copy_images.$(ARCH).o \ $(SRC)/Shutdown.$(ARCH).o \ $(SRC)/join_stars.$(ARCH).o \ Index: /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/include/photdbc.h =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/include/photdbc.h (revision 36744) +++ /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/include/photdbc.h (revision 36745) @@ -20,4 +20,20 @@ int Nmeas; } StatType; + +typedef struct { + double Xc[5]; + double Yc[5]; + double Rc; + double Dc; +} SkyRegionCoords; + +# define MARKTIME(MSG,...) { \ + gettimeofday (&stopTimer, (void *) NULL); \ + float dtime = DTIME (stopTimer, startTimer); \ + fprintf (stderr, MSG, __VA_ARGS__); } + +# define INITTIME \ + struct timeval startTimer, stopTimer; \ + gettimeofday (&startTimer, (void *) NULL); int PARALLEL; @@ -80,7 +96,7 @@ PhotCodeData photcodes; -char *PHOTCODE_DROP_LIST, *PHOTCODE_SKIP_LIST; -int NphotcodesDrop, NphotcodesSkip; -PhotCode **photcodesDrop, **photcodesSkip; +char *PHOTCODE_DROP_LIST, *PHOTCODE_KEEP_LIST; +int NphotcodesDrop, NphotcodesKeep; +PhotCode **photcodesDrop, **photcodesKeep; # define FLAG_AREA 0X0001 @@ -131,6 +147,11 @@ void SetProtect (int mode); int SetSignals (void); -int copy_images (char *outdir); +int copy_images (char *outdir, SkyList *skylist); void usage(); + +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, int UseFullOverlap); int photdbc_catalogs (char *outroot, SkyList *skylist, int hostID); Index: /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/src/args.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/src/args.c (revision 36744) +++ /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/src/args.c (revision 36745) @@ -85,4 +85,5 @@ } + // measurements with these photcodes are not copied to the output PHOTCODE_DROP_LIST = NULL; if ((N = get_argument (argc, argv, "-photcode-drop"))) { @@ -92,8 +93,11 @@ } - PHOTCODE_SKIP_LIST = NULL; - if ((N = get_argument (argc, argv, "-photcode-skip"))) { - remove_argument (N, &argc, argv); - PHOTCODE_SKIP_LIST = strcreate(argv[N]); + // measurements with these photcodes are kept ***regardless of quality*** + // -photcode-keep J will keep all J-band measurements of all kinds + // -photcode-keep GPC1.02.g will keep all g-band measurements from chip XY02 + PHOTCODE_KEEP_LIST = NULL; + if ((N = get_argument (argc, argv, "-photcode-keep"))) { + remove_argument (N, &argc, argv); + PHOTCODE_KEEP_LIST = strcreate(argv[N]); remove_argument (N, &argc, argv); } @@ -254,8 +258,8 @@ } - PHOTCODE_SKIP_LIST = NULL; - if ((N = get_argument (argc, argv, "-photcode-skip"))) { - remove_argument (N, &argc, argv); - PHOTCODE_SKIP_LIST = strcreate(argv[N]); + PHOTCODE_KEEP_LIST = NULL; + if ((N = get_argument (argc, argv, "-photcode-keep"))) { + remove_argument (N, &argc, argv); + PHOTCODE_KEEP_LIST = strcreate(argv[N]); remove_argument (N, &argc, argv); } @@ -290,5 +294,5 @@ fprintf (stderr, " -photcode-drop : remove these photcodes from the output (REF or DEP only)\n"); - fprintf (stderr, " -photcode-skip : ignore these photcodes when assessing the validity (keep unless object is dropped)\n"); + fprintf (stderr, " -photcode-keep : ignore these photcodes when assessing the validity (keep unless object is dropped)\n"); fprintf (stderr, " -instmag (min) (max) : range of valid instrumental magnitudes (or measurements are dropped)\n"); Index: /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/src/copy_images.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/src/copy_images.c (revision 36744) +++ /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/src/copy_images.c (revision 36745) @@ -1,13 +1,14 @@ # include "photdbc.h" -int copy_images (char *outdir) { +int copy_images (char *outdir, SkyList *skylist) { int status; - off_t Nimage; + off_t Nimage, Nsubset; + off_t *LineNumber; char *ImageOut; unsigned int imageID; FITS_DB in; FITS_DB out; - Image *image; + Image *image, *subset; struct stat filestat; char *path; @@ -50,5 +51,8 @@ exit (2); } - dvo_image_addrows (&out, image, Nimage); + + subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset, FALSE); + + dvo_image_addrows (&out, subset, Nsubset); // note that imageID is unsigned int Index: /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/src/initialize.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/src/initialize.c (revision 36744) +++ /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/src/initialize.c (revision 36745) @@ -13,5 +13,5 @@ photcodesDrop = ParsePhotcodeList (PHOTCODE_DROP_LIST, &NphotcodesDrop, FALSE); - photcodesSkip = ParsePhotcodeList (PHOTCODE_SKIP_LIST, &NphotcodesSkip, FALSE); + photcodesKeep = ParsePhotcodeList (PHOTCODE_KEEP_LIST, &NphotcodesKeep, FALSE); if (SHOW_PARAMS) { @@ -54,5 +54,5 @@ photcodesDrop = ParsePhotcodeList (PHOTCODE_DROP_LIST, &NphotcodesDrop, FALSE); - photcodesSkip = ParsePhotcodeList (PHOTCODE_SKIP_LIST, &NphotcodesSkip, FALSE); + photcodesKeep = ParsePhotcodeList (PHOTCODE_KEEP_LIST, &NphotcodesKeep, FALSE); if (SHOW_PARAMS) { Index: /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/src/make_subcatalog.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/src/make_subcatalog.c (revision 36744) +++ /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/src/make_subcatalog.c (revision 36745) @@ -9,5 +9,5 @@ off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm, Nsecfilt; double mag, minMag, minSigma; - int keep, *secSkip; + int keep, *secKeep; PhotCode *photcode; @@ -16,11 +16,11 @@ // set up a list of SEC entries to ignore when evaluating a source - ALLOCATE (secSkip, int, Nsecfilt); + ALLOCATE (secKeep, int, Nsecfilt); for (i = 0; i < Nsecfilt; i++) { - secSkip[i] = FALSE; + secKeep[i] = FALSE; photcode = GetPhotcodebyNsec(i); - for (k = 0; k < NphotcodesSkip; k++) { - if (photcodesSkip[k][0].code != photcode[0].code) continue; - secSkip[i] = TRUE; + for (k = 0; k < NphotcodesKeep; k++) { + if (photcodesKeep[k][0].code != photcode[0].code) continue; + secKeep[i] = TRUE; } } @@ -67,5 +67,5 @@ keep = FALSE; for (j = 0; !keep && (j < Nsecfilt); j++) { - if (secSkip[j]) continue; + if (secKeep[j]) continue; if (catalog[0].secfilt[Nsecfilt*i+j].Ncode >= NCODE_MIN) { keep = TRUE; @@ -101,7 +101,7 @@ if (NphotcodesDrop > 0) { found = FALSE; - for (k = 0; (k < NphotcodesSkip) && !found; k++) { - if (photcodesSkip[k][0].code == catalog[0].measure[offset].photcode) found = TRUE; - if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE; + for (k = 0; (k < NphotcodesDrop) && !found; k++) { + if (photcodesDrop[k][0].code == catalog[0].measure[offset].photcode) found = TRUE; + if (photcodesDrop[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE; } if (found) continue; @@ -109,9 +109,9 @@ // ignore certain photcodes to assess the measurements - if (NphotcodesSkip > 0) { + if (NphotcodesKeep > 0) { found = FALSE; - for (k = 0; (k < NphotcodesSkip) && !found; k++) { - if (photcodesSkip[k][0].code == catalog[0].measure[offset].photcode) found = TRUE; - if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE; + for (k = 0; (k < NphotcodesKeep) && !found; k++) { + if (photcodesKeep[k][0].code == catalog[0].measure[offset].photcode) found = TRUE; + if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE; } if (found) goto keep; Index: /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/src/photdbc.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/src/photdbc.c (revision 36744) +++ /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/src/photdbc.c (revision 36745) @@ -9,11 +9,11 @@ initialize (argc, argv); - // load and copy the image table - copy_images (argv[1]); - // the output catalog needs to inherit the SKY_DEPTH of the input catalog sky = SkyTableLoadOptimal (CATDIR, NULL, NULL, TRUE, 0, VERBOSE); SkyTableSetFilenames (sky, CATDIR, "cpt"); skylist = SkyListByPatch (sky, -1, ®ION); + + // load and copy the image table + copy_images (argv[1], skylist); // hostID is 0 for master program Index: /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/src/photdbc_catalogs.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/src/photdbc_catalogs.c (revision 36744) +++ /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/src/photdbc_catalogs.c (revision 36745) @@ -120,5 +120,5 @@ if (ExcludeByMaxMinMag) { snprintf (tmpline, DVO_MAX_PATH, "%s -maxminmag %f", command, MAX_MIN_MAG); strcpy (command, tmpline); } if (PHOTCODE_DROP_LIST) { snprintf (tmpline, DVO_MAX_PATH, "%s -photcode-drop %s", command, PHOTCODE_DROP_LIST); strcpy (command, tmpline); } - if (PHOTCODE_SKIP_LIST) { snprintf (tmpline, DVO_MAX_PATH, "%s -photcode-skip %s", command, PHOTCODE_SKIP_LIST); strcpy (command, tmpline); } + if (PHOTCODE_KEEP_LIST) { snprintf (tmpline, DVO_MAX_PATH, "%s -photcode-keep %s", command, PHOTCODE_KEEP_LIST); strcpy (command, tmpline); } if (CATFORMAT) { snprintf (tmpline, DVO_MAX_PATH, "%s -set-format %s", command, CATFORMAT); strcpy (command, tmpline); } if (CATMODE) { snprintf (tmpline, DVO_MAX_PATH, "%s -set-mode %s", command, CATMODE); strcpy (command, tmpline); } Index: /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/src/select_images.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/src/select_images.c (revision 36745) +++ /branches/eam_branches/ipp-20140423/Ohana/src/photdbc/src/select_images.c (revision 36745) @@ -0,0 +1,234 @@ +# include "photdbc.h" + +/* this function returns a list of all images which overlap the given SkyList (set of + SkyRegions). All images in the image catalog are tested once, so there is no check that an + image already has been included. LineNum stores the locations in the Image database of the + list of images */ + +Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage, int UseFullOverlap) { + + Image *image; + off_t i, j, k, nimage, NIMAGE, D_NIMAGE; + off_t *line_number; + int found; + double Ri[5], Di[5], Xi[5], Yi[5], dx, dy; + Coords tcoords; + SkyRegionCoords *skycoords; + + double RmaxSkyRegion, RminSkyRegion, RmidSkyRegion, DminSkyRegion, DmaxSkyRegion; + + double *RmaxSky; + off_t *index; + + if (skylist[0].Nregions < 1) { + *Nimage = 0; + *LineNumber = NULL; + if (VERBOSE) fprintf (stderr, "no matching sky regions\n"); + return NULL; + } + + // if no region is selected, we are choosing the whole sky; I should sort-circuit + // this function in that case + + INITTIME; + + // the comparison is made in the catalog local projection. below we set crval1,2 + tcoords.crpix1 = tcoords.crpix2 = 0.0; + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + tcoords.pc1_1 = tcoords.pc2_2 = 1.0; + tcoords.pc1_2 = tcoords.pc2_1 = 0.0; + strcpy (tcoords.ctype, "DEC--TAN"); + + ALLOCATE (skycoords, SkyRegionCoords, skylist[0].Nregions); + + ALLOCATE (RmaxSky, double, skylist[0].Nregions); + ALLOCATE (index, off_t, skylist[0].Nregions); + + RminSkyRegion = +360.0; + RmaxSkyRegion = -360.0; + DminSkyRegion = +90.0; + DmaxSkyRegion = -90.0; + + D_NIMAGE = 6000; + + /* compare with each region file */ + for (i = 0; i < skylist[0].Nregions; i++) { + + /* we make positional comparisons in the projection of catalog */ + skycoords[i].Rc = 0.5*(skylist[0].regions[i][0].Rmax + skylist[0].regions[i][0].Rmin); + skycoords[i].Dc = 0.5*(skylist[0].regions[i][0].Dmax + skylist[0].regions[i][0].Dmin); + tcoords.crval1 = skycoords[i].Rc; + tcoords.crval2 = skycoords[i].Dc; + + /* define catalog corners */ + RD_to_XY (&skycoords[i].Xc[0], &skycoords[i].Yc[0], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmin, &tcoords); + RD_to_XY (&skycoords[i].Xc[1], &skycoords[i].Yc[1], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmin, &tcoords); + RD_to_XY (&skycoords[i].Xc[2], &skycoords[i].Yc[2], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmax, &tcoords); + RD_to_XY (&skycoords[i].Xc[3], &skycoords[i].Yc[3], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmax, &tcoords); + skycoords[i].Xc[4] = skycoords[i].Xc[0]; + 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]); + skycoords[i].Xc[0] -= dx; skycoords[i].Yc[0] -= dy; + skycoords[i].Xc[1] += dx; skycoords[i].Yc[1] -= dy; + skycoords[i].Xc[2] += dx; skycoords[i].Yc[2] += dy; + skycoords[i].Xc[3] -= dx; skycoords[i].Yc[3] += dy; + skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy; + + RminSkyRegion = MIN(RminSkyRegion, skylist[0].regions[i][0].Rmin); + RmaxSkyRegion = MAX(RmaxSkyRegion, skylist[0].regions[i][0].Rmax); + DminSkyRegion = MIN(DminSkyRegion, skylist[0].regions[i][0].Dmin); + DmaxSkyRegion = MAX(DmaxSkyRegion, skylist[0].regions[i][0].Dmax); + } + RmidSkyRegion = 0.5*(RminSkyRegion + RmaxSkyRegion); + MARKTIME("create sky region coords: %f sec\n", dtime); + + dsortindex (RmaxSky, index, skylist[0].Nregions); + MARKTIME("sort sky coords: %f sec\n", dtime); + + if (VERBOSE) fprintf (stderr, "finding images\n"); + BuildChipMatch (timage, Ntimage); + MARKTIME("build chip match: %f sec\n", dtime); + + nimage = 0; + NIMAGE = D_NIMAGE; + ALLOCATE (image, Image, NIMAGE); + ALLOCATE (line_number, off_t, NIMAGE); + + // go through the complete list of images, selecting ones which overlap any region + for (i = 0; i < Ntimage; i++) { + + /* exclude images by photcode, or equiv photcode, if specified */ + if (NphotcodesDrop > 0) { + if (timage[i].photcode == 0) goto keep_image; // PHU image can be immediately passed in here + found = FALSE; // here 'found' refers to a match to a drop photcode + for (k = 0; (k < NphotcodesDrop) && !found; k++) { + if (photcodesDrop[k][0].code == timage[i].photcode) found = TRUE; + if (photcodesDrop[k][0].code == GetPhotcodeEquivCodebyCode(timage[i].photcode)) found = TRUE; + } + if (found) continue; + } +keep_image: + + // this adds 1.3 sec for 3M images + if (!FindMosaicForImage (timage, Ntimage, i)) { + fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i); + continue; + } + + /* define image corners - note the DIS images (mosaic phu) are special */ + if (!strcmp(&timage[i].coords.ctype[4], "-DIS")) { + Xi[0] = -0.5*timage[i].NX; Yi[0] = -0.5*timage[i].NY; + Xi[1] = +0.5*timage[i].NX; Yi[1] = -0.5*timage[i].NY; + Xi[2] = +0.5*timage[i].NX; Yi[2] = +0.5*timage[i].NY; + Xi[3] = -0.5*timage[i].NX; Yi[3] = +0.5*timage[i].NY; + Xi[4] = -0.5*timage[i].NX; Yi[4] = -0.5*timage[i].NY; + } else { + Xi[0] = 0; Yi[0] = 0; + Xi[1] = timage[i].NX; Yi[1] = 0; + Xi[2] = timage[i].NX; Yi[2] = timage[i].NY; + Xi[3] = 0; Yi[3] = timage[i].NY; + Xi[4] = 0; Yi[4] = 0; + } + found = FALSE; + + /* transform corners to ra,dec -- costs ~3sec for 3M images (pikake) */ + double RminImage = RmidSkyRegion + 180.0; + double RmaxImage = RmidSkyRegion - 180.0; + double DminImage = +90.0; + double DmaxImage = -90.0; + for (j = 0; j < 5; j++) { + XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords); + Ri[j] = ohana_normalize_angle_to_midpoint (Ri[j], RmidSkyRegion); + + RminImage = MIN(RminImage, Ri[j]); + RmaxImage = MAX(RmaxImage, Ri[j]); + DminImage = MIN(DminImage, Di[j]); + DmaxImage = MAX(DmaxImage, Di[j]); + } + if (RmaxImage - RminImage > 180.0) { + double tmp = RminImage; + RmaxImage = RminImage; + RminImage = tmp - 360.0; + } + + // check that this image is even in range of the searched region + if (DminImage > DmaxSkyRegion) continue; + if (DmaxImage < DminSkyRegion) continue; + + // the sky region RA is defined to be 0 - 360.0 + if (RminImage > RmaxSkyRegion) continue; + if (RmaxImage < RminSkyRegion) continue; + + // require that the full image be inside the region of interest (only for non PHU + // images) XXX : if we calibrate the mosaic, require that the full mosaic be inside + // the region and all of its chips.. + if (UseFullOverlap && strcmp(&timage[i].coords.ctype[4], "-DIS")) { + if (RmaxImage > REGION.Rmax) continue; + if (RminImage < REGION.Rmin) continue; + if (DmaxImage > REGION.Dmax) continue; + if (DminImage < REGION.Dmin) continue; + } + + image[nimage] = timage[i]; + line_number[nimage] = i; + nimage ++; + if (nimage == NIMAGE) { + NIMAGE += D_NIMAGE; + D_NIMAGE = MAX (100000, D_NIMAGE * 1.5); + REALLOCATE (image, Image, NIMAGE); + REALLOCATE (line_number, off_t, NIMAGE); + } + } + MARKTIME("finish image selection: %f sec\n", dtime); + + if (VERBOSE) fprintf (stderr, "found "OFF_T_FMT" images\n", nimage); + + REALLOCATE (image, Image, MAX (nimage, 1)); + REALLOCATE (line_number, off_t, MAX (nimage, 1)); + free (skycoords); + free (RmaxSky); + free (index); + + *Nimage = nimage; + *LineNumber = line_number; + return (image); +} + +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); +}