IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jun 3, 2010, 9:08:48 PM (16 years ago)
Author:
Paul Price
Message:

Finished rework of ppMops. Not yet tested though.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppTranslate/src/ppMopsMerge.c

    r28209 r28212  
    2222    )
    2323{
    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;
     24    float dx = detections->x->data.F32[index] - detections->naxis1 / 2.0;
     25    float dy = detections->y->data.F32[index] - detections->naxis2 / 2.0;
    2626    return PS_SQR(dx) + PS_SQR(dy);
    2727}
    2828
    2929
    30 ppMopsDetections *ppMopsPurgeDuplicates(const psArray *detections)
     30bool ppMopsPurgeDuplicates(const psArray *detections)
    3131{
    3232    PS_ASSERT_ARRAY_NON_NULL(detections, NULL);
    3333
    3434    int numInputs = detections->n;                // Number of inputs
    35     psTrace("ppMops.merge", 1, "Checking detections from %ld inputs\n", numInputs);
    36 
    37     psArray *dupes = psArrayAlloc(numInputs); // Vector of duplicate bits for each input
     35    psTrace("ppMops.merge", 1, "Checking detections from %d inputs\n", numInputs);
     36
     37    long total = 0;                                // Total number of sources
     38    psArray *duplicates = psArrayAlloc(numInputs); // Vector of duplicate bits for each input
     39    psVector *dupNum = psVectorAlloc(numInputs, PS_TYPE_U32); // Number of duplicates for each input
     40    psVectorInit(dupNum, 0);
    3841    for (int i = 0; i < numInputs; i++) {
    3942        ppMopsDetections *det = detections->data[i]; // Detections from
    40 
    41 
    42 
    43 
    44 
    45         dupes->data[i] = psVector
    46 
    47     ppMopsDetections *merged = NULL;    // Merged list
    48     int num = 1;                                                         // Number of merged files
     43        psVector *dupes = duplicates->data[i] = psVectorAlloc(det->num, PS_TYPE_U8);
     44        psVectorInit(dupes, 0);
     45        total += det->num;
     46    }
     47
     48    psVector *raMerged = psVectorAllocEmpty(total, PS_TYPE_F64); // Merged RAs
     49    psVector *decMerged = psVectorAllocEmpty(total, PS_TYPE_F64); // Merged Decs
     50#if 0
     51    psVector *xMerged = psVectorAllocEmpty(total, PS_TYPE_F32);   // Merged x coords
     52    psVector *yMerged = psVectorAllocEmpty(total, PS_TYPE_F32);   // Merged y coords
     53#endif
     54    psVector *sourceMerged = psVectorAllocEmpty(total, PS_TYPE_U16); // Source image of merged sources
     55    psVector *indexMerged = psVectorAllocEmpty(total, PS_TYPE_U32); // Index of merged sources
     56    long num = 0;                                                   // Number of merged sources
     57
     58    const char *raBoresight = NULL, *decBoresight = NULL; // Boresight coordinates
     59    const char *filter = NULL;                    // Filter name
     60    float airmass = NAN, exptime = NAN;           // Exposure details
     61    double posangle = NAN, alt = NAN, az = NAN; // Telescope pointing
     62    double mjd = NAN;                           // Time of exposure
     63
    4964    for (int i = 0; i < detections->n; i++) {
    5065        ppMopsDetections *det = detections->data[i]; // Detections of interest
    … …  
    5671            continue;
    5772        }
    58         num++;
    59         if (!merged) {
    60             psTrace("ppMops.merge", 3, "Accepting %ld detections from input %d\n", det->num, i);
    61             merged = psMemIncrRefCounter(det);
    62             continue;
    63         }
    64         psTrace("ppMops.merge", 3, "Merging %ld detections from input %d\n", det->num, i);
    65 
    66         // XXX compare exposure properties
    67         if (strcmp(merged->raBoresight, det->raBoresight) != 0) {
    68             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure RA values differ: %s vs %s",
    69                     merged->raBoresight, det->raBoresight);
    70             return NULL;
    71         }
    72         if (strcmp(merged->decBoresight, det->decBoresight) != 0) {
    73             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure Dec values differ: %s vs %s",
    74                     merged->decBoresight, det->decBoresight);
    75             return NULL;
    76         }
    77         if (strcmp(merged->filter, det->filter) != 0) {
    78             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure filter values differ: %s vs %s",
    79                     merged->filter, det->filter);
    80             return NULL;
    81         }
    82 
    83         if (fabsf(merged->airmass - det->airmass) > AIRMASS_TOL) {
    84             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure airmass values differ: %f vs %f",
    85                     merged->airmass, det->airmass);
    86             return NULL;
    87         }
    88         if (fabsf(merged->exptime - det->exptime) > EXPTIME_TOL) {
    89             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure exposure time values differ: %f vs %f",
    90                     merged->exptime, det->exptime);
    91             return NULL;
    92         }
    93         if (fabs(merged->posangle - det->posangle) > POSANGLE_TOL) {
    94             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure position angle values differ: %f vs %f",
    95                     merged->posangle, det->posangle);
    96             return NULL;
    97         }
    98         if (fabs(merged->alt - det->alt) > BORESIGHT_TOL) {
    99             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure altitude values differ: %lf vs %lf",
    100                     merged->alt, det->alt);
    101             return NULL;
    102         }
    103         if (fabs(merged->az - det->az) > BORESIGHT_TOL) {
    104             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure azimuth values differ: %lf vs %lf",
    105                     merged->az, det->az);
    106             return NULL;
    107         }
    108         if (fabs(merged->mjd - det->mjd) > MJD_TOL) {
    109             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure MJD values differ: %lf vs %lf",
    110                     merged->mjd, det->mjd);
    111             return NULL;
    112         }
    113 
    114         merged->seeing += det->seeing;  // Taking average
    115 
    116         psTree *tree = psTreePlant(2, LEAF_SIZE, PS_TREE_SPHERICAL, merged->ra, merged->dec); // kd tree
    117         if (!tree) {
    118             psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate kd tree");
    119             psFree(merged);
    120             return NULL;
    121         }
    122 
     73
     74        // Check exposure characteristics
     75        if (num == 0) {
     76            raBoresight = det->raBoresight;
     77            decBoresight = det->decBoresight;
     78            filter = det->filter;
     79            airmass = det->airmass;
     80            exptime = det->exptime;
     81            posangle = det->posangle;
     82            alt = det->alt;
     83            az = det->az;
     84            continue;
     85        } else {
     86            if (strcmp(raBoresight, det->raBoresight) != 0) {
     87                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure RA values differ: %s vs %s",
     88                        raBoresight, det->raBoresight);
     89                return false;
     90            }
     91            if (strcmp(decBoresight, det->decBoresight) != 0) {
     92                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure Dec values differ: %s vs %s",
     93                        decBoresight, det->decBoresight);
     94                return false;
     95            }
     96            if (strcmp(filter, det->filter) != 0) {
     97                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure filter values differ: %s vs %s",
     98                        filter, det->filter);
     99                return false;
     100            }
     101            if (fabsf(airmass - det->airmass) > AIRMASS_TOL) {
     102                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure airmass values differ: %f vs %f",
     103                        airmass, det->airmass);
     104                return false;
     105            }
     106            if (fabsf(exptime - det->exptime) > EXPTIME_TOL) {
     107                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure exposure time values differ: %f vs %f",
     108                        exptime, det->exptime);
     109                return false;
     110            }
     111            if (fabs(posangle - det->posangle) > POSANGLE_TOL) {
     112                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure position angle values differ: %f vs %f",
     113                        posangle, det->posangle);
     114                return false;
     115            }
     116            if (fabs(alt - det->alt) > BORESIGHT_TOL) {
     117                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure altitude values differ: %lf vs %lf",
     118                        alt, det->alt);
     119                return false;
     120            }
     121            if (fabs(az - det->az) > BORESIGHT_TOL) {
     122                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure azimuth values differ: %lf vs %lf",
     123                        az, det->az);
     124                return false;
     125            }
     126            if (fabs(mjd - det->mjd) > MJD_TOL) {
     127                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure MJD values differ: %lf vs %lf",
     128                        mjd, det->mjd);
     129                return false;
     130            }
     131        }
     132
     133        psTrace("ppMops.merge", 3, "Accepting %ld detections from input %d\n", det->num, i);
     134        memcpy(&raMerged->data.F64[num], det->ra, det->num * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
     135        memcpy(&decMerged->data.F64[num], det->dec, det->num * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
     136        for (long j = 0, k = num; j < det->num; j++, k++) {
     137            sourceMerged->data.U16[k] = i;
     138            indexMerged->data.U32[k] = j;
     139        }
     140        num += det->num;
     141    }
     142
     143    psTrace("ppMops.merge", 3, "Generating kd-tree from %ld sources\n", num);
     144    psTree *tree = psTreePlant(2, LEAF_SIZE, PS_TREE_SPHERICAL, raMerged, decMerged); // kd tree
     145    if (!tree) {
     146        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate kd tree");
     147        return false;
     148    }
     149
     150    for (int i = 0; i < detections->n; i++) {
     151        ppMopsDetections *det = detections->data[i]; // Detections of interest
     152        if (!det) {
     153            continue;
     154        } else if (det->num == 0) {
     155            continue;
     156        }
     157        psTrace("ppMops.merge", 3, "Checking %ld detections from input %d\n", det->num, i);
     158
     159        psVector *dupes = duplicates->data[i];            // Duplicates list
    123160        psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of interest
    124161        for (int j = 0; j < det->num; j++) {
    … …  
    130167                psFree(coords);
    131168                psFree(tree);
    132                 psFree(merged);
    133                 return NULL;
    134             }
    135             if (indices->n == 0) {
    136                 psTrace("ppMops.merge", 9, "No matches for source %d in input %d\n", j, i);
    137                 psFree(indices);
    138                 ppMopsDetectionsCopySingle(merged, det, j);
    139                 continue;
     169                return false;
    140170            }
    141171            psTrace("ppMops.merge", 5, "%ld matches for source %d from input %d\n", indices->n, j, i);
     172            psAssert(indices->n > 0, "Expect at least one match for source %d in input %d", j, i);
    142173
    143174            // Which one do we keep?
    … …  
    145176            long bestIndex = -1;           // Index with best distance
    146177            for (int k = 0; k < indices->n; k++) {
    147                 long index = indices->data.S64[k]; // Index of point
    148                 float distance = mergeDistance(merged, index); // Distance to centre of image
     178                long mergeIndex = indices->data.S64[k]; // Index of point in merged list
     179                int source = sourceMerged->data.U16[mergeIndex]; // Source image
     180                if (source == i) {
     181                    continue;
     182                }
     183                long index = indexMerged->data.U32[mergeIndex];  // Index in source
     184                psVector *dupes = duplicates->data[source];     // Duplicates list
     185                if (dupes->data.U8[index]) {
     186                    continue;
     187                }
     188
     189                float distance = mergeDistance(detections->data[source], index); // Distance to centre of image
    149190                if (distance < bestDistance) {
    150191                    bestDistance = distance;
    … …  
    154195
    155196            float distance = mergeDistance(det, j); // Distance to centre of image
    156             if (distance < bestDistance) {
     197            if (!isfinite(bestDistance) || distance < bestDistance) {
    157198                psTrace("ppMops.merge", 6, "New source clobbers old sources\n");
    158                 // Blow away existing sources
    159199                for (int k = 0; k < indices->n; k++) {
    160                     long index = indices->data.S64[k]; // Index of point
    161                     merged->mask->data.U8[index] = 0xFF;
    162                 }
    163                 ppMopsDetectionsCopySingle(merged, det, j);
     200                    long mergeIndex = indices->data.S64[k]; // Index of point
     201                    int source = sourceMerged->data.U16[mergeIndex]; // Source image
     202                    if (source == i) {
     203                        continue;
     204                    }
     205                    long index = indexMerged->data.U32[mergeIndex];  // Index in source
     206                    psVector *dupes = duplicates->data[source];      // Duplicates list
     207                    if (!dupes->data.U8[index]) {
     208                        dupes->data.U8[index] = 0xFF;
     209                        dupNum->data.U32[source]++;
     210                    }
     211                }
    164212            } else {
    165213                psTrace("ppMops.merge", 6, "Old sources clobber new source\n");
    166             }
    167             psFree(indices);
    168         }
    169 
    170         psTrace("ppMops.merge", 3, "Done merging input %d, %ld merged sources\n", i, merged->num);
    171 
    172         psFree(tree);
    173         ppMopsDetectionsPurge(merged);
    174     }
    175 
    176     psTrace("ppMops.merge", 2, "%ld sources in merged detections list\n", merged->num);
    177 
    178     merged->seeing /= num;
    179 
    180     return merged;
     214                dupes->data.U8[j] = 0xFF;
     215                dupNum->data.U32[i]++;
     216            }
     217        }
     218        psFree(coords);
     219    }
     220    psFree(tree);
     221    psFree(raMerged);
     222    psFree(decMerged);
     223    psFree(sourceMerged);
     224    psFree(indexMerged);
     225
     226    // Remove duplicates
     227    for (int i = 0; i < detections->n; i++) {
     228        ppMopsDetections *det = detections->data[i]; // Detections of interest
     229        if (!det) {
     230            continue;
     231        } else if (det->num == 0) {
     232            continue;
     233        }
     234        psTrace("ppMops.merge", 3, "Purging %d duplicates from input %d\n", dupNum->data.U32[i], i);
     235
     236#define VECTOR_PURGE_CASE(TYPE)                                      \
     237        case PS_TYPE_##TYPE:                                         \
     238          for (int i = 0, j = 0; i < vector->n; i++) {               \
     239              if (!dupes->data.U8[i]) {                              \
     240                  if (i == j) {                                      \
     241                      j++;                                           \
     242                      continue;                                      \
     243                  }                                                  \
     244                  vector->data.TYPE[j] = vector->data.TYPE[i];       \
     245              }                                                      \
     246          }                                                          \
     247          break;
     248
     249        psVector *dupes = duplicates->data[i]; // Duplicates
     250        psArray *table = psListToArray(det->table->list); // Table of data
     251        for (int j = 0; j < table->n; j++) {
     252            psVector *vector = table->data[j]; // Vector to purge
     253            switch (vector->type.type) {
     254                VECTOR_PURGE_CASE(U8);
     255                VECTOR_PURGE_CASE(U16);
     256                VECTOR_PURGE_CASE(U32);
     257                VECTOR_PURGE_CASE(U64);
     258                VECTOR_PURGE_CASE(S8);
     259                VECTOR_PURGE_CASE(S16);
     260                VECTOR_PURGE_CASE(S32);
     261                VECTOR_PURGE_CASE(S64);
     262                VECTOR_PURGE_CASE(F32);
     263                VECTOR_PURGE_CASE(F64);
     264              default:
     265                psAbort("Unrecognised vector type: %x", vector->type.type);
     266            }
     267        }
     268    }
     269
     270    return true;
    181271}
    182272
Note: See TracChangeset for help on using the changeset viewer.