Index: trunk/ppTranslate/src/ppMopsMerge.c
===================================================================
--- trunk/ppTranslate/src/ppMopsMerge.c	(revision 28220)
+++ trunk/ppTranslate/src/ppMopsMerge.c	(revision 28623)
@@ -9,5 +9,5 @@
 #include "ppMops.h"
 
-#define LEAF_SIZE 2                     // Size of leaf
+#define LEAF_SIZE 4                     // Size of leaf
 #define MATCH_RADIUS SEC_TO_RAD(1.0)    // Matching radius
 #define MJD_TOL 1.0/3600.0/24.0         // Tolerance for MJD matching
@@ -17,12 +17,4 @@
 #define AIRMASS_TOL 1.0e-3              // Tolerance for airmass matching
 
-#if 0
-#undef psTrace
-#define psTrace(facil, level, ...)              \
-    if (level <= 5) {                           \
-        fprintf(stderr, __VA_ARGS__);           \
-    }
-#endif
-
 // Get distance from detection to centre of image
 static float mergeDistance(const ppMopsDetections *detections, // Detections of interest
@@ -30,45 +22,19 @@
     )
 {
-    float dx = detections->x->data.F32[index] - detections->naxis1 / 2.0;
-    float dy = detections->y->data.F32[index] - detections->naxis2 / 2.0;
+    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);
 }
 
 
-bool ppMopsPurgeDuplicates(const psArray *detections)
+ppMopsDetections *ppMopsMerge(const psArray *detections)
 {
     PS_ASSERT_ARRAY_NON_NULL(detections, NULL);
 
-    int numInputs = detections->n;                // Number of inputs
-    psTrace("ppMops.merge", 1, "Checking detections from %d inputs\n", numInputs);
+    psTrace("ppMops.merge", 1, "Merging detections from %ld inputs\n", detections->n);
 
-    long total = 0;                                // Total number of sources
-    psArray *duplicates = psArrayAlloc(numInputs); // Vector of duplicate bits for each input
-    psVector *dupNum = psVectorAlloc(numInputs, PS_TYPE_U32); // Number of duplicates for each input
-    psVectorInit(dupNum, 0);
-    for (int i = 0; i < numInputs; i++) {
-        ppMopsDetections *det = detections->data[i]; // Detections from
-        if (!det || det->num == 0) {
-            continue;
-        }
-        psVector *dupes = duplicates->data[i] = psVectorAlloc(det->num, PS_TYPE_U8);
-        psVectorInit(dupes, 0);
-        total += det->num;
-    }
-    psTrace("ppMops.merge", 2, "Total detections: %ld\n", total);
-
-    psVector *raMerged = psVectorAlloc(total, PS_TYPE_F64); // Merged RAs
-    psVector *decMerged = psVectorAlloc(total, PS_TYPE_F64); // Merged Decs
-    psVector *sourceMerged = psVectorAlloc(total, PS_TYPE_U16); // Source image of merged sources
-    psVector *indexMerged = psVectorAlloc(total, PS_TYPE_U32); // Index of merged sources
-    long num = 0;                                                   // Number of merged sources
-
-    const char *raBoresight = NULL, *decBoresight = NULL; // Boresight coordinates
-    const char *filter = NULL;                    // Filter name
-    float airmass = NAN, exptime = NAN;           // Exposure details
-    double posangle = NAN, alt = NAN, az = NAN; // Telescope pointing
-    double mjd = NAN;                           // Time of exposure
-
-    for (int i = 0; i < numInputs; i++) {
+    ppMopsDetections *merged = NULL;    // Merged list
+    int num = 1;                                                         // Number of merged files
+    for (int i = 0; i < detections->n; i++) {
         ppMopsDetections *det = detections->data[i]; // Detections of interest
         if (!det) {
@@ -79,90 +45,69 @@
             continue;
         }
+        num++;
+        if (!merged) {
+            psTrace("ppMops.merge", 3, "Accepting %ld detections from input %d\n", det->num, i);
+            merged = psMemIncrRefCounter(det);
+            continue;
+        }
+        psTrace("ppMops.merge", 3, "Merging %ld detections from input %d\n", det->num, i);
 
-        // Check exposure characteristics
-        if (num == 0) {
-            raBoresight = det->raBoresight;
-            decBoresight = det->decBoresight;
-            filter = det->filter;
-            airmass = det->airmass;
-            exptime = det->exptime;
-            posangle = det->posangle;
-            alt = det->alt;
-            az = det->az;
-        } else {
-            if (strcmp(raBoresight, det->raBoresight) != 0) {
-                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure RA values differ: %s vs %s",
-                        raBoresight, det->raBoresight);
-                return false;
-            }
-            if (strcmp(decBoresight, det->decBoresight) != 0) {
-                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure Dec values differ: %s vs %s",
-                        decBoresight, det->decBoresight);
-                return false;
-            }
-            if (strcmp(filter, det->filter) != 0) {
-                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure filter values differ: %s vs %s",
-                        filter, det->filter);
-                return false;
-            }
-            if (fabsf(airmass - det->airmass) > AIRMASS_TOL) {
-                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure airmass values differ: %f vs %f",
-                        airmass, det->airmass);
-                return false;
-            }
-            if (fabsf(exptime - det->exptime) > EXPTIME_TOL) {
-                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure exposure time values differ: %f vs %f",
-                        exptime, det->exptime);
-                return false;
-            }
-            if (fabs(posangle - det->posangle) > POSANGLE_TOL) {
-                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure position angle values differ: %f vs %f",
-                        posangle, det->posangle);
-                return false;
-            }
-            if (fabs(alt - det->alt) > BORESIGHT_TOL) {
-                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure altitude values differ: %lf vs %lf",
-                        alt, det->alt);
-                return false;
-            }
-            if (fabs(az - det->az) > BORESIGHT_TOL) {
-                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure azimuth values differ: %lf vs %lf",
-                        az, det->az);
-                return false;
-            }
-            if (fabs(mjd - det->mjd) > MJD_TOL) {
-                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure MJD values differ: %lf vs %lf",
-                        mjd, det->mjd);
-                return false;
-            }
+        // 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;
         }
 
-        psTrace("ppMops.merge", 3, "Accepting %ld detections from input %d\n", det->num, i);
-        memcpy(&raMerged->data.F64[num], det->ra->data.F64, det->num * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
-        memcpy(&decMerged->data.F64[num], det->dec->data.F64, det->num * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
-        for (long j = 0, k = num; j < det->num; j++, k++) {
-            sourceMerged->data.U16[k] = i;
-            indexMerged->data.U32[k] = j;
+        if (fabsf(merged->airmass - det->airmass) > AIRMASS_TOL) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure airmass values differ: %f vs %f",
+                    merged->airmass, det->airmass);
+            return NULL;
         }
-        num += det->num;
-    }
+        if (fabsf(merged->exptime - det->exptime) > EXPTIME_TOL) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure exposure time values differ: %f vs %f",
+                    merged->exptime, det->exptime);
+            return NULL;
+        }
+        if (fabs(merged->posangle - det->posangle) > POSANGLE_TOL) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure position angle values differ: %f vs %f",
+                    merged->posangle, det->posangle);
+            return NULL;
+        }
+        if (fabs(merged->alt - det->alt) > BORESIGHT_TOL) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure altitude values differ: %lf vs %lf",
+                    merged->alt, det->alt);
+            return NULL;
+        }
+        if (fabs(merged->az - det->az) > BORESIGHT_TOL) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure azimuth values differ: %lf vs %lf",
+                    merged->az, det->az);
+            return NULL;
+        }
+        if (fabs(merged->mjd - det->mjd) > MJD_TOL) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure MJD values differ: %lf vs %lf",
+                    merged->mjd, det->mjd);
+            return NULL;
+        }
 
-    psTrace("ppMops.merge", 3, "Generating kd-tree from %ld sources\n", num);
-    psTree *tree = psTreePlant(2, LEAF_SIZE, PS_TREE_SPHERICAL, raMerged, decMerged); // kd tree
-    if (!tree) {
-        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate kd tree");
-        return false;
-    }
+        merged->seeing += det->seeing;  // Taking average
 
-    for (int i = 0; i < detections->n; i++) {
-        ppMopsDetections *det = detections->data[i]; // Detections of interest
-        if (!det) {
-            continue;
-        } else if (det->num == 0) {
-            continue;
+        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;
         }
-        psTrace("ppMops.merge", 3, "Checking %ld detections from input %d\n", det->num, i);
 
-        psVector *dupes = duplicates->data[i];            // Duplicates list
         psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of interest
         for (int j = 0; j < det->num; j++) {
@@ -174,14 +119,14 @@
                 psFree(coords);
                 psFree(tree);
-                return false;
+                psFree(merged);
+                return NULL;
+            }
+            if (indices->n == 0) {
+                psTrace("ppMops.merge", 9, "No matches for source %d in input %d\n", j, i);
+                psFree(indices);
+                ppMopsDetectionsCopySingle(merged, det, j);
+                continue;
             }
             psTrace("ppMops.merge", 5, "%ld matches for source %d from input %d\n", indices->n, j, i);
-            psAssert(indices->n > 0, "Expect at least one match for source %d in input %d", j, i);
-
-            if (indices->n == 1 && sourceMerged->data.U16[indices->data.S64[0]] == i) {
-                // It's myself
-                psFree(indices);
-                continue;
-            }
 
             // Which one do we keep?
@@ -189,16 +134,6 @@
             long bestIndex = -1;           // Index with best distance
             for (int k = 0; k < indices->n; k++) {
-                long mergeIndex = indices->data.S64[k]; // Index of point in merged list
-                int source = sourceMerged->data.U16[mergeIndex]; // Source image
-                if (source == i) {
-                    continue;
-                }
-                long index = indexMerged->data.U32[mergeIndex];  // Index in source
-                psVector *dupes = duplicates->data[source];     // Duplicates list
-                if (dupes->data.U8[index]) {
-                    continue;
-                }
-
-                float distance = mergeDistance(detections->data[source], index); // Distance to centre of image
+                long index = indices->data.S64[k]; // Index of point
+                float distance = mergeDistance(merged, index); // Distance to centre of image
                 if (distance < bestDistance) {
                     bestDistance = distance;
@@ -207,84 +142,30 @@
             }
 
-            if (bestIndex >= 0) {
-                float distance = mergeDistance(det, j); // Distance to centre of image
-                if (bestIndex >= 0 && distance < bestDistance) {
-                    psTrace("ppMops.merge", 6, "New source clobbers old sources\n");
-                    for (int k = 0; k < indices->n; k++) {
-                        long mergeIndex = indices->data.S64[k]; // Index of point
-                        int source = sourceMerged->data.U16[mergeIndex]; // Source image
-                        if (source == i) {
-                            continue;
-                        }
-                        long index = indexMerged->data.U32[mergeIndex];  // Index in source
-                        psVector *dupes = duplicates->data[source];      // Duplicates list
-                        if (!dupes->data.U8[index]) {
-                            dupes->data.U8[index] = 0xFF;
-                            dupNum->data.U32[source]++;
-                        }
-                    }
-                } else {
-                    psTrace("ppMops.merge", 6, "Old sources clobber new source\n");
-                    dupes->data.U8[j] = 0xFF;
-                    dupNum->data.U32[i]++;
+            float distance = mergeDistance(det, j); // Distance to centre of image
+            if (distance < bestDistance) {
+                psTrace("ppMops.merge", 6, "New source clobbers old sources\n");
+                // 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);
+            } else {
+                psTrace("ppMops.merge", 6, "Old sources clobber new source\n");
             }
             psFree(indices);
         }
-        psFree(coords);
-    }
-    psFree(tree);
-    psFree(raMerged);
-    psFree(decMerged);
-    psFree(sourceMerged);
-    psFree(indexMerged);
 
-    // Remove duplicates
-    for (int i = 0; i < detections->n; i++) {
-        ppMopsDetections *det = detections->data[i]; // Detections of interest
-        if (!det) {
-            continue;
-        } else if (det->num == 0) {
-            continue;
-        }
-        psTrace("ppMops.merge", 3, "Purging %d duplicates from input %d\n", dupNum->data.U32[i], i);
+        psTrace("ppMops.merge", 3, "Done merging input %d, %ld merged sources\n", i, merged->num);
 
-#define VECTOR_PURGE_CASE(TYPE)                                      \
-        case PS_TYPE_##TYPE:                                         \
-          for (int i = 0, j = 0; i < vector->n; i++) {               \
-              if (!dupes->data.U8[i]) {                              \
-                  if (i == j) {                                      \
-                      j++;                                           \
-                      continue;                                      \
-                  }                                                  \
-                  vector->data.TYPE[j] = vector->data.TYPE[i];       \
-              }                                                      \
-          }                                                          \
-          break;
-
-        psVector *dupes = duplicates->data[i]; // Duplicates
-        psArray *table = psListToArray(det->table->list); // Table of data
-        for (int j = 0; j < table->n; j++) {
-            psMetadataItem *item = table->data[j]; // Table item
-            psAssert(item->type == PS_DATA_VECTOR, "Table column is not a vector: %x", item->type);
-            psVector *vector = item->data.V; // Vector to purge
-            switch (vector->type.type) {
-                VECTOR_PURGE_CASE(U8);
-                VECTOR_PURGE_CASE(U16);
-                VECTOR_PURGE_CASE(U32);
-                VECTOR_PURGE_CASE(U64);
-                VECTOR_PURGE_CASE(S8);
-                VECTOR_PURGE_CASE(S16);
-                VECTOR_PURGE_CASE(S32);
-                VECTOR_PURGE_CASE(S64);
-                VECTOR_PURGE_CASE(F32);
-                VECTOR_PURGE_CASE(F64);
-              default:
-                psAbort("Unrecognised vector type: %x", vector->type.type);
-            }
-        }
+        psFree(tree);
+        ppMopsDetectionsPurge(merged);
     }
 
-    return true;
+    psTrace("ppMops.merge", 2, "%ld sources in merged detections list\n", merged->num);
+
+    merged->seeing /= num;
+
+    return merged;
 }