IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 8, 2011, 2:44:12 PM (15 years ago)
Author:
eugene
Message:

merge changes from trunk

Location:
branches/eam_branches/ipp-20110906
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20110906

  • branches/eam_branches/ipp-20110906/ppTranslate/src/ppMopsMerge.c

    r32175 r32630  
    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.