Changeset 25181 for branches/pap_mops/ppMops/src/ppMopsMerge.c

Timestamp:

Aug 24, 2009, 6:00:00 PM (17 years ago)

Author:

Paul Price

Message:

Wrote functions to merge detections from multiple skycells, and to write detections out. Compiles, not yet tested.

File:

• branches/pap_mops/ppMops/src/ppMopsMerge.c (modified) (1 diff)

branches/pap_mops/ppMops/src/ppMopsMerge.c

@@ -5 +5 @@
 #include <stdio.h>
 #include <pslib.h>
+#include <string.h>
 
 #include "ppMops.h"
 
+#define LEAF_SIZE 4                     // Size of leaf
+#define MATCH_RADIUS 1.0/3600           // Matching radius
+
+// Get distance from detection to centre of image
+static float mergeDistance(const ppMopsDetections *detections, // Detections of interest
+                           long index                          // Index for source of interest
+    )
+{
+    float dx = detections->x->data.F32[index] - detections->naxis1->data.S32[index] / 2.0;
+    float dy = detections->y->data.F32[index] - detections->naxis2->data.S32[index] / 2.0;
+    return PS_SQR(dx) + PS_SQR(dy);
+}
+
+
 ppMopsDetections *ppMopsMerge(const psArray *detections)
 {
+    PS_ASSERT_ARRAY_NON_NULL(detections, NULL);
 
-    return NULL;
+    ppMopsDetections *merged = psMemIncrRefCounter(detections->data[0]); // Merged list
+    int num = 1;                                                         // Number of merged files
+    for (int i = 1; i < detections->n; i++) {
+        ppMopsDetections *det = detections->data[i]; // Detections of interest
+        if (!det) {
+            continue;
+        }
+        num++;
+
+        // XXX compare exposure properties
+        if (strcmp(merged->raBoresight, det->raBoresight) != 0) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure RA values differ: %s vs %s",
+                    merged->raBoresight, det->raBoresight);
+            return NULL;
+        }
+        if (strcmp(merged->decBoresight, det->decBoresight) != 0) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure Dec values differ: %s vs %s",
+                    merged->decBoresight, det->decBoresight);
+            return NULL;
+        }
+        if (strcmp(merged->filter, det->filter) != 0) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure filter values differ: %s vs %s",
+                    merged->filter, det->filter);
+            return NULL;
+        }
+
+        if (merged->airmass != det->airmass) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure airmass values differ: %f vs %f",
+                    merged->airmass, det->airmass);
+            return NULL;
+        }
+        if (merged->exptime != det->exptime) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure exposure time values differ: %f vs %f",
+                    merged->exptime, det->exptime);
+            return NULL;
+        }
+        if (merged->posangle != det->posangle) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure position angle values differ: %f vs %f",
+                    merged->posangle, det->posangle);
+            return NULL;
+        }
+        if (merged->alt != det->alt) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure altitude values differ: %lf vs %lf",
+                    merged->alt, det->alt);
+            return NULL;
+        }
+        if (merged->az != det->az) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure azimuth values differ: %lf vs %lf",
+                    merged->az, det->az);
+            return NULL;
+        }
+        if (merged->mjd != det->mjd) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure MJD values differ: %lf vs %lf",
+                    merged->mjd, det->mjd);
+            return NULL;
+        }
+
+        merged->seeing += det->seeing;  // Taking average
+
+        psTree *tree = psTreePlant(2, LEAF_SIZE, PS_TREE_SPHERICAL, merged->ra, merged->dec); // kd tree
+        if (!tree) {
+            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate kd tree");
+            psFree(merged);
+            return NULL;
+        }
+
+        psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of interest
+        for (int j = 0; j < det->num; j++) {
+            coords->data.F64[0] = det->ra->data.F64[j];
+            coords->data.F64[1] = det->dec->data.F64[j];
+            psVector *indices = psTreeAllWithin(tree, coords, MATCH_RADIUS); // Indices for matching sources
+            if (!indices) {
+                psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to search for matches");
+                psFree(coords);
+                psFree(tree);
+                psFree(merged);
+                return NULL;
+            }
+            if (indices->n == 0) {
+                psFree(indices);
+                continue;
+            }
+
+            // Which one do we keep?
+            float bestDistance = INFINITY; // Best distance to centre
+            long bestIndex = -1;           // Index with best distance
+            for (int k = 0; k < indices->n; k++) {
+                long index = indices->data.S64[k]; // Index of point
+                float distance = mergeDistance(merged, index); // Distance to centre of image
+                if (distance < bestDistance) {
+                    bestDistance = distance;
+                    bestIndex = index;
+                }
+            }
+
+            float distance = mergeDistance(det, j); // Distance to centre of image
+            if (distance < bestDistance) {
+                // Blow away existing sources
+                for (int k = 0; k < indices->n; k++) {
+                    long index = indices->data.S64[k]; // Index of point
+                    merged->mask->data.U8[index] = 0xFF;
+                }
+
+                ppMopsDetectionsCopySingle(merged, det, j);
+            }
+            psFree(indices);
+        }
+
+        psFree(tree);
+        ppMopsDetectionsPurge(merged);
+    }
+
+    merged->seeing /= num;
+
+    return merged;
 }
+