IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 15, 2011, 4:50:36 PM (15 years ago)
Author:
bills
Message:

Modify ppMops to use psFitsReadTableAllColumns and psFitsWriteTableAllColumns. These
functions store the columns in psVectors. This avoids the memory explosion that
results by storing the table as an array of metadatas representing rows.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppTranslate/src/ppMopsMerge.c

    r32175 r32406  
    1717#define AIRMASS_TOL 1.0e-3              // Tolerance for airmass matching
    1818
     19#if 0
     20#undef psTrace
     21#define psTrace(facil, level, ...)              \
     22    if (level <= 5) {                           \
     23        fprintf(stderr, __VA_ARGS__);           \
     24    }
     25#endif
     26
    1927// Get distance from detection to centre of image
    2028static float mergeDistance(const ppMopsDetections *detections, // Detections of interest
    … …  
    2230    )
    2331{
    24   float dx = (float) (detections->x->data.F32[index] - detections->naxis1->data.S32[index] / 2.0);
    25   float dy = (float) (detections->y->data.F32[index] - detections->naxis2->data.S32[index] / 2.0);
    26   return PS_SQR(dx) + PS_SQR(dy);
     32    float dx = detections->x->data.F32[index] - detections->naxis1 / 2.0;
     33    float dy = detections->y->data.F32[index] - detections->naxis2 / 2.0;
     34    return PS_SQR(dx) + PS_SQR(dy);
    2735}
    2836
    2937
    30 ppMopsDetections *ppMopsMerge(const psArray *detections)
     38bool ppMopsPurgeDuplicates(const psArray *detections)
    3139{
    3240    PS_ASSERT_ARRAY_NON_NULL(detections, NULL);
    3341
    34     psTrace("ppMops.merge", 1, "Merging detections from %ld inputs\n", detections->n);
    35 
    36     ppMopsDetections *merged = NULL;    // Merged list
    37     int num = 0;                                                         // Number of merged files
    38     psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of interest
     42    int numInputs = detections->n;                // Number of inputs
     43    psTrace("ppMops.merge", 1, "Checking detections from %d inputs\n", numInputs);
     44
     45    long total = 0;                                // Total number of sources
     46    psArray *duplicates = psArrayAlloc(numInputs); // Vector of duplicate bits for each input
     47    psVector *dupNum = psVectorAlloc(numInputs, PS_TYPE_U32); // Number of duplicates for each input
     48    psVectorInit(dupNum, 0);
     49    for (int i = 0; i < numInputs; i++) {
     50        ppMopsDetections *det = detections->data[i]; // Detections from
     51        if (!det || det->num == 0) {
     52            continue;
     53        }
     54        psVector *dupes = duplicates->data[i] = psVectorAlloc(det->num, PS_TYPE_U8);
     55        psVectorInit(dupes, 0);
     56        total += det->num;
     57    }
     58    psTrace("ppMops.merge", 2, "Total detections: %ld\n", total);
     59
     60    psVector *raMerged = psVectorAlloc(total, PS_TYPE_F64); // Merged RAs
     61    psVector *decMerged = psVectorAlloc(total, PS_TYPE_F64); // Merged Decs
     62    psVector *sourceMerged = psVectorAlloc(total, PS_TYPE_U16); // Source image of merged sources
     63    psVector *indexMerged = psVectorAlloc(total, PS_TYPE_U32); // Index of merged sources
     64    long num = 0;                                                   // Number of merged sources
     65
     66    const char *raBoresight = NULL, *decBoresight = NULL; // Boresight coordinates
     67    const char *filter = NULL;                    // Filter name
     68    float airmass = NAN, exptime = NAN;           // Exposure details
     69    double posangle = NAN, alt = NAN, az = NAN; // Telescope pointing
     70    double mjd = NAN;                           // Time of exposure
     71
     72    for (int i = 0; i < numInputs; i++) {
     73        ppMopsDetections *det = detections->data[i]; // Detections of interest
     74        if (!det) {
     75            psTrace("ppMops.merge", 3, "Ignoring NULL input %d\n", i);
     76            continue;
     77        } else if (det->num == 0) {
     78            psTrace("ppMops.merge", 3, "Ignoring empty input %d\n", i);
     79            continue;
     80        }
     81
     82        // Check exposure characteristics
     83        if (num == 0) {
     84            raBoresight = det->raBoresight;
     85            decBoresight = det->decBoresight;
     86            filter = det->filter;
     87            airmass = det->airmass;
     88            exptime = det->exptime;
     89            posangle = det->posangle;
     90            alt = det->alt;
     91            az = det->az;
     92        } else {
     93            if (strcmp(raBoresight, det->raBoresight) != 0) {
     94                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure RA values differ: %s vs %s",
     95                        raBoresight, det->raBoresight);
     96                return false;
     97            }
     98            if (strcmp(decBoresight, det->decBoresight) != 0) {
     99                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure Dec values differ: %s vs %s",
     100                        decBoresight, det->decBoresight);
     101                return false;
     102            }
     103            if (strcmp(filter, det->filter) != 0) {
     104                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure filter values differ: %s vs %s",
     105                        filter, det->filter);
     106                return false;
     107            }
     108            if (fabsf(airmass - det->airmass) > AIRMASS_TOL) {
     109                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure airmass values differ: %f vs %f",
     110                        airmass, det->airmass);
     111                return false;
     112            }
     113            if (fabsf(exptime - det->exptime) > EXPTIME_TOL) {
     114                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure exposure time values differ: %f vs %f",
     115                        exptime, det->exptime);
     116                return false;
     117            }
     118            if (fabs(posangle - det->posangle) > POSANGLE_TOL) {
     119                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure position angle values differ: %f vs %f",
     120                        posangle, det->posangle);
     121                return false;
     122            }
     123            if (fabs(alt - det->alt) > BORESIGHT_TOL) {
     124                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure altitude values differ: %lf vs %lf",
     125                        alt, det->alt);
     126                return false;
     127            }
     128            if (fabs(az - det->az) > BORESIGHT_TOL) {
     129                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure azimuth values differ: %lf vs %lf",
     130                        az, det->az);
     131                return false;
     132            }
     133            if (fabs(mjd - det->mjd) > MJD_TOL) {
     134                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure MJD values differ: %lf vs %lf",
     135                        mjd, det->mjd);
     136                return false;
     137            }
     138        }
     139
     140        psTrace("ppMops.merge", 3, "Accepting %ld detections from input %d\n", det->num, i);
     141        memcpy(&raMerged->data.F64[num], det->ra->data.F64, det->num * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
     142        memcpy(&decMerged->data.F64[num], det->dec->data.F64, det->num * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
     143        for (long j = 0, k = num; j < det->num; j++, k++) {
     144            sourceMerged->data.U16[k] = i;
     145            indexMerged->data.U32[k] = j;
     146        }
     147        num += det->num;
     148    }
     149
     150    psTrace("ppMops.merge", 3, "Generating kd-tree from %ld sources\n", num);
     151    psTree *tree = psTreePlant(2, LEAF_SIZE, PS_TREE_SPHERICAL, raMerged, decMerged); // kd tree
     152    if (!tree) {
     153        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate kd tree");
     154        return false;
     155    }
     156
    39157    for (int i = 0; i < detections->n; i++) {
    40158        ppMopsDetections *det = detections->data[i]; // Detections of interest
    41159        if (!det) {
    42             psTrace("ppMops.merge", 3, "Ignoring NULL input %d\n", i);
    43160            continue;
    44161        } else if (det->num == 0) {
    45             psTrace("ppMops.merge", 3, "Ignoring empty input %d\n", i);
    46             continue;
    47         }
    48         num++;
    49         if (!merged) {
    50             psTrace("ppMops.merge", 3, "Accepting %ld detections from input %d\n", det->num, i);
    51             merged = psMemIncrRefCounter(det);
    52             continue;
    53         }
    54         psTrace("ppMops.merge", 3, "Merging %ld detections from input %d\n", det->num, i);
    55 
    56         // XXX compare exposure properties
    57         if (strcmp(merged->raBoresight, det->raBoresight) != 0) {
    58             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure RA values differ: %s vs %s",
    59                     merged->raBoresight, det->raBoresight);
    60             return NULL;
    61         }
    62         if (strcmp(merged->decBoresight, det->decBoresight) != 0) {
    63             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure Dec values differ: %s vs %s",
    64                     merged->decBoresight, det->decBoresight);
    65             return NULL;
    66         }
    67         if (strcmp(merged->filter, det->filter) != 0) {
    68             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure filter values differ: %s vs %s",
    69                     merged->filter, det->filter);
    70             return NULL;
    71         }
    72 
    73         if (fabsf(merged->airmass - det->airmass) > AIRMASS_TOL) {
    74             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure airmass values differ: %f vs %f",
    75                     merged->airmass, det->airmass);
    76             return NULL;
    77         }
    78         if (fabsf(merged->exptime - det->exptime) > EXPTIME_TOL) {
    79             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure exposure time values differ: %f vs %f",
    80                     merged->exptime, det->exptime);
    81             return NULL;
    82         }
    83         if (fabs(merged->posangle - det->posangle) > POSANGLE_TOL) {
    84             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure position angle values differ: %f vs %f",
    85                     merged->posangle, det->posangle);
    86             return NULL;
    87         }
    88         if (fabs(merged->alt - det->alt) > BORESIGHT_TOL) {
    89             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure altitude values differ: %lf vs %lf",
    90                     merged->alt, det->alt);
    91             return NULL;
    92         }
    93         if (fabs(merged->az - det->az) > BORESIGHT_TOL) {
    94             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure azimuth values differ: %lf vs %lf",
    95                     merged->az, det->az);
    96             return NULL;
    97         }
    98         if (fabs(merged->mjd - det->mjd) > MJD_TOL) {
    99             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure MJD values differ: %lf vs %lf",
    100                     merged->mjd, det->mjd);
    101             return NULL;
    102         }
    103 
    104         merged->seeing += det->seeing;  // Taking average
    105 
    106         psTree *tree = psTreePlant(2, LEAF_SIZE, PS_TREE_SPHERICAL, merged->ra, merged->dec); // kd tree
    107         if (!tree) {
    108             psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate kd tree");
    109             psFree(merged);
    110             return NULL;
    111         }
    112 
     162            continue;
     163        }
     164        psTrace("ppMops.merge", 3, "Checking %ld detections from input %d\n", det->num, i);
     165
     166        psVector *dupes = duplicates->data[i];            // Duplicates list
     167        psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of interest
    113168        for (int j = 0; j < det->num; j++) {
     169            // XXX: added by Bill
     170            if (det->mask->data.U8[j]) {
     171                // we marked this source marked as bad when we read it, mark it as a duplicate
     172                dupes->data.U8[j] = 0xFF;
     173                continue;
     174            }
     175
    114176            coords->data.F64[0] = det->ra->data.F64[j];
    115177            coords->data.F64[1] = det->dec->data.F64[j];
    … …  
    119181                psFree(coords);
    120182                psFree(tree);
    121                 psFree(merged);
    122                 return NULL;
    123             }
    124             if (indices->n == 0) {
    125                 psTrace("ppMops.merge", 9, "No matches for source %d in input %d\n", j, i);
     183                return false;
     184            }
     185            psTrace("ppMops.merge", 5, "%ld matches for source %d from input %d\n", indices->n, j, i);
     186            psAssert(indices->n > 0, "Expect at least one match for source %d in input %d", j, i);
     187
     188            if (indices->n == 1 && sourceMerged->data.U16[indices->data.S64[0]] == i) {
     189                // It's myself
    126190                psFree(indices);
    127                 ppMopsDetectionsCopySingle(merged, det, j);
    128191                continue;
    129192            }
    130             psTrace("ppMops.merge", 5, "%ld matches for source %d from input %d\n", indices->n, j, i);
    131193
    132194            // Which one do we keep?
    … …  
    134196            long bestIndex = -1;           // Index with best distance
    135197            for (int k = 0; k < indices->n; k++) {
    136                 long index = indices->data.S64[k]; // Index of point
    137                 float distance = mergeDistance(merged, index); // Distance to centre of image
     198                long mergeIndex = indices->data.S64[k]; // Index of point in merged list
     199                int source = sourceMerged->data.U16[mergeIndex]; // Source image
     200                if (source == i) {
     201                    continue;
     202                }
     203                long index = indexMerged->data.U32[mergeIndex];  // Index in source
     204                psVector *dupes = duplicates->data[source];     // Duplicates list
     205                if (dupes->data.U8[index]) {
     206                    continue;
     207                }
     208
     209                float distance = mergeDistance(detections->data[source], index); // Distance to centre of image
    138210                if (distance < bestDistance) {
    139211                    bestDistance = distance;
    … …  
    142214            }
    143215
    144             float distance = mergeDistance(det, j); // Distance to centre of image
    145             if (distance < bestDistance) {
    146                 psTrace("ppMops.merge", 6, "New source clobbers old sources\n");
    147                 // Blow away existing sources
    148                 for (int k = 0; k < indices->n; k++) {
    149                     long index = indices->data.S64[k]; // Index of point
    150                     merged->mask->data.U8[index] = 0xFF;
     216            if (bestIndex >= 0) {
     217                float distance = mergeDistance(det, j); // Distance to centre of image
     218                if (bestIndex >= 0 && distance < bestDistance) {
     219                    psTrace("ppMops.merge", 6, "New source clobbers old sources\n");
     220                    for (int k = 0; k < indices->n; k++) {
     221                        long mergeIndex = indices->data.S64[k]; // Index of point
     222                        int source = sourceMerged->data.U16[mergeIndex]; // Source image
     223                        if (source == i) {
     224                            continue;
     225                        }
     226                        long index = indexMerged->data.U32[mergeIndex];  // Index in source
     227                        psVector *dupes = duplicates->data[source];      // Duplicates list
     228                        if (!dupes->data.U8[index]) {
     229                            dupes->data.U8[index] = 0xFF;
     230                            dupNum->data.U32[source]++;
     231                        }
     232                    }
     233                } else {
     234                    psTrace("ppMops.merge", 6, "Old sources clobber new source\n");
     235                    dupes->data.U8[j] = 0xFF;
     236                    dupNum->data.U32[i]++;
    151237                }
    152                 ppMopsDetectionsCopySingle(merged, det, j);
    153             } else {
    154                 psTrace("ppMops.merge", 6, "Old sources clobber new source\n");
    155238            }
    156239            psFree(indices);
    157240        }
    158 
    159         psTrace("ppMops.merge", 3, "Done merging input %d, %ld merged sources\n", i, merged->num);
    160 
    161         psFree(tree);
    162         ppMopsDetectionsPurge(merged);
    163     }
    164     psFree(coords);
    165 
    166     if (num == 0) {
    167         //All detections were NULL?!
    168         psTrace("ppMops.merge", 3, "All %ld detections were NULL\n", detections->n);
    169         return NULL;
    170     }
    171     psTrace("ppMops.merge", 2, "%ld sources in merged detections list\n", merged->num);
    172 
    173     merged->seeing /= (float) num;
    174 
    175     return merged;
     241        psFree(coords);
     242    }
     243    psFree(tree);
     244    psFree(raMerged);
     245    psFree(decMerged);
     246    psFree(sourceMerged);
     247    psFree(indexMerged);
     248
     249    // Remove duplicates
     250    for (int i = 0; i < detections->n; i++) {
     251        ppMopsDetections *det = detections->data[i]; // Detections of interest
     252        if (!det) {
     253            continue;
     254        } else if (det->num == 0) {
     255            continue;
     256        }
     257        psTrace("ppMops.merge", 3, "Purging %d duplicates from input %d\n", dupNum->data.U32[i], i);
     258
     259#define VECTOR_PURGE_CASE(TYPE)                                     \
     260        case PS_TYPE_##TYPE:    {                                   \
     261          long j = 0;                                               \
     262          for (long i = 0; i < vector->n; i++) {                    \
     263              if (!dupes->data.U8[i]) {                             \
     264                  if (i == j) {                                     \
     265                      j++;                                          \
     266                      continue;                                     \
     267                  }                                                 \
     268                  vector->data.TYPE[j++] = vector->data.TYPE[i];    \
     269              }                                                     \
     270          }                                                         \
     271          vector->n = j;                                            \
     272        }                                                           \
     273        break
     274
     275        psVector *dupes = duplicates->data[i]; // Duplicates
     276        psArray *table = psListToArray(det->table->list); // Table of data
     277        long newLength = -1;
     278        for (int t = 0; t < table->n; t++) {
     279            psMetadataItem *item = table->data[t]; // Table item
     280            psAssert(item->type == PS_DATA_VECTOR, "Table column is not a vector: %x", item->type);
     281            psVector *vector = item->data.V; // Vector to purge
     282            switch (vector->type.type) {
     283                VECTOR_PURGE_CASE(U8);
     284                VECTOR_PURGE_CASE(U16);
     285                VECTOR_PURGE_CASE(U32);
     286                VECTOR_PURGE_CASE(U64);
     287                VECTOR_PURGE_CASE(S8);
     288                VECTOR_PURGE_CASE(S16);
     289                VECTOR_PURGE_CASE(S32);
     290                VECTOR_PURGE_CASE(S64);
     291                VECTOR_PURGE_CASE(F32);
     292                VECTOR_PURGE_CASE(F64);
     293              default:
     294                psAbort("Unrecognised vector type: %x", vector->type.type);
     295            }
     296            if (newLength == -1) {
     297                newLength = vector->n;
     298            } else if (vector->n != newLength) {
     299                psAbort("Unexpected new length found : %ld expected: %ld",
     300                                                vector->n, newLength);
     301            }
     302               
     303        }
     304        // XXX IS this safe? Perhaps save in numGood?
     305        det->num = newLength;
     306        psFree(table);
     307    }
     308    psFree(dupNum);
     309    psFree(duplicates);
     310
     311    return true;
    176312}
    177313
Note: See TracChangeset for help on using the changeset viewer.