IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 11, 2015, 2:14:39 PM (12 years ago)
Author:
eugene
Message:

merge changes from eam branch ipp-20140904

Location:
trunk/Ohana
Files:
36 edited
9 copied

Legend:

Unmodified
Added
Removed

  • trunk/Ohana

  • trunk/Ohana/src/relastro/Makefile

    r37261 r37807  
    6060$(SRC)/CoordOps.$(ARCH).o            \
    6161$(SRC)/extra.$(ARCH).o            \
    62 $(SRC)/FixProblemImages.$(ARCH).o    \
    6362$(SRC)/StarMaps.$(ARCH).o    \
    6463$(SRC)/high_speed_catalogs.$(ARCH).o  \
    … …  
    7574$(SRC)/relastro_parallel_regions.$(ARCH).o \
    7675$(SRC)/MeanPosIO.$(ARCH).o \
     76$(SRC)/MeasPosIO.$(ARCH).o \
    7777$(SRC)/share_mean_pos.$(ARCH).o \
     78$(SRC)/share_meas_pos.$(ARCH).o \
    7879$(SRC)/share_images_pos.$(ARCH).o \
     80$(SRC)/share_icrf_obj.$(ARCH).o \
     81$(SRC)/ICRF.$(ARCH).o \
     82$(SRC)/ICRFobjIO.$(ARCH).o \
     83$(SRC)/FrameCorrection.$(ARCH).o \
     84$(SRC)/FrameCorrectionIO.$(ARCH).o \
     85$(SRC)/FrameCorrectionUtils.$(ARCH).o \
    7986$(SRC)/ImagePosIO.$(ARCH).o \
    8087$(SRC)/ImageTable.$(ARCH).o \
    … …  
    8289$(SRC)/indexCatalogs.$(ARCH).o \
    8390$(SRC)/syncfile.$(ARCH).o \
     91$(SRC)/client_logger.$(ARCH).o \
    8492$(SRC)/relastroVisual.$(ARCH).o
    8593
    … …  
    117125$(SRC)/relastro_objects.$(ARCH).o    \
    118126$(SRC)/save_catalogs.$(ARCH).o       \
    119 $(SRC)/CoordOps.$(ARCH).o            \
    120127$(SRC)/extra.$(ARCH).o            \
    121128$(SRC)/high_speed_catalogs.$(ARCH).o  \
    … …  
    128135$(SRC)/plot_scatter.$(ARCH).o        \
    129136$(SRC)/plotstuff.$(ARCH).o           \
     137$(SRC)/ICRF.$(ARCH).o \
    130138$(SRC)/relastroVisual.$(ARCH).o \
    131139$(SRC)/syncfile.$(ARCH).o \
    132140$(SRC)/BrightCatalog.$(ARCH).o
    133 
    134 #$(SRC)/FitChip.$(ARCH).o             \
    135 #$(SRC)/FitMosaic.$(ARCH).o           \
    136 #$(SRC)/UpdateChips.$(ARCH).o         \
    137 #$(SRC)/UpdateMosaic.$(ARCH).o        \
    138 #$(SRC)/UpdateSimple.$(ARCH).o        \
    139 #$(SRC)/write_coords.$(ARCH).o        \
    140 #$(SRC)/FixProblemImages.$(ARCH).o    \
    141 #$(SRC)/StarMaps.$(ARCH).o    \
    142141
    143142$(RELASTRO_CLIENT): $(INC)/relastro.h

  • trunk/Ohana/src/relastro/doc/pv2.txt

    r37038 r37807  
     1
     22014.10.28
     3
     4Going over the sequence of steps in relastro.  I'm concerned that the
     5iteration is not modifying the measure.R,D values, preventing any
     6actual iterative improvements.
     7
     8* UpdateObjects : takes measure.R,D to calculate average.R,D and
     9  motion terms
     10
     11* ImageOps :
     12  * getImageRaw : takes measure.X,Y, applies image.coords to get (L,M),
     13                  (P,Q), etc, calculates
     14  * getImageRef : takes average.R,D, applies image.coords to get
     15                  (P,Q), (L,M), etc
     16  * FitChip : fit raw.X,Y to ref.L,M, plus residual map (dX,dY)
     17 
     18  *** modified raw.L,M values are then dropped and not used to update measure.R.D
    119
    2202014.07.03

  • trunk/Ohana/src/relastro/include/relastro.h

    r37261 r37807  
    3434
    3535typedef struct {
    36   float R;
    37   float D;
     36  double R;
     37  double D;
    3838  unsigned int objID;
    3939  unsigned int catID;
    4040} MeanPos;
     41
     42typedef struct {
     43  double R;
     44  double D;
     45  unsigned int objID;
     46  unsigned int catID;
     47  unsigned int imageID;
     48} MeasPos;
    4149
    4250typedef struct {
    … …  
    4755  int nFitAstrom;
    4856  int flags;
     57  float refColorBlue;
     58  float refColorRed;
    4959} ImagePos;
    5060
    … …  
    142152} FitStats;
    143153
     154typedef struct {
     155  double scale;
     156  double **Roff;
     157  double **Doff;
     158  double  *dR;
     159  int *Nra;
     160  int Ndec;
     161} FrameCorrectionType;
     162
     163typedef struct {
     164  // AstromOffsetMap *map;
     165  Coords *coords; // carries a pointer to the AstromOffsetMap
     166  FrameCorrectionType *frame;
     167} FrameCorrectionSet;
     168
     169
     170typedef struct {
     171  double *Rave;
     172  double *Dave;
     173  double *dRoff;
     174  double *dDoff;
     175  int Nicrfobj;
     176} ICRFobj;
     177
    144178/* global variables set in parameter file */
    145179# define DVO_MAX_PATH 1024
    … …  
    172206double MIN_ERROR;
    173207
     208int    IMFIT_TOO_FEW; // need more than this number of stars to fit an image
    174209int    IMFIT_CLIP_NITER; // number of clipping iterations to perform in FitChip
    175210double IMFIT_CLIP_NSIGMA; // number of sigma to clip in FitChip
    … …  
    188223
    189224int    USE_FIXED_PIXCOORDS;
     225int    USE_GALAXY_MODEL;
     226int    USE_ICRF_CORRECT;
     227
     228int    USE_ICRF_LOCAL;
     229int    USE_ICRF_SHFIT;
     230int    USE_ICRF_POLE;
    190231
    191232int    RESET;
     233int    RESET_IMAGES;
    192234int    RESET_BAD_IMAGES;
    193235int    NLOOP;
    … …  
    204246int    PLOTDELAY;
    205247int    CHIPORDER;
     248int    CHIPMAP;
    206249
    207250int MaxDensityUse;
    … …  
    242285int    CLIP_THRESH;
    243286int USE_BASIC_CHECK;
     287int APPLY_OFFSETS; // in parallel-regions mode, do not launch UpdateObjectOffsets unless -apply-offsets is called
    244288
    245289int ExcludeBogus;
    … …  
    397441void FitSimple (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
    398442void initObjectData (Catalog *catalog, int Ncatalog);
    399 int UpdateObjects (Catalog *catalog, int Ncatalog);
     443int UpdateObjects (Catalog *catalog, int Ncatalog, int Nloop);
    400444int UpdateSimple (Catalog *catalog, int Ncatalog);
    401445int UpdateChips (Catalog *catalog, int Ncatalog);
    … …  
    417461StarData *getImageRef (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode);
    418462StarData *getImageRaw (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode);
     463int setImageRaw (Catalog *catalog, int Ncatalog, off_t im, StarData *raw, off_t Nraw, CoordMode mode);
    419464off_t getImageByID (off_t ID);
     465int updateImageRaw (Catalog *catalog, int Ncatalog, off_t im);
    420466
    421467Mosaic *getmosaics (off_t *N);
    … …  
    447493int relastroGetVisual(void);
    448494
    449 int FixProblemImages (SkyList *skylist);
    450495int *getCatlist (int *N, off_t im);
    451496
    452 void initCoords (void);
    453 void getOffsets (double *dPos, off_t *nPos, off_t N);
    454 void saveOffsets (double dPos, off_t nPos, off_t N);
    455 void setBadCoords (off_t N);
    456 int badCoords (off_t N);
    457 Coords *getCoords (off_t N);
    458 int saveCoords (Coords *coords, off_t N);
    459 void resetImageRaw (Catalog *catalog, int Ncatalog, off_t im);
     497void SaveCoords (Coords *tgt, Coords *src);
     498void RestoreCoords (Coords *tgt, Coords *src, Image *image);
    460499
    461500int high_speed_catalogs (SkyTable *sky, SkyList *skylist, int hostID, char *hostpath);
    … …  
    468507int checkStarMap(int N);
    469508int createStarMap (Catalog *catalog, int Ncatalog);
     509int printStarMap(int N, char *filename);
    470510
    471511int GetScatterRawRef(float *dLsig, float *dMsig, float *dRsig, int *nKeep, StarData *raw, StarData *ref, int Nstars, float SigmaLimit);
    … …  
    519559MeanPos *MeanPosLoad(char *filename, off_t *nmeanpos);
    520560
     561int share_meas_pos (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop);
     562int slurp_meas_pos (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop);
     563int set_meas_pos (MeasPos *measpos, Average *average, MeasureTiny *measure);
     564MeasPos *merge_meas_pos (MeasPos *target, int *ntarget, MeasPos *source, int Nsource);
     565
     566int MeasPosSave(char *filename, MeasPos *measpos, off_t Nmeaspos);
     567MeasPos *MeasPosLoad(char *filename, off_t *nmeaspos);
     568
    521569int indexCatalogs (Catalog *catalog, int Ncatalog);
    522570int catID_and_objID_to_seq (int catID, int objID, int *catSeq, off_t *objSeq);
    … …  
    528576
    529577int share_image_pos (RegionHostTable *regionHosts, int nloop);
    530 int slurp_image_pos (RegionHostTable *regionHosts, int nloop);
     578int slurp_image_pos (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop);
    531579ImagePos *merge_image_pos (ImagePos *target, int *ntarget, ImagePos *source, int Nsource);
    532580int set_image_pos (ImagePos *image_pos, Image *image);
    … …  
    547595int isGPC1stack (int photcode);
    548596int isGPC1warp (int photcode);
     597
     598int save_astrom_table ();
     599AstromOffsetTable *get_astrom_table ();
     600void put_astrom_table (AstromOffsetTable *myTable);
     601
     602int fit_map (AstromOffsetMap *map, StarData *raw, StarData *ref, int Npts);
     603
     604// ICRF QSOs : these must be marked in the database (flag on average, flag on measure)
     605
     606void ICRFinit ();
     607int ICRFsave (int cat, int ave, int meas);
     608int ICRFdata (int n, int *cat, int *ave, int *meas);
     609int ICRFmax ();
     610int select_catalog_ICRF (Catalog *catalog, int Ncatalog);
     611
     612void ICRFobjFree (ICRFobj *icrfobj);
     613ICRFobj *get_ICRF_data (Catalog *catalog, int Ncatalog);
     614ICRFobj *merge_icrf_obj (ICRFobj *target, ICRFobj *source);
     615ICRFobj *slurp_icrf_obj (RegionHostTable *regionHosts, int nloop);
     616int share_icrf_obj (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop);
     617
     618int ICRFobjSave(char *filename, ICRFobj *icrfobj);
     619ICRFobj *ICRFobjLoad(char *filename);
     620
     621int SHfitWithMask (double *R, double *D, double *value, int *mask, int Npts, SHterms *fit);
     622
     623FrameCorrectionType *FrameCorrectionInit (double scale);
     624void FrameCorrectionFree (FrameCorrectionType *frame);
     625
     626int FrameCorrectionParallelMaster (RegionHostTable *regionHosts);
     627int FrameCorrectionParallelSlave (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop);
     628int FrameCorrectionSerial (Catalog *catalog, int Ncatalog);
     629
     630FrameCorrectionSet *FrameCorrectionMeasure (ICRFobj *icrfobj);
     631
     632int FrameCorrectionFromSH (FrameCorrectionType *frame, SHterms *dR, SHterms *dD);
     633
     634Coords *FrameCorrectionImageToMap (Header *header, Matrix *matrix, Coords *coords, int raDirection);
     635int FrameCorrectionMapToImage (Header *header, Matrix *matrix, Coords *coords, int raDirection);
     636int FrameCorrectionSHtoImage (Header *header, Matrix *matrix, FrameCorrectionType *frame, int raDirection);
     637FrameCorrectionType *FrameCorrectionImageToSH (Header *header, Matrix *matrix, FrameCorrectionType *frame, int raDirection);
     638int FrameCorrectionFromSH (FrameCorrectionType *frame, SHterms *dR, SHterms *dD);
     639void FrameCorrectionFree (FrameCorrectionType *frame);
     640FrameCorrectionType *FrameCorrectionInit (double scale);
     641void FrameCorrectionSetFree (FrameCorrectionSet *set);
     642FrameCorrectionSet *FrameCorrectionSetInit ();
     643
     644int FrameCorrectionSetSave(char *filename, FrameCorrectionSet *set);
     645FrameCorrectionSet *FrameCorrectionSetLoad(char *filename);
     646
     647int FrameCorrectionFitSH (ICRFobj *icrfobj, SHterms *dRc, SHterms *dDc);
     648int FrameCorrectionFitLocal (ICRFobj *icrfobj, Coords *coords);
     649int FrameCorrectionApply (Catalog *catalog, int Ncatalog, FrameCorrectionType *frame, Coords *coords);
     650
     651int dump_stardata_pts (StarData *raw, int Npts, char *filename);
     652void printNcatTotal ();
     653
     654int HarvestRegionHosts (RegionHostTable *regionHosts);
     655
     656void lockUpdateChips ();
     657void unlockUpdateChips ();
     658
     659int client_logger_init (char *dirname);
     660int client_logger_message (char *format,...);

  • trunk/Ohana/src/relastro/src/BrightCatalog.c

    r37038 r37807  
    4343  if (!gfits_fread_header (f, &header)) {
    4444    if (VERBOSE) fprintf (stderr, "can't read image subset header\n");
    45     fclose (f);
    4645    goto escape;
    4746  }
    4847  if (!gfits_fread_matrix (f, &matrix, &header)) {
    4948    if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n");
    50     gfits_free_header (&header);
    51     fclose (f);
    5249    goto escape;
    5350  }
    … …  
    7067    GET_COLUMN(R,         "RA",         double);
    7168    GET_COLUMN(D,         "DEC",        double);
     69    GET_COLUMN(RoffGAL,   "ROFF_GAL",   float);
     70    GET_COLUMN(DoffGAL,   "DOFF_GAL",   float);
    7271    GET_COLUMN(M,         "MAG_SYS",    float);
    7372    GET_COLUMN(Mcal,      "MAG_CAL",    float);
    … …  
    9493      measure[i].R         = R[i];
    9594      measure[i].D         = D[i];
     95      measure[i].RoffGAL   = RoffGAL[i];
     96      measure[i].DoffGAL   = DoffGAL[i];
    9697      measure[i].M         = M[i];
    9798      measure[i].Mcal      = Mcal[i];
    … …  
    115116    free (R       );
    116117    free (D       );
     118    free (RoffGAL );
     119    free (DoffGAL );
    117120    free (M       );
    118121    free (Mcal    );
    … …  
    310313  if (catalog) free (catalog);
    311314
    312   fclose (f);
     315  if (f) fclose (f);
    313316  return NULL;
    314317}
    … …  
    343346    gfits_create_table_header (&theader, "BINTABLE", "MEASURE_TINY");
    344347
    345     gfits_define_bintable_column (&theader, "E", "RA_OFF",   "ra offset",                  "arcsec", 1.0, 0.0);
    346     gfits_define_bintable_column (&theader, "E", "DEC_OFF",  "dec offset",                 "arcsec", 1.0, 0.0);
     348    gfits_define_bintable_column (&theader, "D", "RA",       "ra",                         "degrees", 1.0, 0.0);
     349    gfits_define_bintable_column (&theader, "D", "DEC",      "dec",                        "degrees", 1.0, 0.0);
     350    gfits_define_bintable_column (&theader, "E", "ROFF_GAL", "ra offset",                  "arcsec", 1.0, 0.0);
     351    gfits_define_bintable_column (&theader, "E", "DOFF_GAL", "dec offset",                 "arcsec", 1.0, 0.0);
    347352    gfits_define_bintable_column (&theader, "E", "MAG_SYS",  "magnitude (sys)",             NULL,    1.0, 0.0);
    348353    gfits_define_bintable_column (&theader, "E", "MAG_CAL",  "magnitude (cal)",             NULL,    1.0, 0.0);
    … …  
    368373    double *R         ; ALLOCATE (R        ,  double, catalog->Nmeasure);
    369374    double *D         ; ALLOCATE (D        ,  double, catalog->Nmeasure);
     375    float  *RoffGAL   ; ALLOCATE (RoffGAL  ,  float,  catalog->Nmeasure);
     376    float  *DoffGAL   ; ALLOCATE (DoffGAL  ,  float,  catalog->Nmeasure);
    370377    float  *M         ; ALLOCATE (M        ,  float,  catalog->Nmeasure);
    371378    float  *Mcal      ; ALLOCATE (Mcal     ,  float,  catalog->Nmeasure);
    … …  
    390397      R[i]        = measure[i].R        ;
    391398      D[i]        = measure[i].D        ;
     399      RoffGAL[i]  = measure[i].RoffGAL  ;
     400      DoffGAL[i]  = measure[i].DoffGAL  ;
    392401      M[i]        = measure[i].M        ;
    393402      Mcal[i]     = measure[i].Mcal     ;
    … …  
    411420    gfits_set_bintable_column (&theader, &ftable, "RA",         R,         catalog->Nmeasure);
    412421    gfits_set_bintable_column (&theader, &ftable, "DEC",        D,         catalog->Nmeasure);
     422    gfits_set_bintable_column (&theader, &ftable, "ROFF_GAL",   RoffGAL,   catalog->Nmeasure);
     423    gfits_set_bintable_column (&theader, &ftable, "DOFF_GAL",   DoffGAL,   catalog->Nmeasure);
    413424    gfits_set_bintable_column (&theader, &ftable, "MAG_SYS",    M,         catalog->Nmeasure);
    414425    gfits_set_bintable_column (&theader, &ftable, "MAG_CAL",    Mcal,      catalog->Nmeasure);
    … …  
    430441    free (R       );
    431442    free (D       );
     443    free (RoffGAL );
     444    free (DoffGAL );
    432445    free (M       );
    433446    free (Mcal    );

  • trunk/Ohana/src/relastro/src/ConfigInit.c

    r36630 r37807  
    2525  if (!ScanConfig (config, "RELASTRO_IMFIT_CLIP_NSIGMA",   "%lf", 0, &IMFIT_CLIP_NSIGMA))   IMFIT_CLIP_NSIGMA = 3.0;
    2626  if (!ScanConfig (config, "RELASTRO_IMFIT_SYS_SIGMA_LIM", "%lf", 0, &IMFIT_SYS_SIGMA_LIM)) IMFIT_SYS_SIGMA_LIM = 0.01;
     27  if (!ScanConfig (config, "RELASTRO_IMFIT_TOO_FEW",       "%d",  0, &IMFIT_TOO_FEW))       IMFIT_TOO_FEW = 5;
    2728
    2829  if (!ScanConfig (config, "PM_DT_MIN",              "%lf", 0, &PM_DT_MIN))        PM_DT_MIN = 0.25;   
    … …  
    3637
    3738  if (!ScanConfig (config, "USE_FIXED_PIXCOORDS", "%d", 0, &USE_FIXED_PIXCOORDS))  USE_FIXED_PIXCOORDS = FALSE;
     39  if (!ScanConfig (config, "USE_GALAXY_MODEL",    "%d", 0, &USE_GALAXY_MODEL))     USE_GALAXY_MODEL = FALSE;
     40  if (!ScanConfig (config, "USE_ICRF_CORRECT",    "%d", 0, &USE_ICRF_CORRECT))     USE_ICRF_CORRECT = FALSE;
     41  if (!ScanConfig (config, "USE_ICRF_LOCAL",      "%d", 0, &USE_ICRF_LOCAL))       USE_ICRF_LOCAL = FALSE;
     42  if (!ScanConfig (config, "USE_ICRF_SHFIT",      "%d", 0, &USE_ICRF_SHFIT))       USE_ICRF_SHFIT = FALSE;
     43  if (!ScanConfig (config, "USE_ICRF_POLE",       "%d", 0, &USE_ICRF_POLE))        USE_ICRF_POLE = FALSE;
    3844
    3945  // force CATDIR to be absolute (so parallel mode will work)

  • trunk/Ohana/src/relastro/src/FitChip.c

    r37038 r37807  
    1414
    1515  for (Niter = 0; Niter < IMFIT_CLIP_NITER; Niter ++) {
    16    
     16
    1717    // measure the scatter for the unmarked (good) measurements with dM < limit
    1818    // SIGMA_LIM here is redundant with ImageOps.c:915
    … …  
    4949    // use that count to set order_use below
    5050
     51    // 20140925 : the same logic below is used to constrain the dimensions of the image map
     52
    5153    // figures to assess the fitting process:
    5254    // x vs dx, x vs dy, y vs dx, y vs dy :
    5355    // residual vector field
    5456    // dx vs mag, dy vs mag
     57    if (relastroGetVisual()) {
     58      fprintf (stderr, "fitting chip %s, %d matched stars\n", image->name, Nmatch);
     59    }
    5560    relastroVisualPlotChipFit(raw, ref, dRmax, Nmatch);
     61
     62    // for polynomial fits, I used (5,30,60) = (5,6,7) stars per term
    5663
    5764    // set the maximum order for the polynomial (based on number of stars kept above)
    5865    int order_use = 0;
    59     if (Nkeep >  5) order_use = 1; // 4 stars per polynomial term (per dimension)
    60     if (Nkeep > 30) order_use = 2; // 5 stars per polynomial term (per dimension)
    61     if (Nkeep > 60) order_use = 3; // 6 stars per polynomial term (per dimension)
     66    if (Nkeep >   5) order_use = 1; //  5 stars per cell
     67    if (Nkeep >  24) order_use = 2; //  6 stars per cell
     68    if (Nkeep >  63) order_use = 3; //  7 stars per cell
     69    if (Nkeep > 128) order_use = 4; //  8 stars per cell
     70    if (Nkeep > 225) order_use = 5; //  9 stars per cell
     71    if (Nkeep > 360) order_use = 6; // 10 stars per cell
    6272    if (order_use < 1) {
    6373      if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for linear fit\n", Nkeep);
    … …  
    6575      return FALSE;
    6676    }
    67     image[0].coords.Npolyterms = order_use;
     77    if (CHIPMAP) {
     78      order_use = MIN(MIN(order_use, CHIPMAP), 6); // can only go up to 6th order map (can be user limited)
     79    } else {
     80      order_use = MIN(order_use, 3); // can only go up to 3rd order for polynomials
     81    }
     82
     83    if (VERBOSE2) fprintf (stderr, "using %d for %s\n", order_use, image[0].name);
     84
     85    // when fitting the map, first fit a linear model (below? change Npolyterms to -1)
     86    image[0].coords.Npolyterms = CHIPMAP ? 1 : order_use;
    6887
    6988    if (fit) fit_free (fit);
    … …  
    91110    }
    92111
     112    // apply fit to get the fitted X,Y coordinates.  we need these to fit the residual map below
    93113    for (i = 0; i < Nmatch; i++) {
     114      // we have not yet fitted the map, so Npolyterms needs to be 1 here:
     115      LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[0].coords);
     116    }
     117
     118    if (CHIPMAP) {
     119      if (image[0].coords.offsetMap == NULL) {
     120        // allocate a new table and assign to this image
     121        // need to lock threads here because 'table' is a static with
     122        // modifiable elements
     123        lockUpdateChips ();
     124        AstromOffsetTable *table = get_astrom_table ();
     125        AstromOffsetTableNewMap(table, order_use, order_use, image);
     126        unlockUpdateChips ();
     127      }
     128      // XXX Need to check / update the order of the map here.  if we already have a map, it might have been
     129      // allocated with too small arrays
     130      AstromOffsetMapSetOrder (image[0].coords.offsetMap, order_use, order_use, image);
     131      image[0].coords.offsetMap->keep = TRUE; // if we are trying or re-trying this map, we should keep it unless it fails
     132      fit_map (image[0].coords.offsetMap, raw, ref, Nmatch);
     133      image[0].coords.Npolyterms = -1;
     134    }
     135
     136    if (!strcmp(image[0].name, "o5745g0403o.356777.cm.982487.smf[XY17.hdr]")) {
     137      fprintf (stderr, "test image\n");
     138    }
     139
     140    for (i = 0; i < Nmatch; i++) {
     141      // if we have fitted the map above, Npolyterms needs to be -1 here:
    94142      XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[0].coords);
    95143    }
    96 
    97144  }
    98145
    … …  
    223270   L,M -> X,Y (polynomial transformation : WRP)
    224271
    225 Some details about this function:
    226 
    227 - the initial value of the clipping radius is set based on the distribution of the dR
    228   values for the bright sources.
    229 
    230 - NITER clipping passes are performed
    231 
    232 - the per-star astrometric errors are included in the fit.  The photcode table controls
    233   whether only the reported positional errors are used, or if the magnitudes errors are
    234   scaled to determine the error, and if there is a systematic floor for all sources.
    235 
    236 - input detections are pre-filtered on the basis of the photFlags, etc, in bcatalog
    237 
    238 - outlying detections are clipped in the iterative passes, but the clipped detections are
    239   not recorded
     272   Some details about this function:
     273
     274   - the initial value of the clipping radius is set based on the distribution of the dR
     275   values for the bright sources.
     276
     277   - NITER clipping passes are performed
     278
     279   - the per-star astrometric errors are included in the fit.  The photcode table controls
     280   whether only the reported positional errors are used, or if the magnitudes errors are
     281   scaled to determine the error, and if there is a systematic floor for all sources.
     282
     283   - input detections are pre-filtered on the basis of the photFlags, etc, in bcatalog
     284
     285   - outlying detections are clipped in the iterative passes, but the clipped detections are
     286   not recorded
    240287
    241288*/
    … …  
    243290/* example using fit_apply() :
    244291
    245   f = fopen ("test3.dat", "w");
    246 
    247   // apply new coords to raw (X,Y -> L,M)
    248   for (i = 0; i < Nmatch; i++) {
    249     fprintf (f, "%f %f  %f %f  ", raw[i].X, raw[i].Y, raw[i].L, raw[i].M);
    250     fit_apply (newfit, &L1, &M1, raw[i].X - coords[0].crpix1, raw[i].Y - coords[0].crpix2);
    251     fprintf (f, "%f %f\n", L1, M1);
     292   f = fopen ("test3.dat", "w");
     293
     294   // apply new coords to raw (X,Y -> L,M)
     295   for (i = 0; i < Nmatch; i++) {
     296   fprintf (f, "%f %f  %f %f  ", raw[i].X, raw[i].Y, raw[i].L, raw[i].M);
     297   fit_apply (newfit, &L1, &M1, raw[i].X - coords[0].crpix1, raw[i].Y - coords[0].crpix2);
     298   fprintf (f, "%f %f\n", L1, M1);
     299   }
     300   fclose (f);
     301*/
     302
     303int fit_map (AstromOffsetMap *map, StarData *raw, StarData *ref, int Npts) {
     304
     305  // we are actually fitting the residual after the linear fit has been taken off
     306 
     307  // fit the linear terms as above
     308  // calculate dX (raw.X - ref.X) and dY
     309
     310  int i, N;
     311
     312  float *x, *y, *dX, *dY;
     313  ALLOCATE (x,  float, Npts);
     314  ALLOCATE (y,  float, Npts);
     315  ALLOCATE (dX, float, Npts);
     316  ALLOCATE (dY, float, Npts);
     317
     318  N = 0;
     319  for (i = 0; i < Npts; i++) {
     320    if (raw[i].mask) continue;
     321    x[N] = raw[i].X;
     322    y[N] = raw[i].Y;
     323    dX[N] = ref[i].X - raw[i].X;
     324    dY[N] = ref[i].Y - raw[i].Y;
     325    N++;
     326  }
     327
     328  // in coordsops.c:XY_to_LM, the map is defined to carry dX,dY so that:
     329  // (L,M) = f(X',Y') : (X',Y') = (X,Y) + (dX,dY)
     330
     331  AstromOffsetMapFit (map, x, y, dX, N, TRUE);
     332  AstromOffsetMapFit (map, x, y, dY, N, FALSE);
     333
     334  AstromOffsetMapRepair (map, TRUE);
     335  AstromOffsetMapRepair (map, FALSE);
     336
     337  free (x);
     338  free (y);
     339  free (dX);
     340  free (dY);
     341
     342  return TRUE;
     343}
     344
     345int dump_stardata_pts (StarData *raw, int Npts, char *filename) {
     346
     347  FILE *f = fopen (filename, "w");
     348  if (!f) {
     349    fprintf (stderr, "failed to open file %s\n", filename);
     350    return FALSE;
     351  }
     352
     353  int i;
     354  for (i = 0; i < Npts; i++) {
     355    fprintf (f, "%4d %10.6f %10.6f  %10.6f %10.6f  %9.2f %9.2f  %9.2f %9.2f  %9.4f %3d\n",
     356             i, raw[i].R, raw[i].D, raw[i].P, raw[i].Q, raw[i].L, raw[i].M, raw[i].X, raw[i].Y, raw[i].dPos, raw[i].mask);
    252357  }
    253358  fclose (f);
    254 */
    255 
     359
     360  return TRUE;
     361}

  • trunk/Ohana/src/relastro/src/ImageOps.c

    r37038 r37807  
    2828static off_t        *imageIDs; // list of all image IDs
    2929static off_t        *imageIdx; // list of index for image IDs
    30 # else
     30
     31// as an alternative, we generate imageSeq, which directly maps imageID -> seq
     32static off_t        *imageSeq; // list of index for image IDs
     33
     34// MAX_ID requires 512M to store the image index
     35# define MAX_ID      0x8000000
     36static off_t         minImageID = MAX_ID;
     37static off_t         maxImageID = 0;
     38
     39# else /* not using IMAGE_ID */
    3140static unsigned int *start;
    3241static unsigned int *stop;
    … …  
    8493    imageIdx[i] = i;
    8594    imageIDs[i] = image[i].imageID;
     95    minImageID = MIN(minImageID, image[i].imageID);
     96    maxImageID = MAX(maxImageID, image[i].imageID);
     97    myAssert (image[i].imageID < MAX_ID, "image IDs too large for index memory");
    8698  }
    8799  llsortpair (imageIDs, imageIdx, Nimage);
     100
     101  ALLOCATE (imageSeq, off_t, maxImageID + 1);
     102  for (i = 0; i < maxImageID + 1; i++)  {
     103    imageSeq[i] = -1; // not yet assigned
     104  }
     105
     106  for (i = 0; i < Nimage; i++) {
     107    off_t N = image[i].imageID;
     108    myAssert (imageSeq[N] == -1, "previously assigned");
     109    imageSeq[N] = i;
     110  }
     111
    88112# else
    89113  ALLOCATE (start, unsigned, Nimage);
    … …  
    98122
    99123off_t getImageByID (off_t ID) {
     124# if USE_IMAGE_ID
     125
     126  if (imageSeq) {
     127    if (ID < minImageID) return (-1);
     128    if (ID > maxImageID) return (-1);
     129    off_t N = imageSeq[ID];
     130    return N;
     131  }
    100132
    101133  // we have a pair of vectors (imageIDs, imageIdx) sorted by imageIDs
    102134  // use bisection to find the specified image ID
    103135
    104 # if USE_IMAGE_ID
    105136  off_t Nlo, Nhi, N;
    106137
    … …  
    207238    if (VERBOSE2 || (MATCHCAT && (N_onImage[i] < 20))) {
    208239      name = GetPhotcodeNamebyCode (image[i].photcode);
    209       fprintf (stderr, "image "OFF_T_FMT" (%d, %s) has "OFF_T_FMT" of %d measures (%s, %s)",  i,  image[i].imageID, image[i].name, N_onImage[i], image[i].nstar,
     240      int showExample = (Nfew < 30);
     241      if (showExample) fprintf (stderr, "image "OFF_T_FMT" (%d, %s) has "OFF_T_FMT" of %d measures (%s, %s) ",  i,  image[i].imageID, image[i].name, N_onImage[i], image[i].nstar,
    210242               ohana_sec_to_date(image[i].tzero), name);
    211243      if (N_onImage[i] < 20) {
    212         fprintf (stderr, "*");
     244        if (showExample) fprintf (stderr, "*");
    213245        Nfew ++;
    214246      }
    215247      if (N_onImage[i] < 15) {
    216         fprintf (stderr, "*** warning ***");
     248        if (showExample) fprintf (stderr, "*** warning ***");
    217249        Nbad ++;
    218250      }
    219       fprintf (stderr, "\n");
    220     }
    221   }
    222   fprintf (stderr, OFF_T_FMT" total images, %d with < 20 measurements, %d with < 15\n", Nimage, Nfew, Nbad);
     251      if (showExample) fprintf (stderr, "\n");
     252    }
     253  }
     254  if (MATCHCAT) fprintf (stderr, OFF_T_FMT" total images, %d with < 20 measurements, %d with < 15\n", Nimage, Nfew, Nbad);
    223255}
    224256
    … …  
    245277  }
    246278
     279  if (!ID) return; // detection not associated with an image
     280
    247281  idx = getImageByID (ID);
    248282  if (idx == -1) {
    … …  
    250284    return;
    251285  }
     286  measure->myDet = TRUE; // I 'own' this detection (I am calibrating its image)
    252287
    253288  // index for (catalog, measure) -> image
    … …  
    334369# endif
    335370
    336 /*
    337 Coords *getCoords (off_t meas, int cat) {
    338 
    339   off_t i;
    340 
    341   i = MeasureToImage[cat][meas];
    342   if (i == -1) return (NULL);
    343   return (&image[i].coords);
    344 }
    345 */
    346 
    347371float getColorBlue (off_t meas, int cat) {
    348372
    … …  
    433457  Coords *moscoords, *imcoords;
    434458
    435   // check if this image is bad and should be skipped
    436   if (badCoords(im)) return;
    437 
    438459  // WRP images need to have an associated mosaic
    439460  moscoords = NULL;
    … …  
    491512
    492513    int TESTPT = FALSE;
    493 
    494     TESTPT |= CAT_ID_SRC && OBJ_ID_SRC && (measure->imageID == CAT_ID_SRC) && (measure->detID == OBJ_ID_SRC);
    495     TESTPT |= CAT_ID_DST && OBJ_ID_DST && (measure->imageID == CAT_ID_DST) && (measure->detID == OBJ_ID_DST);
     514    TESTPT |= CAT_ID_SRC && OBJ_ID_SRC && (measure->catID == CAT_ID_SRC) && (measure->objID == OBJ_ID_SRC);
     515    TESTPT |= CAT_ID_DST && OBJ_ID_DST && (measure->catID == CAT_ID_DST) && (measure->objID == OBJ_ID_DST);
    496516    if (TESTPT) {
    497517      fprintf (stderr, "got test det\n");
    … …  
    597617  int Noff = NoffRAave + NoffDECave + NoffRAori + NoffDECori;
    598618  if (VERBOSE2 && (Noff > 0)) fprintf (stderr, "Noff ave RA %d, Noff ave DEC %d, Noff ori RA %d, Noff ori DEC %d\n", NoffRAave, NoffDECave, NoffRAori, NoffDECori);
    599   saveOffsets (dPos, nPos, im);
    600619
    601620  return;
    … …  
    604623void printNcatTotal () {
    605624  fprintf (stderr, "NcatTotal: %d\n", NcatTotal);
    606 }
    607 
    608 // return StarData values for detections in the specified image, converting coordinates from the
    609 // chip positions: X,Y -> L,M -> P,Q -> R,D
    610 void resetImageRaw (Catalog *catalog, int Ncatalog, off_t im) {
    611 
    612   off_t i, m, c, n;
    613   double X, Y, L, M, P, Q, R, D;
    614 
    615   Mosaic *mosaic;
    616   Coords *moscoords, *imcoords, *oldcoords;
    617 
    618   fprintf (stderr, "fix resetImageRaw wrt MeasureTiny\n");
    619   abort();
    620 
    621   // check if this image is bad and should be skipped
    622   if (!badCoords(im)) {
    623     fprintf (stderr, "ERROR: inconsistent result?");
    624     abort();
    625   }
    626 
    627   // replace the current coords with the old coords:
    628   oldcoords = getCoords (im);
    629   memcpy (&image[im].coords, oldcoords, sizeof(Coords));
    630 
    631   // WRP images need to have an associated mosaic
    632   moscoords = NULL;
    633   if (!strcmp(&image[im].coords.ctype[4], "-WRP")) {
    634     mosaic = getMosaicForImage (im);
    635     if (mosaic == NULL) return;  // if we cannot find the associated image, skip it
    636     moscoords = &mosaic[0].coords;
    637   }
    638   imcoords = &image[im].coords;
    639 
    640   for (i = 0; i < N_onImage[im]; i++) {
    641     m = ImageToMeasure[im][i];
    642     c = ImageToCatalog[im][i];
    643 
    644     // XXX unclear if this should use Measure or MeasureTiny; it is only called by FixProblemImages, which has not been updated
    645     Measure *measure = &catalog[c].measure[m];
    646 
    647     X = measure[0].Xccd;
    648     Y = measure[0].Yccd;
    649     if (USE_FIXED_PIXCOORDS) {
    650       if (isfinite(measure[0].Xfix) && isfinite(measure[0].Yfix)) {
    651         float dX = measure[0].Xfix - measure[0].Xccd;
    652         float dY = measure[0].Yfix - measure[0].Yccd;
    653         if (hypot(dX,dY) < 2.0) {
    654           X = measure[0].Xfix;
    655           Y = measure[0].Yfix;
    656         }
    657       }
    658     }
    659     n = measure[0].averef;
    660 
    661     if (moscoords == NULL) {
    662       // this is a Simple image (not a mosaic)
    663       // note that for a Simple image, L,M = P,Q
    664       XY_to_LM (&L, &M, X, Y, imcoords);
    665       LM_to_RD (&R, &D, L, M, imcoords);
    666     } else {
    667       XY_to_LM (&L, &M, X, Y, imcoords);
    668       XY_to_LM (&P, &Q, L, M, moscoords);
    669       LM_to_RD (&R, &D, P, Q, moscoords);
    670     }
    671     measure[0].R = R;
    672     measure[0].D = D;
    673 
    674     float dRoff = dvoOffsetR(measure, &catalog[c].average[n]);
    675 
    676     if (dRoff > +180.0*3600.0) {
    677       // average on high end of boundary, move star up
    678       measure[0].R += 360.0;
    679       dRoff -= 360.0*3600.0;
    680     }
    681     if (dRoff < -180.0*3600.0) {
    682       measure[0].R -= 360.0;
    683       dRoff += 360.0*3600.0;
    684     }
    685   }
    686 
    687   return;
    688625}
    689626
    … …  
    736673
    737674    n = measure[0].averef;
     675
     676    int TESTPT = FALSE;
     677    TESTPT |= CAT_ID_SRC && OBJ_ID_SRC && (catalog[c].average[n].catID == CAT_ID_SRC) && (catalog[c].average[n].objID == OBJ_ID_SRC);
     678    TESTPT |= CAT_ID_DST && OBJ_ID_DST && (catalog[c].average[n].catID == CAT_ID_DST) && (catalog[c].average[n].objID == OBJ_ID_DST);
     679    if (TESTPT) {
     680      fprintf (stderr, "got test det\n");
     681    }
    738682
    739683    // an object with only one detection provides no information about the image calibration
    … …  
    808752    ref[i].D = catalog[c].average[n].D;
    809753
     754    // if we are correcting for the Galaxy Motion Model, we assume the mean R,D is at the J2000 epoch position
     755    if (USE_GALAXY_MODEL) {
     756      myAssert (!isnan(measure[0].RoffGAL), "oops");
     757      myAssert (!isnan(measure[0].DoffGAL), "oops");
     758      ref[i].R += measure[0].RoffGAL / 3600.0;
     759      ref[i].D += measure[0].DoffGAL / 3600.0;
     760    }
     761
    810762    ref[i].Mag  = measure[0].M;
    811763    ref[i].dMag = measure[0].dM;
    … …  
    842794  *Nstars = N_onImage[im];
    843795  return (ref);
     796}
     797
     798// return StarData values for detections in the specified image, converting coordinates from the
     799// chip positions: X,Y -> L,M -> P,Q -> R,D.  This function is used by the image fitting steps, for
     800// which the detections have already been filtered when they were loaded (bcatalog)
     801int setImageRaw (Catalog *catalog, int Ncatalog, off_t im, StarData *raw, off_t Nraw, CoordMode mode) {
     802
     803  off_t i, m, c;
     804
     805  Coords *moscoords;
     806
     807  moscoords = NULL;
     808  if (mode == MODE_MOSAIC) {
     809    moscoords = image[im].coords.mosaic;
     810    myAssert (moscoords, "coords.mosaic not defined for image %s (%d)", image[im].name, (int) im);
     811  }
     812
     813  myAssert (Nraw == N_onImage[im], "impossible!");
     814
     815  for (i = 0; i < N_onImage[im]; i++) {
     816    m = ImageToMeasure[im][i];
     817    c = ImageToCatalog[im][i];
     818
     819    // XXX should I use the raw coords or just measure.X,Y -> R,D?
     820
     821    switch (mode) {
     822      case MODE_SIMPLE:
     823        /* note that for a Simple image, L,M = P,Q */
     824        XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords);
     825        raw[i].P = raw[i].L;
     826        raw[i].Q = raw[i].M;
     827        LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, &image[im].coords);
     828        break;
     829      case MODE_MOSAIC:
     830        XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords);
     831        XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].L, raw[i].M, moscoords);
     832        LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, moscoords);
     833        break;
     834      default:
     835        fprintf (stderr, "error: invalid mode in getImageRaw");
     836        abort ();
     837    }
     838
     839    MeasureTiny *measure = &catalog[c].measureT[m];
     840    measure->R = raw[i].R;
     841    measure->D = raw[i].D;
     842    if (catalog[c].measure) {
     843      catalog[c].measure[m].R = raw[i].R;
     844      catalog[c].measure[m].D = raw[i].D;
     845    }     
     846  }
     847  return TRUE;
     848}
     849
     850// return StarData values for detections in the specified image, converting coordinates from the
     851// chip positions: X,Y -> L,M -> P,Q -> R,D.  This function is used by the image fitting steps, for
     852// which the detections have already been filtered when they were loaded (bcatalog)
     853int updateImageRaw (Catalog *catalog, int Ncatalog, off_t im) {
     854
     855  off_t i, m, c;
     856
     857  for (i = 0; i < N_onImage[im]; i++) {
     858    m = ImageToMeasure[im][i];
     859    c = ImageToCatalog[im][i];
     860
     861    MeasureTiny *measure = &catalog[c].measureT[m];
     862
     863    /* apply the current image transformation or use the current value of R+dR, D+dD? */
     864    double X = measure[0].Xccd;
     865    double Y = measure[0].Yccd;
     866    if (USE_FIXED_PIXCOORDS) {
     867      if (isfinite(measure[0].Xfix) && isfinite(measure[0].Yfix)) {
     868        float dX = measure[0].Xfix - measure[0].Xccd;
     869        float dY = measure[0].Yfix - measure[0].Yccd;
     870        if (hypot(dX,dY) < 2.0) {
     871          X = measure[0].Xfix;
     872          Y = measure[0].Yfix;
     873        }
     874      }
     875    }
     876
     877    double R, D;
     878    XY_to_RD (&R, &D, X, Y, &image[im].coords);
     879    measure->R = R;
     880    measure->D = D;
     881    if (catalog[c].measure) {
     882      catalog[c].measure[m].R = R;
     883      catalog[c].measure[m].D = D;
     884    }     
     885  }
     886  return TRUE;
    844887}
    845888

  • trunk/Ohana/src/relastro/src/ImagePosIO.c

    r36630 r37807  
    7070  GET_COLUMN (polyterms,        "POLYTERMS",         float); // verify that we got 14 columns?
    7171  GET_COLUMN (ctype,            "CTYPE",             char);  // verify that we got 15 columns?
    72   GET_COLUMN (Npolyterms,       "NPOLYTERMS",        char);
     72  GET_COLUMN (Npolyterms,       "NPOLYTERMS",        byte);
    7373  GET_COLUMN (dXpixSys,         "XPIX_SYS_ERR",      float);
    7474  GET_COLUMN (dYpixSys,         "YPIX_SYS_ERR",      float);
     75  GET_COLUMN (refColorBlue,     "REF_COLOR_BLUE",    float);
     76  GET_COLUMN (refColorRed,      "REF_COLOR_RED",     float);
    7577  GET_COLUMN (imageID,          "ID",                int);
    7678  GET_COLUMN (nFitAstrom,       "NFIT",              int);
    … …  
    8587  ALLOCATE (image_pos, ImagePos, Nrow);
    8688  for (i = 0; i < Nrow; i++) {
    87     image_pos[i].coords.crval1             = crval1    [i];
    88     image_pos[i].coords.crval2             = crval2    [i];
    89     image_pos[i].coords.crpix1             = crpix1    [i];
    90     image_pos[i].coords.crpix2             = crpix2    [i];
    91     image_pos[i].coords.cdelt1             = cdelt1    [i];
    92     image_pos[i].coords.cdelt2             = cdelt2    [i];
    93     image_pos[i].coords.pc1_1              = pc1_1     [i];
    94     image_pos[i].coords.pc1_2              = pc1_2     [i];
    95     image_pos[i].coords.pc2_1              = pc2_1     [i];
    96     image_pos[i].coords.pc2_2              = pc2_2     [i];
    97     image_pos[i].coords.Npolyterms         = Npolyterms[i];
    98     image_pos[i].dXpixSys                  = dXpixSys  [i];
    99     image_pos[i].dYpixSys                  = dYpixSys  [i];
    100     image_pos[i].imageID                   = imageID   [i];
    101     image_pos[i].nFitAstrom                = nFitAstrom[i];
    102     image_pos[i].flags                     = flags     [i];
     89    InitCoords(&image_pos[i].coords, NULL);
     90    image_pos[i].coords.crval1             = crval1      [i];
     91    image_pos[i].coords.crval2             = crval2      [i];
     92    image_pos[i].coords.crpix1             = crpix1      [i];
     93    image_pos[i].coords.crpix2             = crpix2      [i];
     94    image_pos[i].coords.cdelt1             = cdelt1      [i];
     95    image_pos[i].coords.cdelt2             = cdelt2      [i];
     96    image_pos[i].coords.pc1_1              = pc1_1       [i];
     97    image_pos[i].coords.pc1_2              = pc1_2       [i];
     98    image_pos[i].coords.pc2_1              = pc2_1       [i];
     99    image_pos[i].coords.pc2_2              = pc2_2       [i];
     100    image_pos[i].coords.Npolyterms         = Npolyterms  [i];
     101    image_pos[i].dXpixSys                  = dXpixSys    [i];
     102    image_pos[i].dYpixSys                  = dYpixSys    [i];
     103    image_pos[i].refColorBlue              = refColorBlue[i];
     104    image_pos[i].refColorRed               = refColorRed [i];
     105    image_pos[i].imageID                   = imageID     [i];
     106    image_pos[i].nFitAstrom                = nFitAstrom  [i];
     107    image_pos[i].flags                     = flags       [i];
    103108
    104109    // polyterms and ctype are a bit different
    … …  
    108113  fprintf (stderr, "loaded data for %lld images\n", (long long) Nrow);
    109114
    110   free (crval1    );
    111   free (crval2    );
    112   free (crpix1    );
    113   free (crpix2    );
    114   free (cdelt1    );
    115   free (cdelt2    );
    116   free (pc1_1     );
    117   free (pc1_2     );
    118   free (pc2_1     );
    119   free (pc2_2     );
    120   free (polyterms );
    121   free (ctype     );
    122   free (Npolyterms);
    123   free (dXpixSys  );
    124   free (dYpixSys  );
    125   free (imageID   );
    126   free (nFitAstrom);
    127   free (flags     );
     115  free (crval1      );
     116  free (crval2      );
     117  free (crpix1      );
     118  free (crpix2      );
     119  free (cdelt1      );
     120  free (cdelt2      );
     121  free (pc1_1       );
     122  free (pc1_2       );
     123  free (pc2_1       );
     124  free (pc2_2       );
     125  free (polyterms   );
     126  free (ctype       );
     127  free (Npolyterms  );
     128  free (dXpixSys    );
     129  free (dYpixSys    );
     130  free (refColorBlue);
     131  free (refColorRed );
     132  free (imageID     );
     133  free (nFitAstrom  );
     134  free (flags       );
    128135
    129136  *nimage_pos = Nrow;
    … …  
    153160  gfits_create_table_header (&theader, "BINTABLE", "IMAGE_POS");
    154161
    155   gfits_define_bintable_column (&theader, "D",   "CRVAL1",       "word", "unit", 1.0, 0.0);
    156   gfits_define_bintable_column (&theader, "D",   "CRVAL2",       "word", "unit", 1.0, 0.0);
    157   gfits_define_bintable_column (&theader, "E",   "CRPIX1",       "word", "unit", 1.0, 0.0);
    158   gfits_define_bintable_column (&theader, "E",   "CRPIX2",       "word", "unit", 1.0, 0.0);
    159   gfits_define_bintable_column (&theader, "E",   "CDELT1",       "word", "unit", 1.0, 0.0);
    160   gfits_define_bintable_column (&theader, "E",   "CDELT2",       "word", "unit", 1.0, 0.0);
    161   gfits_define_bintable_column (&theader, "E",   "PC1_1",        "word", "unit", 1.0, 0.0);
    162   gfits_define_bintable_column (&theader, "E",   "PC1_2",        "word", "unit", 1.0, 0.0);
    163   gfits_define_bintable_column (&theader, "E",   "PC2_1",        "word", "unit", 1.0, 0.0);
    164   gfits_define_bintable_column (&theader, "E",   "PC2_2",        "word", "unit", 1.0, 0.0);
    165   gfits_define_bintable_column (&theader, "14E", "POLYTERMS",    "word", "unit", 1.0, 0.0);
    166   gfits_define_bintable_column (&theader, "15A", "CTYPE",        "word", "unit", 1.0, 0.0);
    167   gfits_define_bintable_column (&theader, "L",   "NPOLYTERMS",   "word", "unit", 1.0, 0.0);
    168   gfits_define_bintable_column (&theader, "E",   "XPIX_SYS_ERR", "word", "unit", 1.0, 0.0);
    169   gfits_define_bintable_column (&theader, "E",   "YPIX_SYS_ERR", "word", "unit", 1.0, 0.0);
    170   gfits_define_bintable_column (&theader, "J",   "ID",           "image ID",               "unitless", 1.0, 0.0);
    171   gfits_define_bintable_column (&theader, "J",   "NFIT",         "number of fitted stars", "unitless", 1.0, 0.0);
    172   gfits_define_bintable_column (&theader, "J",   "FLAGS",        "analysis flags",         "unitless", 1.0, 0.0);
     162  gfits_define_bintable_column (&theader, "D",   "CRVAL1",         "word", "unit", 1.0, 0.0);
     163  gfits_define_bintable_column (&theader, "D",   "CRVAL2",         "word", "unit", 1.0, 0.0);
     164  gfits_define_bintable_column (&theader, "E",   "CRPIX1",         "word", "unit", 1.0, 0.0);
     165  gfits_define_bintable_column (&theader, "E",   "CRPIX2",         "word", "unit", 1.0, 0.0);
     166  gfits_define_bintable_column (&theader, "E",   "CDELT1",         "word", "unit", 1.0, 0.0);
     167  gfits_define_bintable_column (&theader, "E",   "CDELT2",         "word", "unit", 1.0, 0.0);
     168  gfits_define_bintable_column (&theader, "E",   "PC1_1",          "word", "unit", 1.0, 0.0);
     169  gfits_define_bintable_column (&theader, "E",   "PC1_2",          "word", "unit", 1.0, 0.0);
     170  gfits_define_bintable_column (&theader, "E",   "PC2_1",          "word", "unit", 1.0, 0.0);
     171  gfits_define_bintable_column (&theader, "E",   "PC2_2",          "word", "unit", 1.0, 0.0);
     172  gfits_define_bintable_column (&theader, "14E", "POLYTERMS",      "word", "unit", 1.0, 0.0);
     173  gfits_define_bintable_column (&theader, "15A", "CTYPE",          "word", "unit", 1.0, 0.0);
     174  gfits_define_bintable_column (&theader, "L",   "NPOLYTERMS",     "word", "unit", 1.0, 0.0);
     175  gfits_define_bintable_column (&theader, "E",   "XPIX_SYS_ERR",   "word", "unit", 1.0, 0.0);
     176  gfits_define_bintable_column (&theader, "E",   "YPIX_SYS_ERR",   "word", "unit", 1.0, 0.0);
     177  gfits_define_bintable_column (&theader, "E",   "REF_COLOR_BLUE", "word", "unit", 1.0, 0.0);
     178  gfits_define_bintable_column (&theader, "E",   "REF_COLOR_RED",  "word", "unit", 1.0, 0.0);
     179  gfits_define_bintable_column (&theader, "J",   "ID",             "image ID",               "unitless", 1.0, 0.0);
     180  gfits_define_bintable_column (&theader, "J",   "NFIT",           "number of fitted stars", "unitless", 1.0, 0.0);
     181  gfits_define_bintable_column (&theader, "J",   "FLAGS",          "analysis flags",         "unitless", 1.0, 0.0);
    173182
    174183  // generate the output array that carries the data
    175184  gfits_create_table (&theader, &ftable);
    176185
    177   double *crval1    ;
    178   double *crval2    ;
    179   float  *crpix1    ;
    180   float  *crpix2    ;
    181   float  *cdelt1    ;
    182   float  *cdelt2    ;
    183   float  *pc1_1     ;
    184   float  *pc1_2     ;
    185   float  *pc2_1     ;
    186   float  *pc2_2     ;
    187   float  *polyterms ;
    188   char   *ctype     ;
    189   char   *Npolyterms;
    190   float  *dXpixSys  ;
    191   float  *dYpixSys  ;
    192   int    *imageID   ;
    193   int    *nFitAstrom;
    194   int    *flags     ;
     186  double *crval1      ;
     187  double *crval2      ;
     188  float  *crpix1      ;
     189  float  *crpix2      ;
     190  float  *cdelt1      ;
     191  float  *cdelt2      ;
     192  float  *pc1_1       ;
     193  float  *pc1_2       ;
     194  float  *pc2_1       ;
     195  float  *pc2_2       ;
     196  float  *polyterms   ;
     197  char   *ctype       ;
     198  char   *Npolyterms  ;
     199  float  *dXpixSys    ;
     200  float  *dYpixSys    ;
     201  float  *refColorBlue;
     202  float  *refColorRed ;
     203  int    *imageID     ;
     204  int    *nFitAstrom  ;
     205  int    *flags       ;
    195206
    196207  // create intermediate storage arrays
    … …  
    210221  ALLOCATE (dXpixSys  ,         float ,          Nimage_pos);
    211222  ALLOCATE (dYpixSys  ,         float ,          Nimage_pos);
     223  ALLOCATE (refColorBlue,       float ,          Nimage_pos);
     224  ALLOCATE (refColorRed ,       float ,          Nimage_pos);
    212225  ALLOCATE (imageID   ,         int   ,          Nimage_pos);
    213226  ALLOCATE (nFitAstrom,         int   ,          Nimage_pos);
    … …  
    216229  // assign the storage arrays
    217230  for (i = 0; i < Nimage_pos; i++) {
    218     crval1    [i] = image_pos[i].coords.crval1    ;
    219     crval2    [i] = image_pos[i].coords.crval2    ;
    220     crpix1    [i] = image_pos[i].coords.crpix1    ;
    221     crpix2    [i] = image_pos[i].coords.crpix2    ;
    222     cdelt1    [i] = image_pos[i].coords.cdelt1    ;
    223     cdelt2    [i] = image_pos[i].coords.cdelt2    ;
    224     pc1_1     [i] = image_pos[i].coords.pc1_1     ;
    225     pc1_2     [i] = image_pos[i].coords.pc1_2     ;
    226     pc2_1     [i] = image_pos[i].coords.pc2_1     ;
    227     pc2_2     [i] = image_pos[i].coords.pc2_2     ;
    228     Npolyterms[i] = image_pos[i].coords.Npolyterms;
    229     dXpixSys  [i] = image_pos[i].dXpixSys         ;
    230     dYpixSys  [i] = image_pos[i].dYpixSys         ;
    231     imageID   [i] = image_pos[i].imageID          ;
    232     nFitAstrom[i] = image_pos[i].nFitAstrom       ;
    233     flags     [i] = image_pos[i].flags            ;
     231    crval1      [i] = image_pos[i].coords.crval1    ;
     232    crval2      [i] = image_pos[i].coords.crval2    ;
     233    crpix1      [i] = image_pos[i].coords.crpix1    ;
     234    crpix2      [i] = image_pos[i].coords.crpix2    ;
     235    cdelt1      [i] = image_pos[i].coords.cdelt1    ;
     236    cdelt2      [i] = image_pos[i].coords.cdelt2    ;
     237    pc1_1       [i] = image_pos[i].coords.pc1_1     ;
     238    pc1_2       [i] = image_pos[i].coords.pc1_2     ;
     239    pc2_1       [i] = image_pos[i].coords.pc2_1     ;
     240    pc2_2       [i] = image_pos[i].coords.pc2_2     ;
     241    Npolyterms  [i] = image_pos[i].coords.Npolyterms;
     242    dXpixSys    [i] = image_pos[i].dXpixSys         ;
     243    dYpixSys    [i] = image_pos[i].dYpixSys         ;
     244    refColorBlue[i] = image_pos[i].refColorBlue     ;
     245    refColorRed [i] = image_pos[i].refColorRed      ;
     246    imageID     [i] = image_pos[i].imageID          ;
     247    nFitAstrom  [i] = image_pos[i].nFitAstrom       ;
     248    flags       [i] = image_pos[i].flags            ;
    234249
    235250    // polyterms and ctype are a bit different
    … …  
    239254
    240255  // add the columns to the output array
    241   gfits_set_bintable_column (&theader, &ftable, "CRVAL1",      crval1    ,         Nimage_pos);
    242   gfits_set_bintable_column (&theader, &ftable, "CRVAL2",      crval2    ,         Nimage_pos);
    243   gfits_set_bintable_column (&theader, &ftable, "CRPIX1",      crpix1    ,         Nimage_pos);
    244   gfits_set_bintable_column (&theader, &ftable, "CRPIX2",      crpix2    ,         Nimage_pos);
    245   gfits_set_bintable_column (&theader, &ftable, "CDELT1",      cdelt1    ,         Nimage_pos);
    246   gfits_set_bintable_column (&theader, &ftable, "CDELT2",      cdelt2    ,         Nimage_pos);
    247   gfits_set_bintable_column (&theader, &ftable, "PC1_1",       pc1_1     ,         Nimage_pos);
    248   gfits_set_bintable_column (&theader, &ftable, "PC1_2",       pc1_2     ,         Nimage_pos);
    249   gfits_set_bintable_column (&theader, &ftable, "PC2_1",       pc2_1     ,         Nimage_pos);
    250   gfits_set_bintable_column (&theader, &ftable, "PC2_2",       pc2_2     ,         Nimage_pos);
    251   gfits_set_bintable_column (&theader, &ftable, "POLYTERMS",   polyterms ,         Nimage_pos);
    252   gfits_set_bintable_column (&theader, &ftable, "CTYPE",       ctype     ,         Nimage_pos);
    253   gfits_set_bintable_column (&theader, &ftable, "NPOLYTERMS",  Npolyterms,         Nimage_pos);
    254   gfits_set_bintable_column (&theader, &ftable, "XPIX_SYS_ERR",dXpixSys  ,         Nimage_pos);
    255   gfits_set_bintable_column (&theader, &ftable, "YPIX_SYS_ERR",dYpixSys  ,         Nimage_pos);
    256   gfits_set_bintable_column (&theader, &ftable, "ID",          imageID   ,         Nimage_pos);
    257   gfits_set_bintable_column (&theader, &ftable, "NFIT",        nFitAstrom,         Nimage_pos);
    258   gfits_set_bintable_column (&theader, &ftable, "FLAGS",       flags     ,         Nimage_pos);
    259 
    260   free (crval1    );
    261   free (crval2    );
    262   free (crpix1    );
    263   free (crpix2    );
    264   free (cdelt1    );
    265   free (cdelt2    );
    266   free (pc1_1     );
    267   free (pc1_2     );
    268   free (pc2_1     );
    269   free (pc2_2     );
    270   free (polyterms );
    271   free (ctype     );
    272   free (Npolyterms);
    273   free (dXpixSys  );
    274   free (dYpixSys  );
    275   free (imageID   );
    276   free (nFitAstrom);
    277   free (flags     );
     256  gfits_set_bintable_column (&theader, &ftable, "CRVAL1",         crval1    ,         Nimage_pos);
     257  gfits_set_bintable_column (&theader, &ftable, "CRVAL2",         crval2    ,         Nimage_pos);
     258  gfits_set_bintable_column (&theader, &ftable, "CRPIX1",         crpix1    ,         Nimage_pos);
     259  gfits_set_bintable_column (&theader, &ftable, "CRPIX2",         crpix2    ,         Nimage_pos);
     260  gfits_set_bintable_column (&theader, &ftable, "CDELT1",         cdelt1    ,         Nimage_pos);
     261  gfits_set_bintable_column (&theader, &ftable, "CDELT2",         cdelt2    ,         Nimage_pos);
     262  gfits_set_bintable_column (&theader, &ftable, "PC1_1",          pc1_1     ,         Nimage_pos);
     263  gfits_set_bintable_column (&theader, &ftable, "PC1_2",          pc1_2     ,         Nimage_pos);
     264  gfits_set_bintable_column (&theader, &ftable, "PC2_1",          pc2_1     ,         Nimage_pos);
     265  gfits_set_bintable_column (&theader, &ftable, "PC2_2",          pc2_2     ,         Nimage_pos);
     266  gfits_set_bintable_column (&theader, &ftable, "POLYTERMS",      polyterms ,         Nimage_pos);
     267  gfits_set_bintable_column (&theader, &ftable, "CTYPE",          ctype     ,         Nimage_pos);
     268  gfits_set_bintable_column (&theader, &ftable, "NPOLYTERMS",     Npolyterms,         Nimage_pos);
     269  gfits_set_bintable_column (&theader, &ftable, "XPIX_SYS_ERR",   dXpixSys  ,         Nimage_pos);
     270  gfits_set_bintable_column (&theader, &ftable, "YPIX_SYS_ERR",   dYpixSys  ,         Nimage_pos);
     271  gfits_set_bintable_column (&theader, &ftable, "REF_COLOR_BLUE", refColorBlue,       Nimage_pos);
     272  gfits_set_bintable_column (&theader, &ftable, "REF_COLOR_RED",  refColorRed ,       Nimage_pos);
     273  gfits_set_bintable_column (&theader, &ftable, "ID",             imageID   ,         Nimage_pos);
     274  gfits_set_bintable_column (&theader, &ftable, "NFIT",           nFitAstrom,         Nimage_pos);
     275  gfits_set_bintable_column (&theader, &ftable, "FLAGS",          flags     ,         Nimage_pos);
     276
     277  free (crval1      );
     278  free (crval2      );
     279  free (crpix1      );
     280  free (crpix2      );
     281  free (cdelt1      );
     282  free (cdelt2      );
     283  free (pc1_1       );
     284  free (pc1_2       );
     285  free (pc2_1       );
     286  free (pc2_2       );
     287  free (polyterms   );
     288  free (ctype       );
     289  free (Npolyterms  );
     290  free (dXpixSys    );
     291  free (dYpixSys    );
     292  free (refColorBlue);
     293  free (refColorRed );
     294  free (imageID     );
     295  free (nFitAstrom  );
     296  free (flags       );
    278297 
    279298  FILE *f = fopen (filename, "w");

  • trunk/Ohana/src/relastro/src/StarMaps.c

    r36630 r37807  
    176176}
    177177
     178int printStarMap(int N, char *filename) {
     179
     180  off_t i, Nimages;
     181  Image *images = getimages(&Nimages, NULL);
     182  double L, M;
     183
     184  FILE *f = fopen (filename, "w");
     185
     186  for (i = 0; i < starmap[N].Npoints; i++) {
     187    // set the transformed coordinates
     188    XY_to_LM (&L, &M, starmap[N].points[i].X, starmap[N].points[i].Y, &images[N].coords);
     189
     190    fprintf (f, "%d %7.2f %f : %7.2f %7.2f : %7.2f %7.2f : %7.2f %7.2f\n", (int) i, starmap[N].points[i].X, starmap[N].points[i].Y, starmap[N].points[i].L, starmap[N].points[i].M, L, M, starmap[N].points[i].L - L, starmap[N].points[i].M - M);
     191  }
     192  fclose (f);
     193  return (TRUE);
     194}
     195

  • trunk/Ohana/src/relastro/src/UpdateChips.c

    r33652 r37807  
    3636  off_t i, Nraw, Nref, nFitAstr;
    3737  StarData *raw, *ref;
    38   Coords *oldCoords;
    3938  float dXpixSys, dYpixSys;
    4039  double *Ro, *Do;
    … …  
    7473      Nskip ++;
    7574      mode[i] = 0;
     75      continue;
     76    }
     77    if (Nraw <= IMFIT_TOO_FEW) {
     78      Nskip ++;
     79      mode[i] = 0;
     80      free (raw);
    7681      continue;
    7782    }
    … …  
    9095
    9196    // save these in case of failure
    92     saveCoords (&image[i].coords, i);
     97    Coords oldCoords;
     98    SaveCoords (&oldCoords, &image[i].coords);
     99
    93100    dXpixSys = image[i].dXpixSys;
    94101    dYpixSys = image[i].dYpixSys;
    … …  
    100107      if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT" of %d\n", image[i].name,  i,  Nraw, image[i].nstar);
    101108
    102       // restore status quo ante
    103       oldCoords = getCoords (i);
    104       memcpy (&image[i].coords, oldCoords, sizeof(Coords));
    105       image[i].dXpixSys = dXpixSys;
    106       image[i].dYpixSys = dYpixSys;
    107       image[i].nFitAstrom = nFitAstr;
    108 
    109       saveCenter (image, &Ro[i], &Do[i], i);
    110       mode[i] = 1;
    111       NoldFit ++;
    112       free (raw);
    113       free (ref);
    114       continue;
     109      if (1) {
     110        // restore status quo ante (replace truMap with tmpMap)
     111        RestoreCoords (&image[i].coords, &oldCoords, &image[i]);
     112        image[i].dXpixSys = dXpixSys;
     113        image[i].dYpixSys = dYpixSys;
     114        image[i].nFitAstrom = nFitAstr;
     115
     116        saveCenter (image, &Ro[i], &Do[i], i);
     117        mode[i] = 1;
     118        NoldFit ++;
     119        free (raw);
     120        free (ref);
     121        continue;
     122      }
    115123    }
    116124
    117125    if (!checkStarMap (i)) {
    118126      if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name,  i,  Nraw);
    119       // restore status quo ante
    120       oldCoords = getCoords (i);
    121       memcpy (&image[i].coords, oldCoords, sizeof(Coords));
    122       image[i].dXpixSys = dXpixSys;
    123       image[i].dYpixSys = dYpixSys;
    124       image[i].nFitAstrom = nFitAstr;
    125 
    126       saveCenter (image, &Ro[i], &Do[i], i);
    127       mode[i] = 2;
    128       image[i].flags |= ID_IMAGE_ASTROM_POOR;
    129       NoldFit ++;
    130       free (raw);
    131       free (ref);
    132       continue;
     127
     128      if (1) {
     129        // restore status quo ante (replace truMap with tmpMap)
     130        RestoreCoords (&image[i].coords, &oldCoords, &image[i]);
     131        image[i].dXpixSys = dXpixSys;
     132        image[i].dYpixSys = dYpixSys;
     133        image[i].nFitAstrom = nFitAstr;
     134
     135        saveCenter (image, &Ro[i], &Do[i], i);
     136        mode[i] = 2;
     137        image[i].flags |= ID_IMAGE_ASTROM_POOR;
     138        NoldFit ++;
     139        free (raw);
     140        free (ref);
     141        continue;
     142      }
    133143    }
     144
     145    // Apply the modified coords back to the measure.R,D.  Note that raw.R,D, ref.L,M, etc
     146    // are all automatically updated in this block because they are re-generated from
     147    // image.coords on each pass.
     148    setImageRaw (catalog, Ncatalog, i, raw, Nraw, MODE_MOSAIC);
     149    if (USE_GALAXY_MODEL) {
     150      image[i].flags |= ID_IMAGE_ASTROM_GMM;
     151    }
    134152
    135153    saveCenter (image, &Ro[i], &Do[i], i);
    … …  
    235253  off_t Nraw, Nref, nFitAstr;
    236254  StarData *raw, *ref;
    237   Coords *oldCoords;
    238255  float dXpixSys, dYpixSys;
    239256
    … …  
    260277      threadinfo->Nskip ++;
    261278      threadinfo->mode[i] = 0;
     279      continue;
     280    }
     281    if (Nraw <= IMFIT_TOO_FEW) {
     282      threadinfo->Nskip ++;
     283      threadinfo->mode[i] = 0;
     284      free (raw);
    262285      continue;
    263286    }
    … …  
    276299
    277300    // save these in case of failure
    278     saveCoords (&image[i].coords, i);
     301    Coords oldCoords;
     302    SaveCoords (&oldCoords, &image[i].coords);
    279303    dXpixSys = image[i].dXpixSys;
    280304    dYpixSys = image[i].dYpixSys;
    … …  
    286310      if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT" of %d\n", image[i].name,  i,  Nraw, image[i].nstar);
    287311
    288       // restore status quo ante
    289       oldCoords = getCoords (i);
    290       memcpy (&image[i].coords, oldCoords, sizeof(Coords));
     312      // restore status quo ante (replace truMap with tmpMap)
     313      RestoreCoords (&image[i].coords, &oldCoords, &image[i]);
    291314      image[i].dXpixSys = dXpixSys;
    292315      image[i].dYpixSys = dYpixSys;
    … …  
    303326    if (!checkStarMap (i)) {
    304327      if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name,  i,  Nraw);
    305       // restore status quo ante
    306       oldCoords = getCoords (i);
    307       memcpy (&image[i].coords, oldCoords, sizeof(Coords));
     328
     329      // restore status quo ante (replace truMap with tmpMap)
     330      RestoreCoords (&image[i].coords, &oldCoords, &image[i]);
    308331      image[i].dXpixSys = dXpixSys;
    309332      image[i].dYpixSys = dYpixSys;
    … …  
    319342    }
    320343
     344    // apply the modified R,D back to the measures
     345    setImageRaw (threadinfo->catalog, threadinfo->Ncatalog, i, raw, Nraw, MODE_MOSAIC);
     346
    321347    saveCenter (image, &threadinfo->Ro[i], &threadinfo->Do[i], i);
    322348    threadinfo->mode[i] = 3;
    … …  
    333359static pthread_mutex_t UpdateChips_mutex = PTHREAD_MUTEX_INITIALIZER;
    334360
     361void lockUpdateChips () {
     362  pthread_mutex_lock (&UpdateChips_mutex);
     363}
     364void unlockUpdateChips () {
     365  pthread_mutex_unlock (&UpdateChips_mutex);
     366}
     367
    335368// we have an array of chips (image, Nimage).  we need to hand out images one at a time to
    336369// the worker threads as they need
    337370off_t getNextImageForThread () {
    338371
    339   pthread_mutex_lock (&UpdateChips_mutex);
     372  lockUpdateChips ();
    340373  if (nextImage >= Nimage) {
    341     pthread_mutex_unlock (&UpdateChips_mutex);
     374    unlockUpdateChips ();
    342375    return (-1);
    343376  }
    … …  
    345378  nextImage ++;
    346379
    347   pthread_mutex_unlock (&UpdateChips_mutex);
     380  unlockUpdateChips ();
    348381  return (thisImage);
    349382}
    … …  
    494527}
    495528
    496 // XXX if (!FindMosaicForImage (image, Nimage, i)) { }
     529// save the coords and offset map in src at tgt
     530void SaveCoords (Coords *tgt, Coords *src) {
     531  myAssert (tgt, "oops");
     532  myAssert (src, "oops");
     533  CopyCoords (tgt, src); // src retains a pointer to the AstromOffsetMap table, to be saved
     534  tgt->offsetMap = AstromOffsetMapCopy(src->offsetMap); // oldCoords now saves the old values in a new structure
     535}
     536
     537// restore the coords and offset map at tgt from src
     538void RestoreCoords (Coords *tgt, Coords *src, Image *image) {
     539
     540  AstromOffsetMap *truMap = tgt->offsetMap;
     541  AstromOffsetMap *tmpMap = src->offsetMap;
     542
     543  // copy the coords structure data first, then fix the offsetMap pointers
     544  // this modifies tgt.coords, but we need to keep the original pointer
     545  CopyCoords (tgt, src);
     546  tgt->offsetMap = truMap;
     547
     548  // we want to keep the old solution, which is (maybe) a pointer to a stand-alone
     549  // map.  the image has a bad solution, but it is a pointer to the map in the I/O
     550  // table.  those values will be saved on exit.
     551
     552  if (!tmpMap && !truMap) return;
     553
     554  if (tmpMap && truMap) {
     555    AstromOffsetMapSetOrder (truMap, tmpMap->Nx, tmpMap->Ny, image);
     556    AstromOffsetMapCopyData (truMap, tmpMap);
     557    return;
     558  }
     559
     560  if (tmpMap && !truMap) {
     561    // is this even possible?
     562    lockUpdateChips ();
     563    AstromOffsetTable *table = get_astrom_table ();
     564    AstromOffsetTableNewMap(table, tmpMap->Nx, tmpMap->Ny, image);  // registers the map with the image
     565    unlockUpdateChips ();
     566    AstromOffsetMapCopyData (image->coords.offsetMap, tmpMap);
     567    return;
     568  }
     569  if (!tmpMap && truMap) {
     570    truMap->keep = FALSE;
     571    return;
     572  }
     573  myAbort ("oops");
     574}

  • trunk/Ohana/src/relastro/src/UpdateMeasures.c

    r35763 r37807  
    11# include "relastro.h"
     2void dump_measures (Average *average, Measure *measure); // in ImageOps.c
    23
    3 void printNcatTotal ();
     4// this function is called by UpdateObjectOffsets or relastro_images.c for reset
     5int UpdateMeasures (Catalog *catalog, int Ncatalog) {
     6
     7  off_t i, j, Nimage;
     8  Image *image;
     9  double DPOS_MAX_ASEC;
     10
     11  int badImage =
     12    ID_IMAGE_ASTROM_NOCAL |
     13    ID_IMAGE_ASTROM_POOR |
     14    ID_IMAGE_ASTROM_FAIL |
     15    ID_IMAGE_ASTROM_SKIP |
     16    ID_IMAGE_ASTROM_FEW;
     17
     18  if (RESET_BAD_IMAGES) badImage = 0;
     19
     20  image = getimages (&Nimage, NULL);
     21
     22  // track measurements which are far from their original position (poor at pole)
     23  int NoffRAave = 0;  int NoffRAori = 0;
     24  int NoffDECave = 0; int NoffDECori = 0;
     25
     26  for (i = 0; i < Ncatalog; i++) {
     27    myAssert (!catalog[i].Nmeasure || catalog[i].measureT, "programming error");
     28    for (j = 0; j < catalog[i].Nmeasure; j++) {
     29      MeasureTiny *measureT = &catalog[i].measureT[j];
     30      Measure     *measureB = catalog[i].measure ? &catalog[i].measure[j] : NULL;
     31
     32      off_t im = getImageByID (measureT->imageID);
     33      if (im < 0) continue; // detections without imageIDs are not associated with a known image (eg, refs)
     34
     35      // check that we have the right image:
     36      myAssert (measureT->photcode == image[im].photcode, "image photcode mismatch");
     37      myAssert (measureT->t >= image[im].tzero, "image time mismatch (1)");
     38      myAssert (measureT->t <= image[im].tzero + image[im].NX*image[im].trate/5000.0, "image time mismatch (2)");
     39      // I'm allowing the trate*NX to be a factor of 2x too small
     40
     41      // skip measurements on images that have failed solutions (divergent or otherwise)
     42      // XXX apply this or not??
     43      if (FALSE && (image[im].flags & badImage)) continue;
     44
     45      Coords *moscoords = image[im].coords.mosaic;
     46      Coords *imcoords = &image[im].coords;
     47
     48      if (moscoords) {
     49        DPOS_MAX_ASEC = 3600.0*DPOS_MAX*hypot(moscoords->cdelt1, moscoords->cdelt2);
     50      } else {
     51        DPOS_MAX_ASEC = 3600.0*DPOS_MAX*hypot(imcoords->cdelt1, imcoords->cdelt2);
     52      }
     53
     54      double X = measureT->Xccd;
     55      double Y = measureT->Yccd;
     56      if (USE_FIXED_PIXCOORDS) {
     57        if (isfinite(measureT->Xfix) && isfinite(measureT->Yfix)) {
     58          float dX = measureT->Xfix - measureT->Xccd;
     59          float dY = measureT->Yfix - measureT->Yccd;
     60          if (hypot(dX,dY) < 2.0) {
     61            X = measureT->Xfix;
     62            Y = measureT->Yfix;
     63          }
     64        }
     65      }
     66
     67      double R, D;
     68      XY_to_RD (&R, &D, X, Y, imcoords);
     69
     70      int n = measureT[0].averef;
     71      Average *average = &catalog[i].average[n];
     72
     73      double dR = 3600.0*(average[0].R - R);
     74      double dD = 3600.0*(average[0].D - D);
     75
     76      // make sure detection is on the same side of the 0,360 boundary as the average
     77      // this will give some funny results withing ~1 arcsec of the pol
     78      if (dR > +180.0*3600.0) {
     79        // average on high end of boundary, move star up
     80        R += 360.0;
     81        dR = 3600.0*(average[0].R - R);
     82      }
     83      if (dR < -180.0*3600.0) {
     84        // average on low end of boundary, move star down
     85        R -= 360.0;
     86        dR = 3600.0*(average[0].R - R);
     87      }
     88
     89      float csdec = cos(average[0].D * RAD_DEG);
     90
     91      // complain if the new location is far from the average location
     92      // NOTE: This should never happen, or our StarMap tests are not working
     93      if (fabs(dR*csdec) > 3.0*ADDSTAR_RADIUS) {
     94        NoffRAave ++;
     95        if (VERBOSE2 && catalog[i].measure) {
     96          fprintf (stderr, "measurement is far from average location (R): %f %f (%f %f %f)\n", average[0].R, average[0].D, dR, dR*csdec, dD);
     97          dump_measures (&average[0], catalog[i].measure);
     98        }
     99      }
     100      if (fabs(dD) > 3.0*ADDSTAR_RADIUS) {
     101        NoffDECave ++;
     102        if (VERBOSE2 && catalog[i].measure) {
     103          fprintf (stderr, "measurement is far from average location (D): %f %f (%f %f)\n", average[0].R, average[0].D, dR, dD);
     104          dump_measures (&average[0], catalog[i].measure);
     105        }
     106      }
     107
     108      // complain if the new location is far from the old location
     109      if (fabs(csdec*(measureT[0].R - R)) > DPOS_MAX_ASEC) {
     110        NoffRAori ++;
     111        if (VERBOSE2 && catalog[i].measure) {
     112          fprintf (stderr, "measurement moves far from original location (R): %f %f (%f %f %f %f)\n", average[0].R, average[0].D, measureT[0].R, dR, csdec*(measureT[0].R - R), dD);
     113          dump_measures (&average[0], catalog[i].measure);
     114        }
     115      }
     116      if (fabs(measureT[0].D - D) > DPOS_MAX_ASEC) {
     117        NoffDECori ++;
     118        if (VERBOSE2 && catalog[i].measure) {
     119          fprintf (stderr, "measurement moves far from original location (D): %f %f (%f %f %f)\n", average[0].R, average[0].D, dR, measureT[0].D, dD);
     120          dump_measures (&average[0], catalog[i].measure);
     121        }
     122      }
     123      measureT->R = R;
     124      measureT->D = D;
     125      if (measureB) {
     126        measureB->R = R;
     127        measureB->D = D;
     128      }
     129    }
     130  }
     131  // printNcatTotal ();
     132
     133  int Noff = NoffRAave + NoffDECave + NoffRAori + NoffDECori;
     134  if (VERBOSE && (Noff > 0)) fprintf (stderr, "Noff ave RA %d, Noff ave DEC %d, Noff ori RA %d, Noff ori DEC %d\n", NoffRAave, NoffDECave, NoffRAori, NoffDECori);
     135
     136  return (TRUE);
     137}
    4138
    5139// this function operates on Measure, not MeasureTiny
    6 int UpdateMeasures (Catalog *catalog, int Ncatalog) {
     140int UpdateMeasuresOld (Catalog *catalog, int Ncatalog) {
    7141
    8142  off_t i, Nimage;
    … …  
    32166  }
    33167
    34   printNcatTotal ();
     168  // printNcatTotal ();
    35169
    36170  return (TRUE);

  • trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c

    r37038 r37807  
    5959
    6060    // uses MeasureTiny values, but will also update Measure values when populated
    61     UpdateObjects (&catalog, 1);
     61    // if Galaxy Motion Model correction is active, do not apply at this last stage
     62    UpdateObjects (&catalog, 1, FALSE);
    6263
    6364    free_tiny_values(&catalog);
    … …  
    7475    save_catalogs (&catalog, 1);
    7576  }
    76  
     77  printNcatTotal();
    7778  return (TRUE);
    7879}
     80
     81int UpdateObjectOffsets_parallel_group (HostTableGroup *group, SkyList *sky);
     82int UpdateObjectOffsets_parallel_table (HostTable *table, SkyList *sky);
    7983
    8084// CATDIR is supplied globally
    … …  
    9094    exit (1);
    9195  }   
     96
     97# if (0)
     98 
     99  UpdateObjectOffsets_parallel_table (table, sky);
     100
     101# else
     102
     103  int Ngroups;
     104  HostTableGroup *groups = HostTableGroups (table, &Ngroups);
     105  // split the table into Ngroups, each with a unique set of machines (this avoids overloading the remote host machines)
     106
     107  int i;
     108  for (i = 0; i < Ngroups; i++) {
     109    // update only a group of unique machines at a time
     110    UpdateObjectOffsets_parallel_group (&groups[i], sky);
     111  }
     112# endif
     113
     114  return TRUE;
     115}     
     116
     117// CATDIR is supplied globally
     118# define DEBUG 1
     119int UpdateObjectOffsets_parallel_group (HostTableGroup *group, SkyList *sky) {
     120
     121  // launch the relastro_client jobs to the parallel hosts
     122
     123  int i, j;
     124  for (i = 0; i < group->Nhosts; i++) {
     125
     126    if (sky->Nregions < group->Nhosts) {
     127      // do any of the regions want this host?
     128      int wantThisHost = FALSE;
     129      for (j = 0; j < sky->Nregions; j++) {
     130        if (HostTableTestHost (sky->regions[j], group->hosts[i][0].hostID)) {
     131          wantThisHost = TRUE;
     132          break;
     133        }
     134      }
     135      if (!wantThisHost) {
     136        // fprintf (stderr, "skip host %s\n", group->hosts[i][0].hostname);
     137        continue;
     138      }
     139    }
     140
     141    // ensure that the paths are absolute path names
     142    char *tmppath = abspath (group->hosts[i][0].pathname, DVO_MAX_PATH);
     143    free (group->hosts[i][0].pathname);
     144    group->hosts[i][0].pathname = tmppath;
     145
     146    // options / arguments that can affect relastro_client -load:
     147    // VERBOSE, VERBOSE2
     148    // RESET (-reset)
     149    // TimeSelect -time
     150    // (note that psfQF is applied rigidly at 0.85, as is the galaxy test)
     151    // ImagSelect, ImagMin, ImagMax
     152    // MaxDensityUse, MaxDensityValue
     153
     154    char command[1024];
     155    snprintf (command, 1024, "relastro_client -update-offsets -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s -minerror %f",
     156              group->hosts[i][0].hostID, CATDIR, group->hosts[i][0].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, MIN_ERROR);
     157
     158    if (FIT_MODE == FIT_PM_ONLY)         strextend (command, "-pm");
     159    if (FIT_MODE == FIT_PAR_ONLY)        strextend (command, "-par");
     160    if (FIT_MODE == FIT_PM_AND_PAR)      strextend (command, "-pmpar");
     161
     162    if (VERBOSE)       strextend (command, "-v");
     163    if (VERBOSE2)      strextend (command, "-vv");
     164    if (RESET)         strextend (command, "-reset");
     165    if (UPDATE)        strextend (command, "-update");
     166
     167    if (RESET_BAD_IMAGES) strextend (command, "-reset-bad-images");
     168
     169    if (ImagSelect)    strextend (command, "-instmag %f %f", ImagMin, ImagMax);
     170    if (MaxDensityUse) strextend (command, "-max-density %f", MaxDensityValue);
     171   
     172    if (USE_BASIC_CHECK) strextend (command, "-basic-image-search");
     173    if (FlagOutlier)     strextend (command, "-clip %d", CLIP_THRESH);
     174    if (ExcludeBogus)    strextend (command, "-exclude-bogus %f", ExcludeBogusRadius);
     175   
     176    if (USE_FIXED_PIXCOORDS) strextend (command, "-D USE_FIXED_PIXCOORDS 1");
     177
     178    if (PHOTCODE_KEEP_LIST) strextend (command, "+photcode %s", PHOTCODE_KEEP_LIST);
     179    if (PHOTCODE_SKIP_LIST) strextend (command, "-photcode %s", PHOTCODE_SKIP_LIST);
     180    if (PhotFlagSelect)     strextend (command, "+photflags");
     181    if (PhotFlagBad)        strextend (command, "+photflagbad %d", PhotFlagBad);
     182    if (PhotFlagPoor)       strextend (command, "+photflagpoor %d", PhotFlagPoor);
     183    // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values
     184
     185    if (DCR_BLUE_COLOR_POS && DCR_BLUE_COLOR_NEG) {
     186      strextend (command, "-dcr-blue-color %s %s", DCR_BLUE_COLOR_POS, DCR_BLUE_COLOR_NEG);
     187    }
     188    if (DCR_RED_COLOR_POS && DCR_RED_COLOR_NEG) {
     189      strextend (command, "-dcr-red-color %s %s", DCR_RED_COLOR_POS, DCR_RED_COLOR_NEG);
     190    }
     191
     192    if (MinBadQF > 0.0)        strextend (command, "-min-bad-psfqf %f", MinBadQF);
     193    if (MaxMeanOffset != 10.0) strextend (command, "-max-mean-offset  %f", MaxMeanOffset);
     194
     195    if (TimeSelect) {
     196      char *tstart = ohana_sec_to_date (TSTART);
     197      char *tstop  = ohana_sec_to_date (TSTOP);
     198      strextend (command, "-time %s %s", tstart, tstop);
     199      free (tstart);
     200      free (tstop);
     201    }
     202    fprintf (stderr, "command: %s\n", command);
     203
     204    if (PARALLEL_MANUAL) continue;
     205
     206    if (PARALLEL_SERIAL) {
     207      int status = system (command);
     208      if (status) {
     209        fprintf (stderr, "ERROR running relastro_client\n");
     210        exit (2);
     211      }
     212    } else {
     213      // launch the job on the remote machine (no handshake)
     214      int errorInfo = 0;
     215      int pid = rconnect ("ssh", group->hosts[i][0].hostname, command, group->hosts[i][0].stdio, &errorInfo, FALSE);
     216      if (!pid) {
     217        if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", group->hosts[i][0].hostname, errorInfo);
     218        exit (1);
     219      }
     220      group->hosts[i][0].pid = pid; // save for future reference
     221    }
     222  }
     223
     224  if (PARALLEL_MANUAL) {
     225    fprintf (stderr, "run the relastro_client commands above.  when these are done, hit return\n");
     226    getchar();
     227  }
     228  if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) {
     229    HostTableGroupWaitJobsGetIO (group, __FILE__, __LINE__, VERBOSE);
     230  }
     231
     232  return TRUE;
     233}     
     234
     235// CATDIR is supplied globally
     236# define DEBUG 1
     237int UpdateObjectOffsets_parallel_table (HostTable *table, SkyList *sky) {
     238
     239  // launch the relastro_client jobs to the parallel hosts
    92240
    93241  int i, j;

  • trunk/Ohana/src/relastro/src/UpdateObjects.c

    r37261 r37807  
    22# define PAR_TOOFEW 5
    33
    4 int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m);
     4int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m, int applyGalaxyOffset);
    55int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats);
    66
    … …  
    6969
    7070  /* project coordinates to a plane centered on the object with units of arcsec */
    71   coords.crval1 = 0;
    72   coords.crval2 = 0;
    73   coords.crpix1 = 0;
    74   coords.crpix2 = 0;
     71  InitCoords (&coords, "DEC--SIN");
    7572  coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
    76   coords.pc1_1  = coords.pc2_2 = 1.0;
    77   coords.pc1_2  = coords.pc2_1 = 0.0;
    78   coords.Npolyterms = 1;
    79   strcpy (coords.ctype, "DEC--SIN");
    8073
    8174  // use J2000 as a reference time
    … …  
    10699// This function operates on both Measure and MeasureTiny.  In the big stages, this should
    107100// be called with just MeasureTiny set and Measure == NULL
    108 int UpdateObjects (Catalog *catalog, int Ncatalog) {
     101int UpdateObjects (Catalog *catalog, int Ncatalog, int Nloop) {
    109102
    110103  initObjectData (catalog, Ncatalog);
    … …  
    137130
    138131      UpdateObjects_Stack(average, secfilt, measure, measureBig, Nsecfilt, &fitStatsStack);
    139       UpdateObjects_Chips(average, secfilt, measure, measureBig, Nsecfilt, &fitStatsChips, i, m);
    140     }
    141     if (VERBOSE) fprintf (stderr, "catalog %d : chips "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n",  i,  fitStatsChips.Nave,  fitStatsChips.Npm,  fitStatsChips.Npar,  fitStatsChips.Nskip, fitStatsChips.Noffset);
    142     if (VERBOSE) fprintf (stderr, "catalog %d : stack "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n",  i,  fitStatsStack.Nave,  fitStatsStack.Npm,  fitStatsStack.Npar,  fitStatsStack.Nskip, fitStatsStack.Noffset);
     132      UpdateObjects_Chips(average, secfilt, measure, measureBig, Nsecfilt, &fitStatsChips, i, m, Nloop);
     133    }
     134    if (VERBOSE2) fprintf (stderr, "catalog %d : chips "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n",  i,  fitStatsChips.Nave,  fitStatsChips.Npm,  fitStatsChips.Npar,  fitStatsChips.Nskip, fitStatsChips.Noffset);
     135    if (VERBOSE2) fprintf (stderr, "catalog %d : stack "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n",  i,  fitStatsStack.Nave,  fitStatsStack.Npm,  fitStatsStack.Npar,  fitStatsStack.Nskip, fitStatsStack.Noffset);
    143136    sumFitStats (&fitStatsChips, &sumStatsChips);
    144137    sumFitStats (&fitStatsStack, &sumStatsStack);
    … …  
    146139  freeObjectData ();
    147140
    148   if (VERBOSE) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n",  (sumStatsChips.Nave + sumStatsChips.Npm + sumStatsChips.Npar),  sumStatsChips.Nave,  sumStatsChips.Npm,  sumStatsChips.Npar,  sumStatsChips.Nskip, sumStatsChips.Noffset);
    149   if (VERBOSE) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n",  (sumStatsStack.Nave + sumStatsStack.Npm + sumStatsStack.Npar),  sumStatsStack.Nave,  sumStatsStack.Npm,  sumStatsStack.Npar,  sumStatsStack.Nskip, sumStatsStack.Noffset);
     141  if (VERBOSE && (Ncatalog > 1)) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n",  (sumStatsChips.Nave + sumStatsChips.Npm + sumStatsChips.Npar),  sumStatsChips.Nave,  sumStatsChips.Npm,  sumStatsChips.Npar,  sumStatsChips.Nskip, sumStatsChips.Noffset);
     142  if (VERBOSE && (Ncatalog > 1)) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n",  (sumStatsStack.Nave + sumStatsStack.Npm + sumStatsStack.Npar),  sumStatsStack.Nave,  sumStatsStack.Npm,  sumStatsStack.Npar,  sumStatsStack.Nskip, sumStatsStack.Noffset);
    150143  return (TRUE);
    151144}
    … …  
    153146// This function operates on both Measure and MeasureTiny.  In the big stages, this should
    154147// be called with just MeasureTiny set and Measure == NULL
    155 int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m) {
     148int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m, int applyGalaxyOffset) {
    156149
    157150  int setRefColor = areImagesMatched();
    … …  
    215208    double Ri = getMeanR (&measure[k], average, secfilt);
    216209    double Di = getMeanD (&measure[k], average, secfilt);
     210
     211    // if we are correcting for the Galaxy Motion Model, only should apply it here
     212    // (a) when we are working to correct the images (mean R,D assumed to be at J2000) and
     213    // (b) if we think the measure R,D is already at the image epoch position
     214    if (USE_GALAXY_MODEL && applyGalaxyOffset) {
     215      Ri -= measure[k].RoffGAL / 3600.0;
     216      Di -= measure[k].DoffGAL / 3600.0;
     217    }
     218
     219    // XXX add in dR,dD GAL here
    217220
    218221    // mark (as POOR) any measurements which are deviant from the mean by > ExcludeBogusRadius
    … …  
    329332  // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq
    330333
     334  // project all of the R,D coordinates to a plane centered on this coordinate. set
     335  // the times to be relative to Tmean (this is required for parallax as well)
     336  for (k = 0; k < N; k++) {
     337    RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords);
     338    T[k] -= Tmean;
     339    if (XVERB) {
     340      fprintf (stderr, OFF_T_FMT" %f %f %f  %f %f +/- %f %f\n",  k, T[k], R[k], D[k], X[k], Y[k], dX[k], dY[k]);
     341    }
     342  }       
     343
    331344  // *** first fit for the proper motion (skip fit if Trange or Npts is too small) ***
    332345  if ((mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) {
    … …  
    340353    }
    341354
    342     // project all of the R,D coordinates to a plane centered on this coordinate. set
    343     // the times to be relative to Tmean (this is required for parallax as well)
    344     for (k = 0; k < N; k++) {
    345       RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords);
    346       T[k] -= Tmean;
    347       if (XVERB) {
    348         fprintf (stderr, OFF_T_FMT" %f %f %f  %f %f +/- %f %f\n",  k, T[k], R[k], D[k], X[k], Y[k], dX[k], dY[k]);
    349       }
    350     }     
    351 
    352355    FitPM (&fitPM, X, dX, Y, dY, T, N, XVERB);
    353356
    … …  
    357360    XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &coords);
    358361    if (XVERB) fprintf (stderr, "project: %f %f : %f %f : %f\n", fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p);
     362    if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
     363    // XXX : does this make sense at 0,360 boundary?
    359364
    360365    fitPM.p  = fitPM.dp  = 0.0;
    … …  
    409414    fitStats->Npar ++;
    410415
     416    if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
     417
    411418    // XXX a hard-wired hack...
    412419    if ((fabs(fitPAR.uR) > 2.0) || (fabs(fitPAR.uD) > 2.0)) {
    … …  
    419426    // ALWAYS fit the average model
    420427    StatType statsR, statsD;
    421     liststats_pos (R, dR, N, &statsR, XVERB); // WARNING: this function modifies R (do not use after here)
    422     liststats_pos (D, dD, N, &statsD, XVERB); // WARNING: this function modifies D (do not use after here)
    423 
    424     fitAve.Ro = statsR.mean;
    425     fitAve.dRo = 3600.0*statsR.sigma;
    426 
    427     fitAve.Do = statsD.mean;
    428     fitAve.dDo = 3600.0*statsD.sigma;
     428    liststats_pos (X, dX, N, &statsR, XVERB); // WARNING: this function modifies R (do not use after here)
     429    liststats_pos (Y, dY, N, &statsD, XVERB); // WARNING: this function modifies D (do not use after here)
     430
     431    // project Ro, Do back to RA,DEC
     432    XY_to_RD (&fitAve.Ro, &fitAve.Do, statsR.mean, statsD.mean, &coords);
     433    if (XVERB) fprintf (stderr, "average: %f %f\n", fitAve.Ro, fitAve.Do);
     434
     435    fitAve.dRo = statsR.sigma;
     436    fitAve.dDo = statsD.sigma;
    429437
    430438    fitAve.chisq = (N > 1) ? 0.5 * (statsR.chisq + statsD.chisq) : NAN;
    … …  
    466474      break;
    467475  }
    468   if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) plx=(%f +/- %f)\n",
    469                       average[0].R,
    470                       average[0].D,
    471                       fit.Ro, fit.Do,
    472                       3600*(average[0].R - fit.Ro),
    473                       3600*(average[0].D - fit.Do),
    474                       fit.uR, fit.uD, fit.p, fit.dp);
     476  if (XVERB) {
     477    fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) plx=(%f +/- %f)\n",
     478             average[0].R,
     479             average[0].D,
     480             fit.Ro, fit.Do,
     481             3600*(average[0].R - fit.Ro),
     482             3600*(average[0].D - fit.Do),
     483             fit.uR, fit.uD, fit.p, fit.dp);
     484    fprintf (stderr, "-----\n");
     485  }
    475486
    476487  // make sure that the fit succeeded
    … …  
    494505  coords.crval1 = average[0].R;
    495506  coords.crval2 = average[0].D;
     507  if (isnan(coords.crval1)) {
     508    return (FALSE);
     509  }
     510  if (isnan(coords.crval2)) {
     511    return (FALSE);
     512  }
    496513
    497514  double dXoff, dYoff;
    … …  
    505522    return FALSE; // XXX ??
    506523  }
    507 
    508   average[0].R          = fit.Ro; // RA in degrees
    509   average[0].D          = fit.Do; // DEC in degrees
    510   average[0].dR         = fit.dRo; // RA scatter in arcsec
    511   average[0].dD         = fit.dDo; // DEC scatter in arcsec
    512 
    513   average[0].uR         = fit.uR; // RA proper motion in arcsec/year
    514   average[0].uD         = fit.uD; // DEC proper motion in arcsec/year
    515   average[0].duR        = fit.duR; // RA proper motion error in arcsec/year
    516   average[0].duD        = fit.duD; // DEC proper motion error in arcsec/year
    517 
    518   average[0].P          = fit.p; // parallax in arcsec
    519   average[0].dP         = fit.dp; // parallax error in arcsec
    520 
    521   average[0].ChiSqAve   = fitAve.chisq;
    522   average[0].ChiSqPM    = fitPM.chisq;
    523   average[0].ChiSqPar   = fitPAR.chisq;
    524 
    525   average[0].Tmean      = (Tmean * 86400 * 365.25) + T2000;
    526   average[0].Trange     = (Trange * 86400 * 365.25);
    527   average[0].Npos       = fit.Nfit;
    528 
    529   // unset the NO_ASTROM bit (not(NO_ASTROM) == HAVE_ASTROM)
    530   average[0].flags &= ~ID_STAR_NO_ASTROM;
    531524
    532525  if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) chisq=(%f, %f, %f)\n",
    … …  
    540533                      fitAve.chisq, fitPM.chisq, fitPAR.chisq);
    541534
     535  average[0].R          = fit.Ro; // RA in degrees
     536  average[0].D          = fit.Do; // DEC in degrees
     537  average[0].dR         = fit.dRo; // RA scatter in arcsec
     538  average[0].dD         = fit.dDo; // DEC scatter in arcsec
     539
     540  average[0].uR         = fit.uR; // RA proper motion in arcsec/year
     541  average[0].uD         = fit.uD; // DEC proper motion in arcsec/year
     542  average[0].duR        = fit.duR; // RA proper motion error in arcsec/year
     543  average[0].duD        = fit.duD; // DEC proper motion error in arcsec/year
     544
     545  average[0].P          = fit.p; // parallax in arcsec
     546  average[0].dP         = fit.dp; // parallax error in arcsec
     547
     548  average[0].ChiSqAve   = fitAve.chisq;
     549  average[0].ChiSqPM    = fitPM.chisq;
     550  average[0].ChiSqPar   = fitPAR.chisq;
     551
     552  average[0].Tmean      = (Tmean * 86400 * 365.25) + T2000;
     553  average[0].Trange     = (Trange * 86400 * 365.25);
     554  average[0].Npos       = fit.Nfit;
     555
     556  // unset the NO_ASTROM bit (not(NO_ASTROM) == HAVE_ASTROM)
     557  average[0].flags &= ~ID_STAR_NO_ASTROM;
     558
    542559  return (TRUE);
    543560}

  • trunk/Ohana/src/relastro/src/args.c

    r37038 r37807  
    265265    remove_argument (N, &argc, argv);
    266266    MaxDensityUse = TRUE;
     267  }
     268
     269  APPLY_OFFSETS = FALSE;
     270  if ((N = get_argument (argc, argv, "-apply-offsets"))) {
     271    remove_argument (N, &argc, argv);
     272    APPLY_OFFSETS = TRUE;
    267273  }
    268274
    … …  
    372378  }
    373379
     380  CHIPMAP = 0;
     381  if ((N = get_argument (argc, argv, "-chipmap"))) {
     382    remove_argument (N, &argc, argv);
     383    CHIPMAP = atoi(argv[N]);
     384    remove_argument (N, &argc, argv);
     385
     386  }
     387
    374388  SAVEPLOT = FALSE;
    375389  PLOTSTUFF = FALSE;
    … …  
    399413    RESET = TRUE;
    400414  }
    401 
    402415  RESET_BAD_IMAGES = FALSE;
    403416  if ((N = get_argument (argc, argv, "-reset-bad-images"))) {
    404417    remove_argument (N, &argc, argv);
    405418    RESET_BAD_IMAGES = TRUE;
     419  }
     420  RESET_IMAGES = FALSE;
     421  if ((N = get_argument (argc, argv, "-reset-images"))) {
     422    remove_argument (N, &argc, argv);
     423    RESET_IMAGES = TRUE;
    406424  }
    407425
    … …  
    868886  fprintf (stderr, "  -high-speed (code[,code,code]) (code[,code,code]) (radius) (output catdir)\n");
    869887  fprintf (stderr, "  -hpm (radius) (output catdir)\n");
     888  fprintf (stderr, "  -chipmap (MaxOrder)\n");
    870889  fprintf (stderr, "  -merge-source (objID) (catID) into (objID) (catID)\n\n");
    871890
    … …  
    880899  fprintf (stderr, "  -statmode (mode)\n");
    881900  fprintf (stderr, "  -reset\n");
     901  fprintf (stderr, "  -reset-images\n");
    882902  fprintf (stderr, "  -nloop (N) : number of image-fit iterations\n");
    883903  fprintf (stderr, "  -update : apply new fit to database\n");

  • trunk/Ohana/src/relastro/src/assign_images.c

    r36630 r37807  
    2323  BuildChipMatch (image, Nimage);
    2424  MARKTIME("build chip match for %d images: %f sec\n", (int) Nimage, dtime);
     25
     26  char mapfile[DVO_MAX_PATH];
     27  snprintf (mapfile, DVO_MAX_PATH, "%s/AstroMap.fits", CATDIR);
     28  AstromOffsetTable *table = AstromOffsetMapLoad (mapfile, VERBOSE);
     29
     30  // assign images.coords.offsetMap -> table->map[i]
     31  if (table) {
     32    AstromOffsetTableMatchChips (image, Nimage, table);
     33  } else {
     34    table = AstromOffsetTableInit ();
     35  }
     36  put_astrom_table (table);
    2537
    2638  initMosaics (image, Nimage);
    … …  
    89101    if (!strcmp(&image[j].coords.ctype[4], "-DIS")) continue;
    90102
    91     // match the image to its corresponding ASTROMETRIC mosaic (not the same as the Mosaic above)
    92     if (!FindMosaicForImage (image, Nimage, j)) {
    93       if (VERBOSE2) fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i);
    94       continue;
    95     }
    96    
    97103    // Exclude images with crazy astrometry
    98104    // XXX NOTE : this is gpc1-specific
    … …  
    150156
    151157    // this is a bit memory expensive : I am making a complete copy of the image table here
    152     off_t Nsubset = regionHosts->hosts[i].Nimage;
    153     regionHosts->hosts[i].image[Nsubset] = image[j];
    154     regionHosts->hosts[i].imseq[Nsubset] = j;
    155 
    156     regionHosts->hosts[i].Nimage ++;
    157     if (regionHosts->hosts[i].Nimage == regionHosts->hosts[i].NIMAGE) {
    158       regionHosts->hosts[i].NIMAGE += D_NIMAGE;
    159       REALLOCATE (regionHosts->hosts[i].image, Image, regionHosts->hosts[i].NIMAGE);
    160       REALLOCATE (regionHosts->hosts[i].imseq, off_t, regionHosts->hosts[i].NIMAGE);
     158    off_t Nsubset = host->Nimage;
     159    host->image[Nsubset] = image[j];
     160    host->imseq[Nsubset] = j;
     161
     162    host->Nimage ++;
     163    if (host->Nimage == host->NIMAGE) {
     164      host->NIMAGE += D_NIMAGE;
     165      REALLOCATE (host->image, Image, host->NIMAGE);
     166      REALLOCATE (host->imseq, off_t, host->NIMAGE);
     167    }
     168
     169    // save the astrometry maps, where they exist
     170    if (image[j].coords.offsetMap) {
     171      if (!host->astromTable) {
     172        host->astromTable = AstromOffsetTableInit();
     173      }
     174      AstromOffsetTableAddMapFromImage(host->astromTable, &image[j]);
    161175    }
    162176  }
    … …  
    199213}
    200214
    201 int calculate_host_image_bounds (RegionHostTable *regionHosts) {
    202 
    203   int i, n;
    204   off_t j;
    205 
    206   for (i = 0; i < regionHosts->Nhosts; i++) {
    207 
    208     RegionHostInfo *host = &regionHosts->hosts[i];
    209 
    210     // XXX clear the chip match?
    211 
    212     BuildChipMatch (host->image, host->Nimage);
    213 
    214     double Rmin = 360.0;
    215     double Rmax =   0.0;
    216     double Dmin = +90.0;
    217     double Dmax = -90.0;
    218 
    219     for (j = 0; j < host->Nimage; j++) {
    220 
    221       Image *image = &host->image[j];
    222 
    223       if (!FindMosaicForImage (host->image, host->Nimage, j)) {
    224         fprintf (stderr, "missing astrometry? programming error?\n");
    225         abort ();
    226       }
    227      
    228       // define image corners
    229       for (n = 0; n < 4; n++) {
    230         double Xc, Yc, Rc, Dc;
    231         Xc = Xf[n]*image->NX;
    232         Yc = Yf[n]*image->NY;
    233         XY_to_RD (&Rc, &Dc, Xc, Yc, &image->coords);
    234         Rc = ohana_normalize_angle_to_midpoint (Rc, 180.0);
    235      
    236         Rmin = MIN (Rmin, Rc);
    237         Dmin = MIN (Dmin, Rc);
    238        
    239         Rmax = MAX (Rmax, Rc);
    240         Dmax = MAX (Dmax, Rc);
    241       }
    242 
    243     }
    244 
    245     host->RminCat = Rmin;
    246     host->DminCat = Dmin;
    247     host->RmaxCat = Rmax;
    248     host->DmaxCat = Dmax;
    249 
    250     // regionHosts needs to have the full outer boundary
    251     // (so reload_catalogs covers the correct region)
    252     regionHosts->Rmin = MIN(Rmin, regionHosts->Rmin);
    253     regionHosts->Rmax = MAX(Rmax, regionHosts->Rmax);
    254     regionHosts->Dmin = MIN(Dmin, regionHosts->Dmin);
    255     regionHosts->Dmax = MAX(Dmax, regionHosts->Dmax);
    256   }
    257   return TRUE;
    258 }
    259 
    260 
    261215// XXX add a search tree to speed this up?
    262216int find_host_for_coords (RegionHostTable *regionHosts, double Rc, double Dc) {

  • trunk/Ohana/src/relastro/src/bcatalog.c

    r36833 r37807  
    2727
    2828  /* for outlier rejection, project coordinates to a plane centered on the object with units of arcsec */
    29   coords.crval1 = 0;
    30   coords.crval2 = 0;
    31   coords.crpix1 = 0;
    32   coords.crpix2 = 0;
     29  InitCoords (&coords, "DEC--SIN");
    3330  coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
    34   coords.pc1_1  = coords.pc2_2 = 1.0;
    35   coords.pc1_2  = coords.pc2_1 = 0.0;
    36   coords.Npolyterms = 1;
    37   strcpy (coords.ctype, "DEC--SIN");
    3831
    3932  // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
    … …  
    7265    if (catalog[0].average[i].Nmeasure <= SRC_MEAS_TOOFEW) {
    7366      continue;
    74       // XXX can't (easily) count this for a given chip
    7567    }
    7668   
    … …  
    8476    if (RESET) {
    8577      subcatalog[0].average[Naverage].flags &= ~averageBits;
     78    }
     79
     80    // if we want to do static correction to ICRF, keep the ICRF QSOs (regardless of other data quality)
     81    if (USE_ICRF_CORRECT && (subcatalog[0].average[Naverage].flags & ID_ICRF_QSO)) {
     82      Nm = 0;
     83      for (j = 0; j < catalog[0].average[i].Nmeasure; j++) {
     84         
     85        offset = catalog[0].average[i].measureOffset + j;
     86         
     87        CopyMeasureToTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]);
     88        subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_ASTROM;
     89        subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL;
     90        subcatalog[0].measureT[Nmeasure].averef   = Naverage;
     91        if (RESET) {
     92          subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM;
     93          subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_AREA;
     94        }
     95         
     96        Nmeasure ++;
     97        Nm ++;
     98        if (Nmeasure == NMEASURE) {
     99          NMEASURE += 1000;
     100          REALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE);
     101        }
     102      }
     103      subcatalog[0].average[Naverage].Nmeasure = Nm;
     104      Naverage ++;
     105      if (Naverage == NAVERAGE) {
     106        NAVERAGE += 50;
     107        REALLOCATE (subcatalog[0].average, Average, NAVERAGE);
     108        REALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt);
     109      }
     110      continue;
     111    }
     112
     113    // if we want to correct to the galay model, we need to use objects with distances determined by Green & Schlafly
     114    if (USE_GALAXY_MODEL && !catalog[0].average[i].Nstarpar) {
     115      // we may want to add some starpar quality filters...
     116      continue;
    86117    }
    87118
    … …  
    116147      // exclude bogus
    117148      if (ExcludeBogus) {
    118           double Ri = getMeanR_Big (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]);
    119           double Di = getMeanD_Big (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]);
    120           coords.crval1 = catalog[0].average[i].R;
    121           coords.crval2 = catalog[0].average[i].D;
    122           double Xi, Yi;
    123           RD_to_XY (&Xi, &Yi, Ri, Di, &coords);
    124           double radius = hypot(Xi, Yi);
    125           if (radius > ExcludeBogusRadius) {
    126               NskipBogus ++;
    127               if (VERBOSE2) {
    128                 FILE *foutput = fbogus ? fbogus : stderr;
    129                 char *date = ohana_sec_to_date(catalog[0].measure[offset].t);
    130                 fprintf (foutput, "exclude bogus: %10.6f %10.6f : %10.6f %10.6f : %6.2f %6.2f : %6.2f : %5d %s\n", catalog[0].average[i].R, catalog[0].average[i].D, Ri, Di, Xi, Yi, radius, catalog[0].measure[offset].photcode, date);
    131                 free (date);
    132               }
    133               continue;
     149        double Ri = getMeanR_Big (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]);
     150        double Di = getMeanD_Big (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]);
     151        coords.crval1 = catalog[0].average[i].R;
     152        coords.crval2 = catalog[0].average[i].D;
     153        double Xi, Yi;
     154        RD_to_XY (&Xi, &Yi, Ri, Di, &coords);
     155        double radius = hypot(Xi, Yi);
     156        if (radius > ExcludeBogusRadius) {
     157          NskipBogus ++;
     158          if (VERBOSE2) {
     159            FILE *foutput = fbogus ? fbogus : stderr;
     160            char *date = ohana_sec_to_date(catalog[0].measure[offset].t);
     161            fprintf (foutput, "exclude bogus: %10.6f %10.6f : %10.6f %10.6f : %6.2f %6.2f : %6.2f : %5d %s\n", catalog[0].average[i].R, catalog[0].average[i].D, Ri, Di, Xi, Yi, radius, catalog[0].measure[offset].photcode, date);
     162            free (date);
    134163          }
     164          continue;
     165        }
    135166      }
    136167
    … …  
    157188      subcatalog[0].measureT[Nmeasure].averef   = Naverage;
    158189      if (RESET) {
    159         subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM;
    160         subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_AREA;
     190        subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM;
     191        subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_AREA;
    161192      }
    162193
    … …  
    164195      Nm ++;
    165196      if (Nmeasure == NMEASURE) {
    166         NMEASURE += 1000;
    167         REALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE);
    168       }
     197        NMEASURE += 1000;
     198        REALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE);
     199      }
     200    }
     201    // we have filtered this object.  do not keep it if we now have too few meas
     202    if (Nm <= SRC_MEAS_TOOFEW) {
     203      Nmeasure -= Nm;
     204      continue;
    169205    }
    170206    subcatalog[0].average[Naverage].Nmeasure = Nm;
    … …  
    186222  assert (Nsecfilt == catalog[0].Nsecfilt);
    187223
    188   // limit the total number of stars in the catalog
     224// limit the total number of stars in the catalog
    189225  if (MaxDensityUse) {
    190226    LimitDensityCatalog_ByNmeasureGrid (subcatalog, catalog);
    … …  
    217253
    218254# define SWAPFUNC(A,B){ int tmpI; off_t tmpT;   \
    219   tmpI = X[A]; X[A] = X[B]; X[B] = tmpI; \
    220   tmpT = Y[A]; Y[A] = Y[B]; Y[B] = tmpT; \
    221 }
     255    tmpI = X[A]; X[A] = X[B]; X[B] = tmpI;      \
     256    tmpT = Y[A]; Y[A] = Y[B]; Y[B] = tmpT;      \
     257  }
    222258# define COMPARE(A,B)(X[A] > X[B])
    223259
    … …  
    231267/* generate a grid in a locally-projected coordinate system, try to select average entries
    232268   from each grid cell in decending Nmeasure order.
    233  */
     269*/
    234270int LimitDensityCatalog_ByNmeasureGrid (Catalog *subcatalog, Catalog *oldcatalog) {
    235271
    … …  
    269305  /* project coordinates to a plane centered on the object with units of arcsec */
    270306  Coords coords;
     307  InitCoords (&coords, "DEC--SIN");
    271308  coords.crval1 = Rc;
    272309  coords.crval2 = Dc;
    273   coords.crpix1 = 0;
    274   coords.crpix2 = 0;
    275310  coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
    276   coords.pc1_1  = coords.pc2_2 = 1.0;
    277   coords.pc1_2  = coords.pc2_1 = 0.0;
    278   coords.Npolyterms = 1;
    279   strcpy (coords.ctype, "DEC--SIN");
    280311
    281312  // convert all average R,D values to X,Y:

  • trunk/Ohana/src/relastro/src/fitpoly.c

    r33652 r37807  
    7878}
    7979 
     80int array_print (double **matrix, int Nx, int Ny) {
     81
     82  int i, j;
     83  for (i = 0; i < Nx; i++) {
     84    for (j = 0; j < Ny; j++) {
     85      fprintf (stderr, "%10.3e ", matrix[i][j]);
     86    }
     87    fprintf (stderr, "\n");
     88  }
     89  return TRUE;
     90}
     91
     92int fit_print (CoordFit *fit) {
     93
     94  fprintf (stderr, "Npts: %d, Norder: %d, Nterms: %d, Nsums: %d, Nelems: %d\n",
     95           fit[0].Npts, fit[0].Norder, fit[0].Nterms, fit[0].Nsums, fit[0].Nelems);
     96
     97  fprintf (stderr, "sums: \n");
     98  array_print (fit[0].sum, fit[0].Nsums, fit[0].Nsums);
     99
     100  fprintf (stderr, "xsums: \n");
     101  array_print (fit[0].xsum, fit[0].Nterms, fit[0].Nterms);
     102
     103  fprintf (stderr, "ysums: \n");
     104  array_print (fit[0].ysum, fit[0].Nterms, fit[0].Nterms);
     105
     106  fprintf (stderr, "xfits: \n");
     107  array_print (fit[0].xfit, fit[0].Nterms, fit[0].Nterms);
     108
     109  fprintf (stderr, "yfits: \n");
     110  array_print (fit[0].yfit, fit[0].Nterms, fit[0].Nterms);
     111 
     112  return TRUE;
     113}
     114
    80115// XXX use implicit masks as below or explicit masks (with function to set?)
    81116// XXX eg, add a global mask to this file and
    … …  
    134169      for (j = 0; j < fit[0].Nelems; j++) {
    135170        matrix[i][j] = 0.0;
     171        matrix[j][i] = 0.0;
    136172      }
    137173      vector[i][0] = 0.0;
    … …  
    148184  }     
    149185 
    150   if (!dgaussjordan (matrix, vector, fit[0].Nelems, 2)) {
     186  if (!dgaussjordan_pivot (matrix, vector, fit[0].Nelems, 2, 1e-16)) {
    151187    array_free (matrix, fit[0].Nelems);
    152188    array_free (vector, fit[0].Nelems);
    … …  
    277313/* this should only apply to the polynomial, not the projection terms */
    278314/* compare with psastro supporting code */
     315// this code will work for linear (Npolyterm == 0 or 1) and linear + map fits
    279316int fit_apply_coords (CoordFit *fit, Coords *coords, int keepRef) {
    280317
    … …  
    352389  /* keep the order and type from initial values */
    353390 
     391  if (isnan(coords[0].crval1)) {
     392    return FALSE;
     393  }
     394  if (isnan(coords[0].crval2)) {
     395    return FALSE;
     396  }
     397  if (isnan(coords[0].crpix1)) {
     398    return FALSE;
     399  }
     400  if (isnan(coords[0].crpix2)) {
     401    return FALSE;
     402  }
     403
    354404  return (TRUE);
    355405}

  • trunk/Ohana/src/relastro/src/high_speed_objects.c

    r37038 r37807  
    149149
    150150  // define a local projection
     151  InitCoords (&tcoords, "DEC--ARC");
    151152  tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax);
    152153  if (region[0].Dmax < 90) {
    … …  
    155156    tcoords.crval2 = 90.0;
    156157  }
    157   tcoords.crpix1 = 0;
    158   tcoords.crpix2 = 0;
    159158  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
    160   tcoords.pc1_1 = tcoords.pc2_2 = 1.0;
    161   tcoords.pc1_2 = tcoords.pc2_1 = 0.0;
    162   tcoords.Npolyterms = 1;
    163   strcpy (tcoords.ctype, "DEC--ARC");
    164159
    165160  /* build spatial index (RA sort) referencing input array sequence */
    … …  
    271266
    272267  populate_tiny_values (&catalogOut, DVO_TV_MEASURE);
    273   UpdateObjects (&catalogOut, 1);
     268  UpdateObjects (&catalogOut, 1, 0);
    274269  free_tiny_values(&catalogOut);
    275270

  • trunk/Ohana/src/relastro/src/high_speed_utils.c

    r33652 r37807  
    141141    off_t n;
    142142    for (n = 0; n < NfieldsA; n++) {
    143         valuesA[n] = dbExtractAverages (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], &catalog[0].measure[m], &fieldsA[n]);
     143      valuesA[n] = dbExtractAverages (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], &catalog[0].measure[m], NULL, NULL, &fieldsA[n]);
    144144    }
    145145    return dbBooleanCond(stackA, NstackA, valuesA);
    … …  
    155155    off_t n;
    156156    for (n = 0; n < NfieldsB; n++) {
    157         valuesB[n] = dbExtractAverages (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], &catalog[0].measure[m], &fieldsB[n]);
     157        valuesB[n] = dbExtractAverages (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], &catalog[0].measure[m], NULL, NULL, &fieldsB[n]);
    158158    }
    159159    return dbBooleanCond(stackB, NstackB, valuesB);

  • trunk/Ohana/src/relastro/src/hpm_objects.c

    r37038 r37807  
    126126
    127127  // define a local projection
     128  InitCoords (&tcoords, "DEC--ARC");
    128129  tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax);
    129130  if (region[0].Dmax < 90) {
    … …  
    132133    tcoords.crval2 = 90.0;
    133134  }
    134   tcoords.crpix1 = 0;
    135   tcoords.crpix2 = 0;
    136135  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
    137   tcoords.pc1_1 = tcoords.pc2_2 = 1.0;
    138   tcoords.pc1_2 = tcoords.pc2_1 = 0.0;
    139   tcoords.Npolyterms = 1;
    140136  strcpy (tcoords.ctype, "DEC--ARC");
    141137
    … …  
    256252    // distinguishing a bad fit from a fit with 1 or 2 bad points
    257253    FIT_MODE = FIT_PM_ONLY;
    258     UpdateObjects (&testcat, 1);
     254    UpdateObjects (&testcat, 1, 0);
    259255
    260256    free_tiny_values(&testcat);
    … …  
    286282
    287283  populate_tiny_values (&catalogOut, DVO_TV_MEASURE);
    288   UpdateObjects (&catalogOut, 1);
     284  UpdateObjects (&catalogOut, 1, 0);
    289285  free_tiny_values(&catalogOut);
    290286

  • trunk/Ohana/src/relastro/src/initialize.c

    r37038 r37807  
    99
    1010  if (DCR_BLUE_COLOR_POS)  fprintf (stderr, "DCR_BLUE_COLOR_POS:  %s - %s\n", DCR_BLUE_COLOR_POS, DCR_BLUE_COLOR_NEG);
    11   if (DCR_BLUE_COLOR_POS)  fprintf (stderr, "DCR_BLUE_COLOR_POS:  %s - %s\n", DCR_BLUE_COLOR_POS, DCR_BLUE_COLOR_NEG);
     11  if (DCR_RED_COLOR_POS)   fprintf (stderr, "DCR_RED_COLOR_POS:   %s - %s\n", DCR_RED_COLOR_POS, DCR_RED_COLOR_NEG);
    1212
    1313  if (PHOTCODE_KEEP_LIST)  fprintf (stderr, "PHOTCODE_KEEP_LIST:  %s\n", PHOTCODE_KEEP_LIST);
    … …  
    6363  initstats (STATMODE);
    6464
     65  if (USE_ICRF_CORRECT) {
     66    if (!USE_ICRF_LOCAL && !USE_ICRF_SHFIT) {
     67      fprintf (stderr, "no ICRF correction method chosen\n");
     68      exit (2);
     69    }
     70  }
     71  if (USE_ICRF_CORRECT) ICRFinit();
     72
    6573  /* XXX drop irrelevant entries */
    6674  if (SHOW_PARAMS) {

  • trunk/Ohana/src/relastro/src/launch_region_hosts.c

    r37038 r37807  
    66  int i;
    77
    8   // clear the I/O files
    9   for (i = 0; i < regionHosts->Nhosts; i++) {
    10     char *syncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "meanmags.sync");
    11     if (truncate (syncfile, 0)) fprintf (stderr, "trouble clearing syncfile %s\n", syncfile);
    12     free (syncfile);
     8  // do not remote the sync and fits files if we do a manual run -- user must clear if needed
     9  if (!PARALLEL_REGIONS_MANUAL) {
     10    // clear the I/O files
     11    for (i = 0; i < regionHosts->Nhosts; i++) {
     12      char *meansync = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "meanpos.sync");
     13      if (truncate (meansync, 0)) fprintf (stderr, "trouble clearing meansync %s\n", meansync);
     14      free (meansync);
    1315
    14     char *fitsfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "meanmags.fits");
    15     if (truncate (fitsfile, 0)) fprintf (stderr, "trouble clearing fitsfile %s\n", fitsfile);
    16     free (fitsfile);
     16      char *meanfits = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "meanpos.fits");
     17      if (truncate (meanfits, 0)) fprintf (stderr, "trouble clearing meanfits %s\n", meanfits);
     18      free (meanfits);
    1719
    18     char *imsyncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "imagemags.sync");
    19     if (truncate (imsyncfile, 0)) fprintf (stderr, "trouble clearing imsyncfile %s\n", imsyncfile);
    20     free (imsyncfile);
     20      char *meassync = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "measpos.sync");
     21      if (truncate (meassync, 0)) fprintf (stderr, "trouble clearing meassync %s\n", meassync);
     22      free (meassync);
    2123
    22     char *imfitsfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "imagemags.fits");
    23     if (truncate (imfitsfile, 0)) fprintf (stderr, "trouble clearing imfitsfile %s\n", imfitsfile);
    24     free (imfitsfile);
     24      char *measfits = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "measpos.fits");
     25      if (truncate (measfits, 0)) fprintf (stderr, "trouble clearing measfits %s\n", measfits);
     26      free (measfits);
    2527
    26     char *loopsyncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "loop.sync");
    27     if (truncate (loopsyncfile, 0)) fprintf (stderr, "trouble clearing loopsyncfile %s\n", loopsyncfile);
    28     free (loopsyncfile);
     28      char *icrfsync = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "icrfobj.sync");
     29      if (truncate (icrfsync, 0)) fprintf (stderr, "trouble clearing icrfsync %s\n", icrfsync);
     30      free (icrfsync);
     31
     32      char *icrffits = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "icrfobj.fits");
     33      if (truncate (icrffits, 0)) fprintf (stderr, "trouble clearing icrffits %s\n", icrffits);
     34      free (icrffits);
     35
     36      char *imsyncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "imagepos.sync");
     37      if (truncate (imsyncfile, 0)) fprintf (stderr, "trouble clearing imsyncfile %s\n", imsyncfile);
     38      free (imsyncfile);
     39
     40      char *imfitsfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "imagepos.fits");
     41      if (truncate (imfitsfile, 0)) fprintf (stderr, "trouble clearing imfitsfile %s\n", imfitsfile);
     42      free (imfitsfile);
     43
     44      char *loopsyncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "loop.sync");
     45      if (truncate (loopsyncfile, 0)) fprintf (stderr, "trouble clearing loopsyncfile %s\n", loopsyncfile);
     46      free (loopsyncfile);
     47    }
     48
     49    char *framesync = make_filename (CATDIR, "master", 0, "frame.corr.sync");
     50    if (truncate (framesync, 0)) fprintf (stderr, "trouble clearing framesync %s\n", framesync);
     51    free (framesync);
     52 
     53    char *framefits = make_filename (CATDIR, "master", 0, "frame.corr.fits");
     54    if (truncate (framefits, 0)) fprintf (stderr, "trouble clearing framefits %s\n", framefits);
     55    free (framefits);
    2956  }
    3057
    … …  
    3764    char filename[1024];
    3865    snprintf (filename, 1024, "%s/Image.%d.fits", CATDIR, host->hostID);
     66    if (unlink (filename)) fprintf (stderr, "trouble clearing image %s\n", filename);
    3967
    4068    // write the image subset for this host
    4169    ImageTableSave (filename, host->image, host->Nimage);
     70
     71    if (host->astromTable) {
     72      char mapname[1024];
     73      snprintf (mapname, 1024, "%s/AstroMap.%d.fits", CATDIR, host->hostID);
     74
     75      // write the image subset for this host
     76      AstromOffsetMapSave (host->astromTable, mapname);
     77    }
    4278
    4379    char command[1024];
    … …  
    92128    if (PhotFlagPoor)       strextend (command, "+photflagpoor %d", PhotFlagPoor);
    93129
     130    if (CHIPMAP)            strextend (command, "-chipmap %d", CHIPMAP);
     131    if (RESET_IMAGES)       strextend (command, "-reset-images");
     132
    94133    if (MinBadQF > 0.0)        strextend (command, "-min-bad-psfqf %f", MinBadQF);
    95134    if (MaxMeanOffset != 10.0) strextend (command, "-max-mean-offset  %f", MaxMeanOffset);
    … …  
    99138    strextend (command, "-D RELASTRO_DPOS_MAX %f", DPOS_MAX);
    100139    strextend (command, "-D ADDSTAR_RADIUS %f", ADDSTAR_RADIUS);
     140
     141    strextend (command, "-D USE_GALAXY_MODEL %d", USE_GALAXY_MODEL);
     142
     143    strextend (command, "-D USE_ICRF_CORRECT %d", USE_ICRF_CORRECT);
     144    strextend (command, "-D USE_ICRF_LOCAL %d",   USE_ICRF_LOCAL);
     145    strextend (command, "-D USE_ICRF_SHFIT %d",   USE_ICRF_SHFIT);
     146    strextend (command, "-D USE_ICRF_POLE %d",    USE_ICRF_POLE);
    101147
    102148    if (USE_FIXED_PIXCOORDS) strextend (command, "-D USE_FIXED_PIXCOORDS 1");
    … …  
    133179  }
    134180
     181  if (USE_ICRF_CORRECT) return TRUE;
     182
     183  int status = HarvestRegionHosts (regionHosts);
     184  return status;
     185}
     186
     187int HarvestRegionHosts (RegionHostTable *regionHosts) {
     188
     189  int status = TRUE;
     190
    135191  if (PARALLEL_REGIONS_MANUAL) {
    136192    fprintf (stderr, "run the relastro_client commands above.  when these are done, hit return\n");
    137193    getchar();
    138194  } else {
    139     RegionHostTableWaitJobsGetIO (regionHosts, __FILE__, __LINE__, VERBOSE);
     195    status = RegionHostTableWaitJobsGetIO (regionHosts, __FILE__, __LINE__, VERBOSE);
    140196  }
    141  
    142   int status = TRUE;
    143   for (i = 0; i < regionHosts->Nhosts; i++) {
    144     status = status && (regionHosts->hosts[i].status == 0);
    145   }
    146 
    147197  return status;
    148198}

  • trunk/Ohana/src/relastro/src/liststats.c

    r34576 r37807  
    252252    X2 += M / dM;
    253253    dS += M;
    254     if (XVERB) fprintf (stderr, "chisq contrib : %f : %f : %f : %f : %f\n", Mo, 3600.0*(value[i] - Mo), 3600.0*dvalue[i], (value[i] - Mo) / dvalue[i], X2);
     254    if (XVERB) fprintf (stderr, "chisq contrib : %f : %f : %f : %f : %f\n", Mo, (value[i] - Mo), dvalue[i], (value[i] - Mo) / dvalue[i], X2);
    255255  }
    256256  X2 = X2 / (Nm - 1);

  • trunk/Ohana/src/relastro/src/load_catalogs.c

    r37038 r37807  
    66Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath, char *syncfile) {
    77
    8   int i, j, Nstar;
     8  int i, j;
    99  // int k, m;
    1010  Catalog *catalog, *pcatalog, tcatalog;
    … …  
    8888
    8989  bcatalog_show_skips();
    90 
    91   Nstar = 0;
     90  fprintf (stderr, "included %d ICRF QSOs\n", ICRFmax());
     91
     92  int Nstar = 0;
     93  int Nmeas = 0;
    9294  for (i = 0; i < Ncat; i++) {
    9395    Nstar += catalog[i].Naverage;
     96    Nmeas += catalog[i].Nmeasure;
    9497    if ((sizeof(IDX_T) == 8) && (catalog[i].Nmeasure > 0xffffffff)) {
    9598      fprintf (stderr, "ERROR: using small-sized IDX_T on data with more than 2G detections per table will cause errors\n");
    … …  
    101104    fprintf (stderr, "warning: insufficient stars %d\n", Nstar);
    102105  }
     106  fprintf (stderr, "using %d stars, %d measurements to calibrate\n", Nstar, Nmeas);
    103107
    104108  // if we are running with parallel_images but not a parallel database, we need to
    … …  
    170174    if (FlagOutlier)   strextend (command, "-clip %d", CLIP_THRESH);
    171175    if (ExcludeBogus)  strextend (command, "-exclude-bogus %f", ExcludeBogusRadius);
     176
     177    if (USE_ICRF_CORRECT) strextend (command, "-D USE_ICRF_CORRECT %d", USE_ICRF_CORRECT);
     178    if (USE_GALAXY_MODEL) strextend (command, "-D USE_GALAXY_MODEL %d", USE_GALAXY_MODEL);
    172179   
    173180    if (DCR_BLUE_COLOR_POS && DCR_BLUE_COLOR_NEG) {

  • trunk/Ohana/src/relastro/src/load_images.c

    r36630 r37807  
    11# include "relastro.h"
     2
     3static AstromOffsetTable *table = NULL;
     4
     5/* AstromOffsetMap is the map of Astrometry Offsets for non-polynomial astrometric
     6 * corrections for each chip. 
     7 */
    28
    39int load_images (FITS_DB *db, SkyList *skylist, int UseFullOverlap) {
    … …  
    1723  MARKTIME("  convert image table: %f sec\n", dtime);
    1824
     25  // assign image->parent and image->coords.mosaic
     26  BuildChipMatch (image, Nimage);
     27  MARKTIME("build chip match: %f sec\n", dtime);
     28
     29  char mapfile[DVO_MAX_PATH];
     30  snprintf (mapfile, DVO_MAX_PATH, "%s/AstroMap.fits", CATDIR);
     31  table = AstromOffsetMapLoad (mapfile, VERBOSE);
     32
     33  // assign images.coords.offsetMap -> table->map[i]
     34  if (table) {
     35    AstromOffsetTableMatchChips (image, Nimage, table);
     36  } else {
     37    table = AstromOffsetTableInit ();
     38  }
     39
    1940  // select the images which overlap the selected sky regions
    2041  subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset, UseFullOverlap);
    2142  MARKTIME("  select images: %f sec\n", dtime);
     43
     44  if (Nsubset == Nimage) {
     45    free (subset);
     46    subset = image;
     47  }
    2248
    2349  initImages (subset, LineNumber, Nsubset);
    … …  
    5177  return TRUE;
    5278}
     79
     80int save_astrom_table () {
     81
     82  char mapfile[DVO_MAX_PATH];
     83  snprintf (mapfile, DVO_MAX_PATH, "%s/AstroMap.fits", CATDIR);
     84  AstromOffsetMapSave (table, mapfile);
     85
     86  return TRUE;
     87}
     88
     89AstromOffsetTable *get_astrom_table () {
     90  return table;
     91}
     92
     93void put_astrom_table (AstromOffsetTable *myTable) {
     94  table = myTable;
     95}
     96

  • trunk/Ohana/src/relastro/src/markObjects.c

    r37038 r37807  
    1414      int m = catalog[i].average[j].measureOffset;
    1515      for (n = 0; n < catalog[i].average[j].Nmeasure; n++) {
    16         if (!catalog[i].measureT[m+n].myDet) continue;
    17         nOwn ++;
     16        // this is a detection I own (it is on one of my images)
     17        if (catalog[i].measureT[m+n].myDet) {
     18          nOwn ++;
     19          continue;
     20        }
     21        // 2MASS detections are counted as 'owned' (I completely control this object is nOwn_t == Nmeasure; no sharing needed)
     22        if ((catalog[i].measureT[m+n].photcode == 2011) ||
     23            (catalog[i].measureT[m+n].photcode == 2012) ||
     24            (catalog[i].measureT[m+n].photcode == 2013)) {
     25          nOwn ++;
     26          continue;
     27        }
    1828      }
    1929      catalog[i].nOwn_t[j] = nOwn;

  • trunk/Ohana/src/relastro/src/relastroVisual.c

    r37038 r37807  
    2424/** Spawn a kapa window */
    2525static int initWindow(int *kapid) {
     26  if (NTHREADS) {
     27    fprintf (stderr, "-visual and -threads are incompatible, quitting\n");
     28    exit (2);
     29  }
    2630  if (*kapid == -1) {
    2731    *kapid = KapaOpenNamedSocket("kapa", "relastro_plots");
    … …  
    348352
    349353  int i;
    350   for(i = 0; i < numObj; i++) {
     354  for (i = 0; i < numObj; i++) {
    351355    if (raw[i].mask) continue; // XXX
    352356

  • trunk/Ohana/src/relastro/src/relastro_client.c

    r36833 r37807  
    7979      load_images (&db, skylist, FALSE);
    8080     
     81      // I can free the database after I have loaded the images...
     82
    8183      // // load the image subset table from the specified location
    8284      // off_t Nimage;

  • trunk/Ohana/src/relastro/src/relastro_images.c

    r36833 r37807  
    2323  MARKTIME("load images: %f sec\n", dtime);
    2424
    25   initCoords();
    26 
    2725  /* load catalog data from region files : subselect high-quality measurements */
    28   // XXX pass in the image table
    29   // XXX who carries the image grid?
    3026
    3127  // photcodesKeep is used here to allow measurements from the images being calibrated
    … …  
    3329  catalog = load_catalogs (skylist, &Ncatalog, TRUE, 0, NULL, NULL);
    3430  MARKTIME("load catalog data: %f sec\n", dtime);
     31
     32  // find ICRF QSOs for reference downstream (only if USE_ICRF_CORRECT)
     33  select_catalog_ICRF (catalog, Ncatalog);
    3534
    3635  if (photcodesReset) {
    … …  
    6564    case TARGET_SIMPLE:
    6665      for (i = 0; i < NLOOP; i++) {
    67         UpdateObjects (catalog, Ncatalog);
     66        UpdateObjects (catalog, Ncatalog, i);
    6867        UpdateSimple (catalog, Ncatalog);
    6968      }
    … …  
    7170
    7271    case TARGET_CHIPS:
     72      if (RESET_IMAGES) UpdateMeasures (catalog, Ncatalog);
    7373      for (i = 0; i < NLOOP; i++) {
    74         UpdateObjects (catalog, Ncatalog);
    75         UpdateChips (catalog, Ncatalog);
     74        UpdateObjects (catalog, Ncatalog, (i > 0)); // calculate <R>,<D>; if (i > 0), apply Galaxy Motion Model (if desired)
     75        if ((i > 1) || !USE_GALAXY_MODEL) {
     76          // if GALAXY_MODEL is selected, we want to delay the frame correction until we have
     77          // applied the galaxy model a couple of times.
     78          FrameCorrectionSerial (catalog, Ncatalog);
     79        }
     80        UpdateChips (catalog, Ncatalog);   // measure.X,Y -> R,D, fit image.coords
    7681        MARKTIME("update chips: %f sec\n", dtime);
    7782      }
    … …  
    8287    case TARGET_MOSAICS:
    8388      for (i = 0; i < NLOOP; i++) {
    84         UpdateObjects (catalog, Ncatalog);
     89        UpdateObjects (catalog, Ncatalog, i);
    8590        UpdateMosaic (catalog, Ncatalog);
    8691      }
    … …  
    102107  if (PARALLEL) {
    103108    // save the updated image parameters
     109    // need to also save the image map table...
     110    save_astrom_table ();
    104111    dvo_image_update (&db, VERBOSE);
    105112    dvo_image_unlock (&db);
    … …  
    109116  UpdateObjectOffsets (skylist, 0, NULL);
    110117
    111   // iterate over catalogs to make detection coordinates consistant
    112   // FixProblemImages (skylist);
    113 
    114118  if (!PARALLEL) {
    115119    // save the updated image parameters
     120    save_astrom_table ();
    116121    dvo_image_update (&db, VERBOSE);
    117122    dvo_image_unlock (&db);

  • trunk/Ohana/src/relastro/src/relastro_objects.c

    r36833 r37807  
    7777    populate_tiny_values(&catalog, DVO_TV_MEASURE);
    7878
    79     UpdateObjects (&catalog, 1);
     79    UpdateObjects (&catalog, 1, 0);
    8080
    8181    free_tiny_values(&catalog);

  • trunk/Ohana/src/relastro/src/relastro_parallel_images.c

    r36833 r37807  
    1717  INITTIME;
    1818
     19  client_logger_init (CATDIR);
     20
    1921  // load the RegionTable (UserRegion should not be used at this level)
    2022  RegionHostTable *regionHosts = RegionHostTableLoad (CATDIR, REGION_FILE);
    2123  int myHost = regionHosts->index[REGION_HOST_ID];
    2224  RegionHostFindNeighbors (regionHosts, myHost);
     25  client_logger_message ("started parallel images on %s\n", regionHosts->hosts[myHost].hostname);
    2326
    2427  // load the subset images belonging to this host
    … …  
    2932    exit (2);
    3033  }
     34
     35  // assign image->parent and image->coords.mosaic
     36  BuildChipMatch (image, Nimage);
     37
     38  char mapfile[DVO_MAX_PATH];
     39  snprintf (mapfile, DVO_MAX_PATH, "%s/AstroMap.%d.fits", CATDIR, REGION_HOST_ID);
     40  AstromOffsetTable *table = AstromOffsetMapLoad (mapfile, VERBOSE);
     41
     42  // assign images.coords.offsetMap -> table->map[i]
     43  if (table) {
     44    AstromOffsetTableMatchChips (image, Nimage, table);
     45  } else {
     46    table = AstromOffsetTableInit ();
     47  }
     48  put_astrom_table (table);
     49  client_logger_message ("loaded images on %s\n", regionHosts->hosts[myHost].hostname);
    3150
    3251  // once we have read this table, we should remove it for repeat runs
    … …  
    4261  SkyList *skylist = SkyListByBounds (sky, -1, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax);
    4362
    44   initCoords();
    45 
    4663  /* load catalog data from region files (hostID is 0 since we are not a client */
    4764  char *syncfile = make_filename (CATDIR, regionHosts->hosts[myHost].hostname, REGION_HOST_ID, "loadcat.sync");
    4865  catalog = load_catalogs (skylist, &Ncatalog, TRUE, 0, NULL, syncfile);
    49   MARKTIME("-- load catalog data: %f sec\n", dtime);
     66  MARKTIME("-- load catalog data, host %d: %f sec\n", REGION_HOST_ID, dtime);
    5067  free (syncfile);
    5168 
     69  client_logger_message ("loaded catalog data %s\n", regionHosts->hosts[myHost].hostname);
     70
     71  // find ICRF QSOs for reference downstream (only if USE_ICRF_CORRECT)
     72  select_catalog_ICRF (catalog, Ncatalog);
     73
    5274  if (photcodesReset) {
    5375    photcodesKeep  = photcodesReset;
    … …  
    6082  /* match measurements with images, mosaics */
    6183  initImageBins  (catalog, Ncatalog, TRUE);
    62   MARKTIME("-- make image bins: %f sec\n", dtime);
     84  MARKTIME("-- make image bins, host %d: %f sec\n", REGION_HOST_ID, dtime);
    6385
    6486  findImages (catalog, Ncatalog, TRUE);
    65   MARKTIME("-- set up image indexes: %f sec\n", dtime);
     87  MARKTIME("-- set up image indexes, host %d: %f sec\n", REGION_HOST_ID, dtime);
    6688
    6789  // set test points based on the starmap
    6890  createStarMap (catalog, Ncatalog);
     91
     92  client_logger_message ("set up image indexes on %s\n", regionHosts->hosts[myHost].hostname);
    6993
    7094  markObjects (catalog, Ncatalog);
    … …  
    7296  SAVEPLOT = FALSE;
    7397
     98  client_logger_message ("starting the loops: %s\n", regionHosts->hosts[myHost].hostname);
     99
    74100  /* major modes */
    75101  switch (FIT_TARGET) {
    76102    case TARGET_SIMPLE:
    77103      for (i = 0; i < NLOOP; i++) {
    78         UpdateObjects (catalog, Ncatalog);
     104        UpdateObjects (catalog, Ncatalog, i);
    79105        share_mean_pos (catalog, Ncatalog, regionHosts, i);
    80106        slurp_mean_pos (catalog, Ncatalog, regionHosts, i);
    81107        UpdateSimple (catalog, Ncatalog);
    82108        share_image_pos (regionHosts, i);
    83         slurp_image_pos (regionHosts, i);
     109        slurp_image_pos (catalog, Ncatalog, regionHosts, i);
    84110      }
    85111      break;
    86112
    87113    case TARGET_CHIPS:
     114      if (RESET_IMAGES) {
     115        UpdateMeasures (catalog, Ncatalog);
     116        MARKTIME("reset measures, host %d: %f sec\n", REGION_HOST_ID, dtime);
     117        client_logger_message ("reset measures: %s\n", regionHosts->hosts[myHost].hostname);
     118      }
    88119      for (i = 0; i < NLOOP; i++) {
    89         UpdateObjects (catalog, Ncatalog);
     120        UpdateObjects (catalog, Ncatalog, (i > 0));
     121        MARKTIME("update objects loop %d, host %d: %f sec\n", i, REGION_HOST_ID, dtime);
     122        client_logger_message ("update objects loop %d: %s\n", i, regionHosts->hosts[myHost].hostname);
     123        if ((i > 1) || !USE_GALAXY_MODEL) {
     124          // if GALAXY_MODEL is selected, we want to delay the frame correction until we have
     125          // applied the galaxy model a couple of times.
     126          FrameCorrectionParallelSlave (catalog, Ncatalog, regionHosts, i);
     127          MARKTIME("frame correction loop %d, host %d: %f sec\n", i, REGION_HOST_ID, dtime);
     128          client_logger_message ("frame correction loop %d: %s\n", i, regionHosts->hosts[myHost].hostname);
     129        }
    90130        share_mean_pos (catalog, Ncatalog, regionHosts, i);
    91131        slurp_mean_pos (catalog, Ncatalog, regionHosts, i);
     132        client_logger_message ("slurp mags loop %d: %s\n", i, regionHosts->hosts[myHost].hostname);
     133
    92134        UpdateChips (catalog, Ncatalog);
     135        MARKTIME("update chips loop %d, host %d: %f sec\n", i, REGION_HOST_ID, dtime);
     136        client_logger_message ("update chips loop %d: %s\n", i, regionHosts->hosts[myHost].hostname);
     137
    93138        share_image_pos (regionHosts, i);
    94         slurp_image_pos (regionHosts, i);
    95         MARKTIME("update chips: %f sec\n", dtime);
     139        slurp_image_pos (catalog, Ncatalog, regionHosts, i);
     140        client_logger_message ("slurp image pos loop %d: %s\n", i, regionHosts->hosts[myHost].hostname);
     141
     142        share_meas_pos (catalog, Ncatalog, regionHosts, i);
     143        slurp_meas_pos (catalog, Ncatalog, regionHosts, i);
     144        MARKTIME("done exchange loop %d, host %d: %f sec\n", i, REGION_HOST_ID, dtime);
     145        client_logger_message ("slurp meas loop %d: %s\n", i, regionHosts->hosts[myHost].hostname);
    96146      }
    97147      // create summary plots of the process
    … …  
    101151    case TARGET_MOSAICS:
    102152      for (i = 0; i < NLOOP; i++) {
    103         UpdateObjects (catalog, Ncatalog);
     153        UpdateObjects (catalog, Ncatalog, i);
    104154        share_mean_pos (catalog, Ncatalog, regionHosts, i);
    105155        slurp_mean_pos (catalog, Ncatalog, regionHosts, i);
    106156        UpdateMosaic (catalog, Ncatalog);
    107157        share_image_pos (regionHosts, i);
    108         slurp_image_pos (regionHosts, i);
     158        slurp_image_pos (catalog, Ncatalog, regionHosts, i);
    109159      }
    110160      break;
    … …  
    128178
    129179  share_image_pos (regionHosts, -1);
     180  client_logger_message ("share image pos loop %d: %s\n", -1, regionHosts->hosts[myHost].hostname);
    130181
    131182  exit (0);

  • trunk/Ohana/src/relastro/src/relastro_parallel_regions.c

    r36630 r37807  
    4444  if (!launch_region_hosts (regionHosts)) Shutdown ("error launching region hosts");
    4545
     46  // If we are doing FrameCorrection, we have to loop here waiting for ICRF objects from the region
     47  // hosts.  we then measure the frame correction and send back the result
     48  if (!FrameCorrectionParallelMaster (regionHosts)) Shutdown ("error running parallel frame correction");
     49
    4650  // retrieve updated image parameters from the remote hosts (also set Image.mcal)
    47   if (!slurp_image_pos (regionHosts, -1)) Shutdown ("error loading image updates");
     51  if (!slurp_image_pos (NULL, 0, regionHosts, -1)) Shutdown ("error loading image updates");
    4852
    4953  if (!UPDATE) {
    … …  
    5963
    6064  if (PARALLEL) {
    61     // save the updated image parameters
     65    // save the updated image parameters.  if we are NOT PARALLEL, we have to wait
     66    // until after UpdateObjectOffsets, or the image parameters are swapped and/or freed
     67    save_astrom_table ();
    6268    dvo_image_save (&db, VERBOSE);
    6369    dvo_image_unlock (&db);
    … …  
    6571
    6672  // iterate over catalogs to make detection coordinates consistant
    67   UpdateObjectOffsets (skylist, 0, NULL);
     73  if (APPLY_OFFSETS) {
     74    UpdateObjectOffsets (skylist, 0, NULL);
     75  }
    6876
    6977  if (!PARALLEL) {
    7078    // save the changes to the image parameters
     79    save_astrom_table ();
    7180    dvo_image_save (&db, VERBOSE);
    7281    dvo_image_unlock (&db);

  • trunk/Ohana/src/relastro/src/select_images.c

    r36630 r37807  
    2323  int InRange, found;
    2424  double Ri[5], Di[5], Xi[5], Yi[5], dx, dy;
    25   Coords tcoords;
    2625  SkyRegionCoords *skycoords;
    2726 
    … …  
    4645
    4746  // the comparison is made in the catalog local projection. below we set crval1,2
    48   tcoords.crpix1 = tcoords.crpix2 = 0.0;
     47  Coords tcoords;
     48  InitCoords (&tcoords, "DEC--TAN");
    4949  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
    50   tcoords.pc1_1  = tcoords.pc2_2 = 1.0;
    51   tcoords.pc1_2  = tcoords.pc2_1 = 0.0;
    52   strcpy (tcoords.ctype, "DEC--TAN");
    5350
    5451  ALLOCATE (skycoords, SkyRegionCoords, skylist[0].Nregions);
    … …  
    104101
    105102  if (VERBOSE) fprintf (stderr, "finding images\n");
    106   BuildChipMatch (timage, Ntimage);
    107   MARKTIME("build chip match: %f sec\n", dtime);
    108103
    109104  nimage = 0;
    … …  
    148143    }
    149144   
    150     // this adds 1.3 sec for 3M images
    151     if (!FindMosaicForImage (timage, Ntimage, i)) {
    152       fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i);
    153       continue;
    154     }
    155 
    156145    /* define image corners - note the DIS images (mosaic phu) are special */
    157146    if (!strcmp(&timage[i].coords.ctype[4], "-DIS")) {

  • trunk/Ohana/src/relastro/src/share_images_pos.c

    r36630 r37807  
    3434  free (iposfile);
    3535
     36  // save the per-chip residual table
     37  AstromOffsetTable *table = get_astrom_table ();
     38  if (table) {
     39    char mapname[1024];
     40    snprintf (mapname, 1024, "%s/AstroMapUpdate.%d.fits", CATDIR, REGION_HOST_ID);
     41   
     42    // write the image subset for this host
     43    AstromOffsetMapSave (table, mapname);
     44  }
     45
    3646  char *syncfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "imagepos.sync");
    3747  update_sync_file (syncfile, nloop);
    … …  
    4151}
    4252
    43 int slurp_image_pos (RegionHostTable *regionHosts, int nloop) {
     53int slurp_image_pos (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop) {
    4454
    45   off_t Nimage, i;
     55  off_t Nimage, i, j;
    4656  Image *images = getimages (&Nimage, NULL);
    4757
    … …  
    7484      continue;
    7585    }
     86    Coords *moscoords = images[seq].coords.mosaic;
    7687    images[seq].coords     = image_pos[i].coords    ;
    77     images[seq].dXpixSys   = image_pos[i].dXpixSys  ;
    78     images[seq].dYpixSys   = image_pos[i].dYpixSys  ;
    79     images[seq].imageID    = image_pos[i].imageID   ;
    80     images[seq].nFitAstrom = image_pos[i].nFitAstrom;
    81     images[seq].flags      = image_pos[i].flags     ;
     88    images[seq].coords.mosaic = moscoords;
     89
     90    images[seq].dXpixSys     = image_pos[i].dXpixSys  ;
     91    images[seq].dYpixSys     = image_pos[i].dYpixSys  ;
     92    images[seq].refColorBlue = image_pos[i].refColorBlue;
     93    images[seq].refColorRed  = image_pos[i].refColorRed;
     94    images[seq].imageID      = image_pos[i].imageID   ;
     95    images[seq].nFitAstrom   = image_pos[i].nFitAstrom;
     96    images[seq].flags        = image_pos[i].flags     ;
     97  }
     98
     99  // load the astrometry offset maps, if they exist, and apply
     100  for (i = 0; i < regionHosts->Nhosts; i++) {
     101    if (regionHosts->hosts[i].hostID == REGION_HOST_ID) continue;
     102    if (REGION_HOST_ID && !regionHosts->hosts[i].isNeighbor) continue;
     103
     104    char mapname[1024];
     105    snprintf (mapname, 1024, "%s/AstroMapUpdate.%d.fits", CATDIR, regionHosts->hosts[i].hostID);
     106   
     107    AstromOffsetTable *table = AstromOffsetMapLoad (mapname, VERBOSE);
     108 
     109    // apply table entries here to existing images
     110    for (j = 0; j < table->Nmap; j++) {
     111      off_t seq = getImageByID (table->map[j][0].imageID);
     112      // I do not necessarily own all images listed in the table.  skip if not found
     113      if (seq < 0) continue;
     114
     115      AstromOffsetMap *oldMap = images[seq].coords.offsetMap;
     116      AstromOffsetMap *newMap = table->map[j];
     117
     118      if (oldMap) {
     119        AstromOffsetMapSetOrder (oldMap, newMap->Nx, newMap->Ny, &images[seq]);
     120        AstromOffsetMapCopyData (oldMap, newMap);
     121      } else {
     122        lockUpdateChips ();
     123        AstromOffsetTable *FullTable = get_astrom_table ();
     124        AstromOffsetTableNewMap(FullTable, newMap->Nx, newMap->Ny, &images[seq]);  // registers the map with the image
     125        unlockUpdateChips ();
     126        AstromOffsetMapCopyData (images[seq].coords.offsetMap, newMap); // oldMap is now not valid
     127      }
     128    }
     129    AstromOffsetTableFree(table);
     130  }
     131
     132  // apply the modified image parameters to my detections
     133  for (i = 0; FALSE && (Ncatalog > 0) && (i < Nimage_pos); i++) {
     134    off_t seq = getImageByID (image_pos[i].imageID);
     135    if (seq < 0) continue;
     136    updateImageRaw (catalog, Ncatalog, seq);
    82137  }
    83138  free (image_pos);
    … …  
    93148  image_pos->dXpixSys     = image->dXpixSys;
    94149  image_pos->dYpixSys     = image->dYpixSys;
     150  image_pos->refColorBlue = image->refColorBlue;
     151  image_pos->refColorRed  = image->refColorRed;
    95152  image_pos->imageID      = image->imageID;
    96153  image_pos->nFitAstrom   = image->nFitAstrom;

  • trunk/Ohana/src/relastro/src/testparallax.c

    r36482 r37807  
    4949
    5050  /* project coordinates to a plane centered on the object with units of arcsec */
    51   coords.crval1 = 0;
    52   coords.crval2 = 0;
    53   coords.crpix1 = 0;
    54   coords.crpix2 = 0;
     51  InitCoords (&coords, "DEC--SIN");
    5552  coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
    56   coords.pc1_1  = coords.pc2_2 = 1.0;
    57   coords.pc1_2  = coords.pc2_1 = 0.0;
    58   coords.Npolyterms = 1;
    59   strcpy (coords.ctype, "DEC--SIN");
    6053
    6154  // use one point as a local reference
Note: See TracChangeset for help on using the changeset viewer.