IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 2, 2012, 10:32:09 AM (14 years ago)
Author:
Serge CHASTEL
Message:

MOPS trails fitting integration

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppTranslate/src/ppMopsRead.c

    r32406 r34279  
    88#include "ppMops.h"
    99
     10static psMetadata* sort_by_increasing_idet(psMetadata *fittedTrail, long length, long fillLength);
     11
    1012/*
    1113  ppMopsRead possibly modifies the args->version if the user did not
    1214  set it explicitely.
    13  */
     15*/
    1416psArray *ppMopsRead(ppMopsArguments *args)
    1517{
    16     psTrace("ppMops.read", 1, "Reading input detections\n");
    17 
    18     psArray *inNames = args->input;          // Input names
    19     long num = inNames->n;                   // Number of inputs
    20     psArray *detections = psArrayAlloc(num); // Array of detections, to return
    21     for (int i = 0; i < num; i++) {
    22         const char *name = inNames->data[i];
    23 
    24         psFits *fits = psFitsOpen(name,  "r"); // FITS file
    25         if (!fits) {
    26             psError(PS_ERR_IO, false, "Unable to open input %d", i);
    27             return false;
    28         }
    29 
    30         psMetadata *header = psFitsReadHeader(NULL, fits); // Primary header
    31         if (!header) {
    32             psError(PS_ERR_IO, false, "Unable to read header %d", i);
    33             return false;
    34         }
    35 
    36         psS64 diffSkyfileId = psMetadataLookupS64(NULL, header, "IMAGEID"); // Identifier for image
    37         if (diffSkyfileId == 0) {
    38             psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find identifier for image %d", i);
    39             return false;
    40         }
    41 
    42         if (!psFitsMoveExtName(fits, "SkyChip.psf")) {
    43             psError(PS_ERR_IO, false, "Unable to move to HDU with detections");
    44             return false;
    45         }
    46         int skyChipPsfVersion = ppMopsGetSkyChipPsfVersion(fits);
    47         if (args->version == 0) {
    48           psTrace("ppMops.read", 1, "Changing args->version to %d\n", skyChipPsfVersion);
    49           args->version = (unsigned short) skyChipPsfVersion;
     18  psTrace("ppMops.read", 1, "Reading input detections\n");
     19
     20  psArray *inNames = args->input;          // Input names
     21  long num = inNames->n;                   // Number of inputs
     22  psArray *detections = psArrayAlloc(num); // Array of detections, to return
     23  for (int i = 0; i < num; i++) {
     24    const char *name = inNames->data[i];
     25
     26    psFits *fits = psFitsOpen(name,  "r"); // FITS file
     27    if (!fits) {
     28      psError(PS_ERR_IO, false, "Unable to open input %d", i);
     29      return false;
     30    }
     31
     32    psMetadata *header = psFitsReadHeader(NULL, fits); // Primary header
     33    if (!header) {
     34      psError(PS_ERR_IO, false, "Unable to read header %d", i);
     35      return false;
     36    }
     37
     38    psS64 diffSkyfileId = psMetadataLookupS64(NULL, header, "IMAGEID"); // Identifier for image
     39    if (diffSkyfileId == 0) {
     40      psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find identifier for image %d", i);
     41      return false;
     42    }
     43
     44    if (!psFitsMoveExtName(fits, "SkyChip.psf")) {
     45      psError(PS_ERR_IO, false, "Unable to move to HDU with detections");
     46      return false;
     47    }
     48    int skyChipPsfVersion = ppMopsGetSkyChipPsfVersion(fits);
     49    if (args->version == 0) {
     50      psTrace("ppMops.read", 1, "Changing args->version to %d\n", skyChipPsfVersion);
     51      args->version = (unsigned short) skyChipPsfVersion;
     52    }
     53    if (skyChipPsfVersion == 0) {
     54      // Try to read with the user specified version?
     55      skyChipPsfVersion = args->version;
     56    }
     57    /* Display a warning message if there are version
     58       inconsistencies between the file and the flag (note that
     59       those inconsistencies might be wanted) */
     60    if (skyChipPsfVersion != args->version) {
     61      if (skyChipPsfVersion > args->version) {
     62        psWarning("The FITS data will be downgraded from PS1_DV%d to PS1_DV%d\n",
     63                  skyChipPsfVersion, args->version);
     64      } else { // Necessarily: skyChipPsfVersion > args->version
     65        psWarning("The FITS data will be upgraded from PS1_DV%d to PS1_DV%d (new values set to default 0, NaN...)\n",
     66                  skyChipPsfVersion, args->version);
     67      }
     68    }
     69
     70    long size = psFitsTableSize(fits); // Size of table
     71    if (size <= 0) {
     72      psErrorStackPrint(stderr, "Unable to determine size of table %d", i);
     73      psErrorClear();
     74      psWarning("Ignoring input %d", i);
     75      psFree(header);
     76      psFitsClose(fits);
     77      continue;
     78    }
     79    ppMopsDetections *det = ppMopsDetectionsAlloc(size);
     80    detections->data[i] = det;
     81    det->component = psStringNCopy(name, strrchr(name, '.') - name); // Strip off extension
     82    det->num = size;
     83    det->diffSkyfileId = diffSkyfileId;
     84
     85    psTrace("ppMops.read", 3, "Reading %ld rows from %s\n", size, (const char*)inNames->data[i]);
     86
     87    det->raBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.RA"));
     88    det->decBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.DEC"));
     89    det->filter = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.FILTER"));
     90    det->airmass = psMetadataLookupF32(NULL, header, "AIRMASS");
     91    det->exptime = psMetadataLookupF32(NULL, header, "EXPTIME");
     92    det->posangle = psMetadataLookupF64(NULL, header, "FPA.POSANGLE");
     93    det->alt = psMetadataLookupF64(NULL, header, "FPA.ALT");
     94    det->az = psMetadataLookupF64(NULL, header, "FPA.AZ");
     95    det->mjd = psMetadataLookupF64(NULL, header, "MJD-OBS") + det->exptime / 2.0 / 3600 / 24;
     96
     97    det->seeing = (float) 0.5 * (psMetadataLookupF32(NULL, header, "FWHM_MAJ") +
     98                                 psMetadataLookupF32(NULL, header, "FWHM_MIN"));
     99
     100    det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns
     101    det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows
     102
     103    psFree(header);
     104
     105    psMetadata *table = psFitsReadTableAllColumns(fits); // Table of interest
     106    if (!table) {
     107      psError(PS_ERR_IO, false, "Unable to read table %d", i);
     108      return NULL;
     109    }
     110    //Sort table by idet if necessary
     111    psMetadata* sortedTable = sort_by_increasing_idet(table, det->num, det->num);
     112    if (sortedTable != table) {
     113      psFree(table);
     114      table = sortedTable;
     115    } // else det->table is sorted
     116    det->table = table;
     117    psTrace("ppMops.read", 9, "PSF table:\n%s\n", psMetadataConfigFormat(det->table));
     118    psTrace("ppMops.read", 9, "End of PSF table\n");
     119
     120    //fittedTrail extension (Supported for releases 2 and above)
     121    if (skyChipPsfVersion >= 2) {
     122      if (!psFitsMoveExtName(fits, "SkyChip.xfit")) {
     123        psTrace("ppMops.read", 3, "No fitted trails extension");
     124      } else {
     125        psTrace("ppMops.read", 3, "Fitted trails extension found\n");
     126        det->fittedTrailsSize = psFitsTableSize(fits);
     127        if (det->fittedTrailsSize <= 0) {
     128          psErrorStackPrint(stderr, "Unable to determine size of fitted trails extension table %d", i);
     129          psTrace("ppMops.read", 3, "No entry in fitted trails extension!!!!\n");
     130          psErrorClear();
     131          det->fittedTrails = NULL;
     132        } else {
     133          det->fittedTrails = psFitsReadTableAllColumns(fits); // Table of interest
     134          if (!det->fittedTrails) {
     135            psError(PS_ERR_IO, false, "Unable to read fittedTrails table %d", i);
     136            return NULL;
     137          }
     138          psMetadata* sortedFittedTrails = sort_by_increasing_idet(det->fittedTrails, det->fittedTrailsSize, det->num);
     139          if (sortedFittedTrails != det->fittedTrails) {
     140            psFree(det->fittedTrails);
     141            det->fittedTrails = sortedFittedTrails;
     142          } // else det->fittedTrails is sorted
     143          psTrace("ppMops.read", 9, "Fitted trail:\n%s\n", psMetadataConfigFormat(det->fittedTrails));
     144          psTrace("ppMops.read", 9, "End of Fitted trail\n");
    50145        }
    51         if (skyChipPsfVersion == 0) {
    52           // Try to read with the user specified version?
    53           skyChipPsfVersion = args->version;
    54         }
    55         /* Display a warning message if there are version
    56            inconsistencies between the file and the flag (note that
    57            those inconsistencies might be wanted) */
    58         if (skyChipPsfVersion != args->version) {
    59           if (skyChipPsfVersion > args->version) {
    60             psWarning("The FITS data will be downgraded from PS1_DV%d to PS1_DV%d\n",
    61                       skyChipPsfVersion, args->version);
    62           } else { // Necessarily: skyChipPsfVersion > args->version
    63             psWarning("The FITS data will be upgraded from PS1_DV%d to PS1_DV%d (new values set to default 0, NaN...)\n",
    64                       skyChipPsfVersion, args->version);           
    65           }
    66         }
    67 
    68         long size = psFitsTableSize(fits); // Size of table
    69         if (size <= 0) {
    70             psErrorStackPrint(stderr, "Unable to determine size of table %d", i);
    71             psErrorClear();
    72             psWarning("Ignoring input %d", i);
    73             psFree(header);
    74             psFitsClose(fits);
    75             continue;
    76         }
    77         ppMopsDetections *det = ppMopsDetectionsAlloc(size);
    78         detections->data[i] = det;
    79         det->component = psStringNCopy(name, strrchr(name, '.') - name); // Strip off extension
    80         det->num = size;
    81         det->diffSkyfileId = diffSkyfileId;
    82 
    83         psTrace("ppMops.read", 3, "Reading %ld rows from %s\n", size, (const char*)inNames->data[i]);
    84 
    85         det->raBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.RA"));
    86         det->decBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.DEC"));
    87         det->filter = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.FILTER"));
    88         det->airmass = psMetadataLookupF32(NULL, header, "AIRMASS");
    89         det->exptime = psMetadataLookupF32(NULL, header, "EXPTIME");
    90         det->posangle = psMetadataLookupF64(NULL, header, "FPA.POSANGLE");
    91         det->alt = psMetadataLookupF64(NULL, header, "FPA.ALT");
    92         det->az = psMetadataLookupF64(NULL, header, "FPA.AZ");
    93         det->mjd = psMetadataLookupF64(NULL, header, "MJD-OBS") + det->exptime / 2.0 / 3600 / 24;
    94 
    95         det->seeing = (float) 0.5 * (psMetadataLookupF32(NULL, header, "FWHM_MAJ") +
    96                                      psMetadataLookupF32(NULL, header, "FWHM_MIN"));
    97 
    98         det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns
    99         det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows
    100 
    101         psFree(header);
    102 
    103         psMetadata *table = psFitsReadTableAllColumns(fits); // Table of interest
    104         if (!table) {
    105             psError(PS_ERR_IO, false, "Unable to read table %d", i);
    106             return NULL;
    107         }
    108         det->table = table;
    109         psFitsClose(fits);
    110         if (args->version == 0) {
    111           if (skyChipPsfVersion < 2) {
    112             // XXX: TODO: Do we need to add dummy vectors for the missing columns?
    113            }
    114         }
    115 
    116         psVector *ra = psMetadataLookupVector(NULL, table, "RA_PSF");
    117         psVector *dec = psMetadataLookupVector(NULL, table, "DEC_PSF");
    118 
    119         det->raErr = psVectorAlloc(size, PS_TYPE_F64);
    120         det->decErr = psVectorAlloc(size, PS_TYPE_F64);
    121         det->mask = psVectorAlloc(size, PS_TYPE_U8);
    122 
    123         // convert ra and dec to radians for use in the purge duplicates function
    124         det->ra = (psVector*)psBinaryOp(NULL, ra, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
    125         det->dec = (psVector*)psBinaryOp(NULL, dec, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
    126         det->x = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "X_PSF"));
    127         det->y = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "Y_PSF"));
    128         if (!det->ra || !det->dec || !det->x || !det->y) {
    129             psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find all of RA, Dec, X and Y columns");
    130             return NULL;
    131         }
    132 
    133         // Add our new vectors to the table so that duplicates and masked items may be purged
    134         psMetadataAddVector(table, PS_LIST_HEAD, "DEC_ERR", 0, NULL, det->decErr);
    135         psMetadataAddVector(table, PS_LIST_HEAD, "RA_ERR", 0, NULL, det->raErr);
    136 
    137         psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->num, det->component);
    138 
    139         psVector *mag    = psMetadataLookupVector(NULL, table, "PSF_INST_MAG");
    140         psVector *magErr = psMetadataLookupVector(NULL, table, "PSF_INST_MAG_SIG");
    141         psVector *xErrV  = psMetadataLookupVector(NULL, table, "X_PSF_SIG");
    142         psVector *yErrV  = psMetadataLookupVector(NULL, table, "Y_PSF_SIG");
    143         psVector *scaleV = psMetadataLookupVector(NULL, table, "PLTSCALE");
    144         psVector *angleV = psMetadataLookupVector(NULL, table, "POSANGLE");
    145         psVector *flagsV = psMetadataLookupVector(NULL, table, "FLAGS");
    146 
    147         double plateScale = 0.0;        // Plate scale
    148         long numGood = 0;               // Number of good rows
    149         for (long row = 0; row < size; row++) {
    150 
    151             psU32 flags = flagsV->data.U32[row]; // psFitsTableGetU32(NULL, table, row, "FLAGS");
    152             if (flags & SOURCE_MASK) {
    153                 psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", row, i, flags);
    154                 det->mask->data.U8[row] = 0xFF;
    155                 continue;
    156             }
    157 
    158             // Calculate error in RA, Dec
    159 
    160             double xErr = xErrV->data.F32[row];
    161             double yErr = yErrV->data.F32[row];
    162             double scale = scaleV->data.F32[row];
    163             double angle = angleV->data.F32[row];
    164 
    165             if (!isfinite(det->x->data.F32[row]) || !isfinite(det->y->data.F32[row]) ||
    166                 !isfinite(det->ra->data.F64[row]) || !isfinite(det->dec->data.F64[row]) ||
    167                 !isfinite(mag->data.F32[row]) || !isfinite(magErr->data.F32[row]) ||
    168                 !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) {
    169                 psTrace("ppMops.read", 10,
    170                         "Discarding row %ld from input %d because of non-finite values: "
    171                         "%f %f %lf %lf %f %f %f %f %f %f",
    172                         row, i,
    173                         det->x->data.F32[row], det->y->data.F32[row],
    174                         det->ra->data.F64[row], det->dec->data.F64[row],
    175                         mag->data.F32[row], magErr->data.F32[row],
    176                         xErr, yErr, scale, angle);
    177                 det->mask->data.U8[row] = 0xFF;
    178                 continue;
    179             }
    180 
    181             // XXX Not at all sure I've got the angles around the right way here...
    182             double cosAngle = cos(angle), sinAngle = sin(angle);
    183             double cosAngle2 = PS_SQR(cosAngle), sinAngle2 = PS_SQR(sinAngle);
    184             double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr);
    185             double errScale = scale / 3600.0;
    186             det->raErr->data.F64[row] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);
    187             det->decErr->data.F64[row] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);
    188 
    189             det->mask->data.U8[row] = 0;
    190             plateScale += scale;
    191             numGood++;
    192         }
    193         det->seeing *= ((float) plateScale) / ((float) numGood);
    194 
    195         // Are we using numGood for anything outside of this function?
    196         det->numGood = numGood;
    197 
    198         if (isfinite(args->zp) && numGood > 0) {
    199             psBinaryOp(mag, mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32));
    200         }
    201 
    202         psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)name);
    203     }
    204 
    205     psTrace("ppMops.read", 1, "Done reading input detections\n");
    206 
    207     return detections;
    208 }
     146      }
     147    }
     148
     149    psFitsClose(fits);
     150
     151    if (args->version == 0) {
     152      if (skyChipPsfVersion < 2) {
     153        // XXX: TODO: Do we need to add dummy vectors for the missing columns?
     154      }
     155    }
     156
     157    psVector *ra = psMetadataLookupVector(NULL, table, "RA_PSF");
     158    psVector *dec = psMetadataLookupVector(NULL, table, "DEC_PSF");
     159
     160    det->raErr = psVectorAlloc(size, PS_TYPE_F64);
     161    det->decErr = psVectorAlloc(size, PS_TYPE_F64);
     162    det->mask = psVectorAlloc(size, PS_TYPE_U8);
     163
     164    // convert ra and dec to radians for use in the purge duplicates function
     165    det->ra = (psVector*)psBinaryOp(NULL, ra, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
     166    det->dec = (psVector*)psBinaryOp(NULL, dec, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
     167    det->x = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "X_PSF"));
     168    det->y = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "Y_PSF"));
     169    if (!det->ra || !det->dec || !det->x || !det->y) {
     170      psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find all of RA, Dec, X and Y columns");
     171      return NULL;
     172    }
     173
     174    // Add our new vectors to the table so that duplicates and masked items may be purged
     175    psMetadataAddVector(table, PS_LIST_HEAD, "DEC_ERR", 0, NULL, det->decErr);
     176    psMetadataAddVector(table, PS_LIST_HEAD, "RA_ERR", 0, NULL, det->raErr);
     177
     178    psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->num, det->component);
     179
     180    psVector *mag    = psMetadataLookupVector(NULL, table, "PSF_INST_MAG");
     181    psVector *magErr = psMetadataLookupVector(NULL, table, "PSF_INST_MAG_SIG");
     182    psVector *xErrV  = psMetadataLookupVector(NULL, table, "X_PSF_SIG");
     183    psVector *yErrV  = psMetadataLookupVector(NULL, table, "Y_PSF_SIG");
     184    psVector *scaleV = psMetadataLookupVector(NULL, table, "PLTSCALE");
     185    psVector *angleV = psMetadataLookupVector(NULL, table, "POSANGLE");
     186    psVector *flagsV = psMetadataLookupVector(NULL, table, "FLAGS");
     187
     188    double plateScale = 0.0;        // Plate scale
     189    long numGood = 0;               // Number of good rows
     190    for (long row = 0; row < size; row++) {
     191
     192      psU32 flags = flagsV->data.U32[row]; // psFitsTableGetU32(NULL, table, row, "FLAGS");
     193      if (flags & SOURCE_MASK) {
     194        psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", row, i, flags);
     195        det->mask->data.U8[row] = 0xFF;
     196        continue;
     197      }
     198
     199      // Calculate error in RA, Dec
     200
     201      double xErr = xErrV->data.F32[row];
     202      double yErr = yErrV->data.F32[row];
     203      double scale = scaleV->data.F32[row];
     204      double angle = angleV->data.F32[row];
     205
     206      if (!isfinite(det->x->data.F32[row]) || !isfinite(det->y->data.F32[row]) ||
     207          !isfinite(det->ra->data.F64[row]) || !isfinite(det->dec->data.F64[row]) ||
     208          !isfinite(mag->data.F32[row]) || !isfinite(magErr->data.F32[row]) ||
     209          !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) {
     210        psTrace("ppMops.read", 10,
     211                "Discarding row %ld from input %d because of non-finite values: "
     212                "%f %f %lf %lf %f %f %f %f %f %f",
     213                row, i,
     214                det->x->data.F32[row], det->y->data.F32[row],
     215                det->ra->data.F64[row], det->dec->data.F64[row],
     216                mag->data.F32[row], magErr->data.F32[row],
     217                xErr, yErr, scale, angle);
     218        det->mask->data.U8[row] = 0xFF;
     219        continue;
     220      }
     221
     222      // XXX Not at all sure I've got the angles around the right way here...
     223      double cosAngle = cos(angle), sinAngle = sin(angle);
     224      double cosAngle2 = PS_SQR(cosAngle), sinAngle2 = PS_SQR(sinAngle);
     225      double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr);
     226      double errScale = scale / 3600.0;
     227      det->raErr->data.F64[row] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);
     228      det->decErr->data.F64[row] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);
     229      det->mask->data.U8[row] = 0;
     230      plateScale += scale;
     231      numGood++;
     232    }
     233    det->seeing *= ((float) plateScale) / ((float) numGood);
     234
     235    // Are we using numGood for anything outside of this function?
     236    det->numGood = numGood;
     237
     238    if (isfinite(args->zp) && numGood > 0) {
     239      psBinaryOp(mag, mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32));
     240    }
     241
     242    psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)name);
     243  }
     244
     245  psTrace("ppMops.read", 1, "Done reading input detections\n");
     246
     247  return detections;
     248}
     249
     250static long* sorted_idets_indexes(psVector *idets);
     251static int compare (const void *a, const void *b);
     252static psVector* fill_F32_vector(psVector* in, int size, long* indexes, int indexes_length);
     253
     254// We want the entries in table/fittedTrails to be sorted by increasing idet
     255static psMetadata* sort_by_increasing_idet(psMetadata *aTable, long length, long fillLength) {
     256  if (aTable == NULL) {
     257    return aTable;
     258  }
     259  psVector *idets = psMetadataLookupVector(NULL, aTable, "IPP_IDET");
     260  if (!idets) {
     261    psError(PS_ERR_PROGRAMMING, true, "Failed to find IPP_IDET column");
     262    return NULL;
     263  }
     264  bool is_sorted = true;
     265  int count = 0;
     266  long previous_value = -1;
     267  while (is_sorted && (count<length) ){
     268    is_sorted = (previous_value < idets->data.S64[count]);
     269    previous_value = idets->data.S64[count];
     270    count += 1;
     271  }
     272  if (is_sorted && (length==fillLength)) {
     273    psTrace("ppMops.read", 5, "Values are already sorted by idet and the table is complete\n");
     274    return aTable;
     275  }
     276  //otherwise sort (even if it should not happen) and fill
     277  long* sortedIndexVector = sorted_idets_indexes(idets);
     278  psMetadata* sortedTable = psMetadataAlloc();
     279  char* column_name = "X_EXT";
     280  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL,
     281                       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
     282  column_name = "Y_EXT";
     283  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL,
     284                       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
     285  column_name = "X_EXT_SIG";
     286  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL,
     287                       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
     288  column_name = "Y_EXT_SIG";
     289  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL,
     290                       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
     291  column_name = "EXT_INST_MAG";
     292  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL,
     293                       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
     294  column_name = "EXT_INST_MAG_SIG";
     295  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL,
     296                       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
     297  column_name = "NPARAMS";
     298  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL,
     299                       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
     300  column_name = "EXT_WIDTH_MAJ";
     301  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL,
     302                       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
     303  column_name = "EXT_WIDTH_MIN";
     304  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL,
     305                       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
     306  column_name = "EXT_THETA";
     307  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL,
     308                       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
     309  column_name = "EXT_WIDTH_MAJ_ERR";
     310  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL,
     311                       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
     312  column_name = "EXT_WIDTH_MIN_ERR";
     313  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL,
     314                       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
     315  column_name = "EXT_THETA_ERR";
     316  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL,
     317                       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
     318  return sortedTable;
     319}
     320
     321static psVector* fill_F32_vector(psVector* in, int length, long* indexes, int indexes_length) {
     322  if (in == NULL) {
     323    return NULL;
     324  }
     325  psVector* out = psVectorAlloc(length, in->type.type);
     326  int current_index = 0;
     327  for (int i = 0; i < indexes_length; i++) {
     328    while ( (current_index < indexes[i]) && (current_index < length) ) {
     329      out->data.F32[current_index] = NAN;
     330      current_index += 1;
     331    }
     332    out->data.F32[current_index] = in->data.F32[i];
     333    current_index += 1;
     334  }
     335  while ( current_index < length ) {
     336    out->data.F32[current_index] = NAN;
     337    current_index += 1;
     338  }
     339  return out;
     340}
     341
     342static int compare (const void *a, const void *b) {
     343  return ( *((long*) a)-*((long*) b));
     344}
     345
     346static long* sorted_idets_indexes(psVector *idets) {
     347  long* sortedIndexes = (long*) malloc(idets->n * sizeof(long));
     348  for (int i=0;i<idets->n;i++) {
     349    sortedIndexes[i] = idets->data.S64[i];
     350  }
     351  qsort(sortedIndexes, idets->n, sizeof(long), (void *) compare);
     352  return sortedIndexes;
     353}
Note: See TracChangeset for help on using the changeset viewer.