IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 37445


Ignore:
Timestamp:
Sep 29, 2014, 8:41:28 AM (12 years ago)
Author:
eugene
Message:

initial testing of the astrom map IO / application is now working in relastro

Location:
branches/eam_branches/ipp-20140904/Ohana/src/libdvo
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20140904/Ohana/src/libdvo/Makefile

    r37386 r37445  
    8888$(SRC)/AstromOffsetMapIO.$(ARCH).o    \
    8989$(SRC)/AstromOffsetMapOps.$(ARCH).o \
     90$(SRC)/AstromOffsetMapUtils.$(ARCH).o \
    9091$(SRC)/dvo_set_skyregion.$(ARCH).o \
    9192$(SRC)/dvo_set_catdir.$(ARCH).o \

  • branches/eam_branches/ipp-20140904/Ohana/src/libdvo/include/libdvo_astro.h

    r37421 r37445  
    5757  float **dXv;
    5858  float **dYv;
    59   unsigned int ID;
     59  unsigned int tableID;
    6060  unsigned int imageID;
    6161} AstromOffsetMap;
    … …  
    6363typedef struct {
    6464  int Nmap;
    65   AstromOffsetMap *map;
     65  AstromOffsetMap **map;
    6666  int *IDtoSeq;
    67   int MaxID;
     67  int MaxImageID;
     68  int MaxTableID;
    6869} AstromOffsetTable;
    6970
    … …  
    164165int AstromOffsetMapFit (AstromOffsetMap *map, float *x, float *y, float *f, int Npts, int xdir);
    165166
     167/* in AstromOffsetMapUtils.c */
     168int AstromOffsetTableNewMap (AstromOffsetTable *table, int order, Image *image);
     169int AstromOffsetTableMatchChips (Image *images, off_t Nimages, AstromOffsetTable *table);
     170AstromOffsetTable *AstromOffsetTableInit();
     171
    166172# endif

  • branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/AstromOffsetMapIO.c

    r37389 r37445  
    11# include "dvo.h"
    2 
    3 /* AstromOffsetMap is the map of Astrometry Offsets for non-polynomial astrometric
    4  * corrections for each chip. 
    5  */
    62
    73AstromOffsetTable *AstromOffsetMapLoad (char *filename, int VERBOSE) {
    … …  
    4945  }
    5046
    51   // convert the disk table format to internal format:
     47  // AstromOffsetMap_Disk_6x6 *map_disk is an external format stored as an array of maps.
     48  // Convert the disk array of maps to then internal format in a rich structure:
    5249  AstromOffsetTable *table = AstromOffsetMapToTable (map_disk, Nmap);
    5350  free (map_disk);
    … …  
    7774  gfits_create_matrix (&header, &matrix);
    7875
     76  // AstromOffsetMap_Disk_6x6 *map_disk is an external format stored as an array of maps.
     77  // Convert the internal format in a rich structure into a disk array of maps.
    7978  AstromOffsetMap_Disk_6x6 *map_disk = AstromOffsetTableToMap (table, &Nmap);
    8079
    … …  
    105104
    106105  table->Nmap = Nmap;
    107   ALLOCATE (table->map, AstromOffsetMap, Nmap);
     106  ALLOCATE (table->map, AstromOffsetMap *, Nmap);
    108107
    109108  // find the max value of imageID
    110   int MaxID = 0;
     109  int MaxTableID = 0;
     110  int MaxImageID = 0;
    111111  for (i = 0; i < Nmap; i++) {
    112     MaxID = MAX(map_disk[i].imageID, MaxID);
     112    MaxTableID = MAX(map_disk[i].tableID, MaxTableID);
     113    MaxImageID = MAX(map_disk[i].imageID, MaxImageID);
    113114  }
    114   MaxID ++; // we want the outer bound, not the last value
     115  table->MaxTableID = MaxTableID;
     116  table->MaxImageID = MaxImageID;
    115117
    116118  // generate the index and init values to -1
    117   ALLOCATE (table->IDtoSeq, int, MaxID);
    118   for (i = 0; i < MaxID; i++) {
     119  ALLOCATE (table->IDtoSeq, int, MaxImageID + 1);
     120  for (i = 0; i <= MaxImageID; i++) {
    119121    table->IDtoSeq[i] = -1;
    120122  }
    121   table->MaxID = MaxID - 1;
    122123
    123124  // assign the ID values
    124125  for (i = 0; i < Nmap; i++) {
    125     int ID = map_disk[i].imageID;
    126     myAssert (table->IDtoSeq[ID] == -1, "oops, duplicate image IDs");
    127     table->IDtoSeq[ID] = i;
     126    int ImageID = map_disk[i].imageID;
     127    myAssert (table->IDtoSeq[ImageID] == -1, "oops, duplicate image IDs");
     128    table->IDtoSeq[ImageID] = i;
    128129  }
    129130
    … …  
    131132  // full 6x6, saving some memory while doing the analysis)
    132133  for (i = 0; i < Nmap; i++) {
    133     table->map[i].Nx      = map_disk[i].Nx;
    134     table->map[i].Ny      = map_disk[i].Ny;
    135     table->map[i].ID      = map_disk[i].ID;
    136     table->map[i].imageID = map_disk[i].imageID;
     134    ALLOCATE (table->map[i], AstromOffsetMap, 1);
     135
     136    table->map[i][0].Nx           = map_disk[i].Nx;
     137    table->map[i][0].Ny           = map_disk[i].Ny;
     138    table->map[i][0].tableID = map_disk[i].tableID;
     139    table->map[i][0].imageID = map_disk[i].imageID;
    137140   
    138141    // STORE THESE VALUES?
    139     table->map[i].dY      = map_disk[i].dX;
    140     table->map[i].dX      = map_disk[i].dY;
     142    table->map[i][0].dY           = map_disk[i].dX;
     143    table->map[i][0].dX           = map_disk[i].dY;
    141144
    142     ALLOCATE (table->map[i].dXv, float *, map_disk[i].Nx);
    143     ALLOCATE (table->map[i].dYv, float *, map_disk[i].Nx);
     145    ALLOCATE (table->map[i][0].dXv, float *, map_disk[i].Nx);
     146    ALLOCATE (table->map[i][0].dYv, float *, map_disk[i].Nx);
    144147
    145148    for (j = 0; j < map_disk[i].Nx; j++) {
    146       ALLOCATE (table->map[i].dXv[j], float, map_disk[i].Ny);
    147       ALLOCATE (table->map[i].dYv[j], float, map_disk[i].Ny);
     149      ALLOCATE (table->map[i][0].dXv[j], float, map_disk[i].Ny);
     150      ALLOCATE (table->map[i][0].dYv[j], float, map_disk[i].Ny);
    148151
    149152      for (k = 0; k < map_disk[i].Ny; k++) {
    150         table->map[i].dXv[j][k] = map_disk[i].dXv[j][k];
    151         table->map[i].dYv[j][k] = map_disk[i].dYv[j][k];
     153        table->map[i][0].dXv[j][k] = map_disk[i].dXv[j][k];
     154        table->map[i][0].dYv[j][k] = map_disk[i].dYv[j][k];
    152155      }
    153156    }
    154157  }
    155 
    156158  return table;
    157159}
    … …  
    167169  // full 6x6, saving some memory while doing the analysis)
    168170  for (i = 0; i < table->Nmap; i++) {
    169     map_disk[i].Nx       = table->map[i].Nx;
    170     map_disk[i].Ny       = table->map[i].Ny;
    171     map_disk[i].ID       = table->map[i].ID;
    172     map_disk[i].imageID  = table->map[i].imageID;
     171    map_disk[i].Nx       = table->map[i][0].Nx;
     172    map_disk[i].Ny       = table->map[i][0].Ny;
     173    map_disk[i].tableID  = table->map[i][0].tableID;
     174    map_disk[i].imageID  = table->map[i][0].imageID;
    173175   
    174     map_disk[i].dX       = table->map[i].dX;
    175     map_disk[i].dY       = table->map[i].dY;
     176    map_disk[i].dX       = table->map[i][0].dX;
     177    map_disk[i].dY       = table->map[i][0].dY;
    176178
    177179    for (j = 0; j < map_disk[i].Nx; j++) {
    178180      for (k = 0; k < map_disk[i].Ny; k++) {
    179         map_disk[i].dXv[j][k] = table->map[i].dXv[j][k];
    180         map_disk[i].dYv[j][k] = table->map[i].dYv[j][k];
     181        map_disk[i].dXv[j][k] = table->map[i][0].dXv[j][k];
     182        map_disk[i].dYv[j][k] = table->map[i][0].dYv[j][k];
    181183      }
    182184    }

  • branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/AstromOffsetMapOps.c

    r37389 r37445  
    1717    float ymo = y * map->dY - 0.5;
    1818    int   ymi = floor(ymo);
    19     ymi = MAX(0,MIN(map->Ny - 1, ymi)); // force range of ymi to be 0,Ny-1
     19    ymi = MAX(0,MIN(map->Ny - 2, ymi)); // force range of ymi to be 0,Ny-2 (Ny must be > 1)
    2020
    2121    float ymf = ymo - ymi;
    … …  
    2929    float xmo = x * map->dX - 0.5;
    3030    int   xmi = floor(xmo);
    31     xmi = MAX(0,MIN(map->Nx - 1, xmi)); // force range of ymi to be 0,Ny-1
     31    xmi = MAX(0,MIN(map->Nx - 2, xmi)); // force range of ymi to be 0,Nx-2  (Nx must be > 1)
    3232
    3333    float xmf = xmo - xmi;
    … …  
    4141  float xmo = x * map->dX - 0.5;
    4242  int   xmi = floor(xmo);
    43   xmi = MAX(0,MIN(map->Nx - 1, xmi)); // force range of ymi to be 0,Ny-1
     43  xmi = MAX(0,MIN(map->Nx - 2, xmi)); // force range of ymi to be 0,Nx-1
    4444  float xmf = xmo - xmi;
    4545
    4646  float ymo = y * map->dY - 0.5;
    4747  int   ymi = floor(ymo);
    48   ymi = MAX(0,MIN(map->Ny - 1, ymi)); // force range of ymi to be 0,Ny-1
     48  ymi = MAX(0,MIN(map->Ny - 2, ymi)); // force range of ymi to be 0,Ny-2
    4949  float ymf = ymo - ymi;
    5050
    … …  
    113113  for (i = 0; i < Nx*Ny; i++) {
    114114    ALLOCATE (A[i], float, Nx*Ny);
     115    memset (A[i], 0, sizeof(float)*Nx*Ny);
    115116    ALLOCATE (B[i], float, 1);
     117    memset (B[i], 0, sizeof(float));
    116118  }   
    117119
    … …  
    120122  // for (int n = 1; n < Nx - 1; n++) {
    121123  // for (int m = 1; m < Ny - 1; m++) {
     124
     125  if (0) {
     126    FILE *fd = fopen ("stats.dump.txt", "w");
     127    for (i = 0; i < Npts; i++) {
     128      fprintf (fd, "%d %f %f %f\n", i, x[i], y[i], f[i]);
     129    }
     130    fclose(fd);
     131  }
    122132
    123133  // float Total = 0.0;
    … …  
    267277          int J = (ix + jx) + Nx * (iy + jy);
    268278          A[J][I] = sA[jx][jy];
    269           // fprintf (stderr, "A %d %d (%d %d : %d %d): %f\n", I, J, n, m, n + jn, m + jm, sA[jn][jm]);
     279          if (abs(A[J][I]) > 1000) {
     280            fprintf (stderr, "A %d %d (%d %d : %d %d): %f\n", I, J, ix, iy, ix + jx, iy + jy, sA[jx][jy]);
     281          }
    270282          // Sum += sA[jn][jm];
    271283        }
    … …  
    283295  for (i = 0; i < Nx*Ny; i++) {
    284296    Empty[i] = 0;
    285     if (A[i][i] == 0.0) {
     297    if (fabs(A[i][i]) < 1.e-2) {
    286298      Empty[i] = 1;
    287299      for (j = 0; j < Nx*Ny; j++) {
    … …  
    292304      B[i][0] = 0.0;
    293305    }
     306  }
     307
     308  if (0) {
     309    FILE *fd = fopen ("matrix.dat", "w");
     310    for (i = 0; i < Nx*Ny; i++) {
     311      for (j = 0; j < Nx*Ny; j++) {
     312        fprintf (fd, "%10.4f ", A[i][j]);
     313      }
     314      fprintf (fd, " : %10.4f\n", B[i][0]);
     315    }
     316    fclose (fd);
    294317  }
    295318

  • branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/astrom_maps.c

    r37438 r37445  
    2424}
    2525
    26 int
     26int AstromOffsetTableNewMap (AstromOffsetTable *table, int order, Image *image) {
     27
     28  int Nmap = table->Nmap;
     29
     30  table->Nmap++;
     31  REALLOCATE (table->map, AstromOffsetMap, table->Nmap);
     32
     33  int Nx = order;
     34  int Ny = order;
     35
     36  if (image->imageID > table->MaxID) {
     37    int oldMaxID = table->MaxID;
     38    table->MaxID = image->imageID;
     39    REALLOCATE (table->IDtoSeq, int, table->MaxID + 1);
     40    for (i = oldMaxID + 1; i < table->MaxID + 1; i++) {
     41      table->IDtoSeq[i] = -1;
     42    }
     43  }
     44  myAssert (table->IDtoSeq[image->imageID] == -1, "table IDtoSeq not initiazed or image collision");
     45  table->IDtoSeq[image->imageID] = Nmap;
     46 
     47  table->map[Nmap].Nx      = Nx;
     48  table->map[Nmap].Ny      = Ny;
     49  table->map[Nmap].ID      = table->maxID; table->maxID ++;
     50  table->map[Nmap].imageID = image->imageID;
     51 
     52  table->map[Nmap].dX = Nx / image->Nx;
     53  table->map[Nmap].dY = Ny / image->Ny;
     54
     55  ALLOCATE (table->map[Nmap].dXv, float *, Nx);
     56  ALLOCATE (table->map[Nmap].dYv, float *, Nx);
     57
     58  for (j = 0; Nx; j++) {
     59    ALLOCATE (table->map[i].dXv[j], float, Ny);
     60    ALLOCATE (table->map[i].dYv[j], float, Ny);
     61
     62    for (k = 0; k < Ny; k++) {
     63      table->map[i].dXv[j][k] = 0.0;
     64      table->map[i].dYv[j][k] = 0.0;
     65    }
     66  }
     67  image[0].coords.offsetMap = &table->map[Nmap];
     68  return TRUE;   
     69}
     70

  • branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/mosaic_astrom.c

    r37377 r37445  
    207207# endif
    208208
     209// XXX Note that this links (images[i].coords.mosaic, images[i].parent) will break
     210// if we reallocate the Image array in the middle of program
    209211int BuildChipMatch (Image *images, off_t Nimages) {
    210212
Note: See TracChangeset for help on using the changeset viewer.