IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 7, 2010, 10:44:29 AM (16 years ago)
Author:
Paul Price
Message:

Reverse merging to the 'new old version' of ppMops (r28043), since this is what MOPS wants right now. The 'new version' is at /branches/ppTranslate/

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppTranslate/src/ppMopsMerge.c

    r28220 r28623  
    99#include "ppMops.h"
    1010
    11 #define LEAF_SIZE 2                     // Size of leaf
     11#define LEAF_SIZE 4                     // Size of leaf
    1212#define MATCH_RADIUS SEC_TO_RAD(1.0)    // Matching radius
    1313#define MJD_TOL 1.0/3600.0/24.0         // Tolerance for MJD matching
    … …  
    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 
    2719// Get distance from detection to centre of image
    2820static float mergeDistance(const ppMopsDetections *detections, // Detections of interest
    … …  
    3022    )
    3123{
    32     float dx = detections->x->data.F32[index] - detections->naxis1 / 2.0;
    33     float dy = detections->y->data.F32[index] - detections->naxis2 / 2.0;
     24    float dx = detections->x->data.F32[index] - detections->naxis1->data.S32[index] / 2.0;
     25    float dy = detections->y->data.F32[index] - detections->naxis2->data.S32[index] / 2.0;
    3426    return PS_SQR(dx) + PS_SQR(dy);
    3527}
    3628
    3729
    38 bool ppMopsPurgeDuplicates(const psArray *detections)
     30ppMopsDetections *ppMopsMerge(const psArray *detections)
    3931{
    4032    PS_ASSERT_ARRAY_NON_NULL(detections, NULL);
    4133
    42     int numInputs = detections->n;                // Number of inputs
    43     psTrace("ppMops.merge", 1, "Checking detections from %d inputs\n", numInputs);
     34    psTrace("ppMops.merge", 1, "Merging detections from %ld inputs\n", detections->n);
    4435
    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++) {
     36    ppMopsDetections *merged = NULL;    // Merged list
     37    int num = 1;                                                         // Number of merged files
     38    for (int i = 0; i < detections->n; i++) {
    7339        ppMopsDetections *det = detections->data[i]; // Detections of interest
    7440        if (!det) {
    … …  
    7945            continue;
    8046        }
     47        num++;
     48        if (!merged) {
     49            psTrace("ppMops.merge", 3, "Accepting %ld detections from input %d\n", det->num, i);
     50            merged = psMemIncrRefCounter(det);
     51            continue;
     52        }
     53        psTrace("ppMops.merge", 3, "Merging %ld detections from input %d\n", det->num, i);
    8154
    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             }
     55        // XXX compare exposure properties
     56        if (strcmp(merged->raBoresight, det->raBoresight) != 0) {
     57            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure RA values differ: %s vs %s",
     58                    merged->raBoresight, det->raBoresight);
     59            return NULL;
     60        }
     61        if (strcmp(merged->decBoresight, det->decBoresight) != 0) {
     62            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure Dec values differ: %s vs %s",
     63                    merged->decBoresight, det->decBoresight);
     64            return NULL;
     65        }
     66        if (strcmp(merged->filter, det->filter) != 0) {
     67            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure filter values differ: %s vs %s",
     68                    merged->filter, det->filter);
     69            return NULL;
    13870        }
    13971
    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;
     72        if (fabsf(merged->airmass - det->airmass) > AIRMASS_TOL) {
     73            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure airmass values differ: %f vs %f",
     74                    merged->airmass, det->airmass);
     75            return NULL;
    14676        }
    147         num += det->num;
    148     }
     77        if (fabsf(merged->exptime - det->exptime) > EXPTIME_TOL) {
     78            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure exposure time values differ: %f vs %f",
     79                    merged->exptime, det->exptime);
     80            return NULL;
     81        }
     82        if (fabs(merged->posangle - det->posangle) > POSANGLE_TOL) {
     83            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure position angle values differ: %f vs %f",
     84                    merged->posangle, det->posangle);
     85            return NULL;
     86        }
     87        if (fabs(merged->alt - det->alt) > BORESIGHT_TOL) {
     88            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure altitude values differ: %lf vs %lf",
     89                    merged->alt, det->alt);
     90            return NULL;
     91        }
     92        if (fabs(merged->az - det->az) > BORESIGHT_TOL) {
     93            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure azimuth values differ: %lf vs %lf",
     94                    merged->az, det->az);
     95            return NULL;
     96        }
     97        if (fabs(merged->mjd - det->mjd) > MJD_TOL) {
     98            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure MJD values differ: %lf vs %lf",
     99                    merged->mjd, det->mjd);
     100            return NULL;
     101        }
    149102
    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     }
     103        merged->seeing += det->seeing;  // Taking average
    156104
    157     for (int i = 0; i < detections->n; i++) {
    158         ppMopsDetections *det = detections->data[i]; // Detections of interest
    159         if (!det) {
    160             continue;
    161         } else if (det->num == 0) {
    162             continue;
     105        psTree *tree = psTreePlant(2, LEAF_SIZE, PS_TREE_SPHERICAL, merged->ra, merged->dec); // kd tree
     106        if (!tree) {
     107            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate kd tree");
     108            psFree(merged);
     109            return NULL;
    163110        }
    164         psTrace("ppMops.merge", 3, "Checking %ld detections from input %d\n", det->num, i);
    165111
    166         psVector *dupes = duplicates->data[i];            // Duplicates list
    167112        psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of interest
    168113        for (int j = 0; j < det->num; j++) {
    … …  
    174119                psFree(coords);
    175120                psFree(tree);
    176                 return false;
     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);
     126                psFree(indices);
     127                ppMopsDetectionsCopySingle(merged, det, j);
     128                continue;
    177129            }
    178130            psTrace("ppMops.merge", 5, "%ld matches for source %d from input %d\n", indices->n, j, i);
    179             psAssert(indices->n > 0, "Expect at least one match for source %d in input %d", j, i);
    180 
    181             if (indices->n == 1 && sourceMerged->data.U16[indices->data.S64[0]] == i) {
    182                 // It's myself
    183                 psFree(indices);
    184                 continue;
    185             }
    186131
    187132            // Which one do we keep?
    … …  
    189134            long bestIndex = -1;           // Index with best distance
    190135            for (int k = 0; k < indices->n; k++) {
    191                 long mergeIndex = indices->data.S64[k]; // Index of point in merged list
    192                 int source = sourceMerged->data.U16[mergeIndex]; // Source image
    193                 if (source == i) {
    194                     continue;
    195                 }
    196                 long index = indexMerged->data.U32[mergeIndex];  // Index in source
    197                 psVector *dupes = duplicates->data[source];     // Duplicates list
    198                 if (dupes->data.U8[index]) {
    199                     continue;
    200                 }
    201 
    202                 float distance = mergeDistance(detections->data[source], index); // Distance to centre of image
     136                long index = indices->data.S64[k]; // Index of point
     137                float distance = mergeDistance(merged, index); // Distance to centre of image
    203138                if (distance < bestDistance) {
    204139                    bestDistance = distance;
    … …  
    207142            }
    208143
    209             if (bestIndex >= 0) {
    210                 float distance = mergeDistance(det, j); // Distance to centre of image
    211                 if (bestIndex >= 0 && distance < bestDistance) {
    212                     psTrace("ppMops.merge", 6, "New source clobbers old sources\n");
    213                     for (int k = 0; k < indices->n; k++) {
    214                         long mergeIndex = indices->data.S64[k]; // Index of point
    215                         int source = sourceMerged->data.U16[mergeIndex]; // Source image
    216                         if (source == i) {
    217                             continue;
    218                         }
    219                         long index = indexMerged->data.U32[mergeIndex];  // Index in source
    220                         psVector *dupes = duplicates->data[source];      // Duplicates list
    221                         if (!dupes->data.U8[index]) {
    222                             dupes->data.U8[index] = 0xFF;
    223                             dupNum->data.U32[source]++;
    224                         }
    225                     }
    226                 } else {
    227                     psTrace("ppMops.merge", 6, "Old sources clobber new source\n");
    228                     dupes->data.U8[j] = 0xFF;
    229                     dupNum->data.U32[i]++;
     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;
    230151                }
     152                ppMopsDetectionsCopySingle(merged, det, j);
     153            } else {
     154                psTrace("ppMops.merge", 6, "Old sources clobber new source\n");
    231155            }
    232156            psFree(indices);
    233157        }
    234         psFree(coords);
    235     }
    236     psFree(tree);
    237     psFree(raMerged);
    238     psFree(decMerged);
    239     psFree(sourceMerged);
    240     psFree(indexMerged);
    241158
    242     // Remove duplicates
    243     for (int i = 0; i < detections->n; i++) {
    244         ppMopsDetections *det = detections->data[i]; // Detections of interest
    245         if (!det) {
    246             continue;
    247         } else if (det->num == 0) {
    248             continue;
    249         }
    250         psTrace("ppMops.merge", 3, "Purging %d duplicates from input %d\n", dupNum->data.U32[i], i);
     159        psTrace("ppMops.merge", 3, "Done merging input %d, %ld merged sources\n", i, merged->num);
    251160
    252 #define VECTOR_PURGE_CASE(TYPE)                                      \
    253         case PS_TYPE_##TYPE:                                         \
    254           for (int i = 0, j = 0; i < vector->n; i++) {               \
    255               if (!dupes->data.U8[i]) {                              \
    256                   if (i == j) {                                      \
    257                       j++;                                           \
    258                       continue;                                      \
    259                   }                                                  \
    260                   vector->data.TYPE[j] = vector->data.TYPE[i];       \
    261               }                                                      \
    262           }                                                          \
    263           break;
    264 
    265         psVector *dupes = duplicates->data[i]; // Duplicates
    266         psArray *table = psListToArray(det->table->list); // Table of data
    267         for (int j = 0; j < table->n; j++) {
    268             psMetadataItem *item = table->data[j]; // Table item
    269             psAssert(item->type == PS_DATA_VECTOR, "Table column is not a vector: %x", item->type);
    270             psVector *vector = item->data.V; // Vector to purge
    271             switch (vector->type.type) {
    272                 VECTOR_PURGE_CASE(U8);
    273                 VECTOR_PURGE_CASE(U16);
    274                 VECTOR_PURGE_CASE(U32);
    275                 VECTOR_PURGE_CASE(U64);
    276                 VECTOR_PURGE_CASE(S8);
    277                 VECTOR_PURGE_CASE(S16);
    278                 VECTOR_PURGE_CASE(S32);
    279                 VECTOR_PURGE_CASE(S64);
    280                 VECTOR_PURGE_CASE(F32);
    281                 VECTOR_PURGE_CASE(F64);
    282               default:
    283                 psAbort("Unrecognised vector type: %x", vector->type.type);
    284             }
    285         }
     161        psFree(tree);
     162        ppMopsDetectionsPurge(merged);
    286163    }
    287164
    288     return true;
     165    psTrace("ppMops.merge", 2, "%ld sources in merged detections list\n", merged->num);
     166
     167    merged->seeing /= num;
     168
     169    return merged;
    289170}
    290171
Note: See TracChangeset for help on using the changeset viewer.