Changeset 28484 for branches/pap/ppTranslate/src/ppMopsMerge.c
- Timestamp:
- Jun 24, 2010, 2:59:09 PM (16 years ago)
- Location:
- branches/pap
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
ppTranslate/src/ppMopsMerge.c (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/pap
- Property svn:mergeinfo changed
-
branches/pap/ppTranslate/src/ppMopsMerge.c
r25256 r28484 9 9 #include "ppMops.h" 10 10 11 #define LEAF_SIZE 4// Size of leaf11 #define LEAF_SIZE 2 // Size of leaf 12 12 #define MATCH_RADIUS SEC_TO_RAD(1.0) // Matching radius 13 13 #define MJD_TOL 1.0/3600.0/24.0 // Tolerance for MJD matching … … 17 17 #define AIRMASS_TOL 1.0e-3 // Tolerance for airmass matching 18 18 19 #if 0 20 #undef psTrace 21 #define psTrace(facil, level, ...) \ 22 if (level <= 5) { \ 23 fprintf(stderr, __VA_ARGS__); \ 24 } 25 #endif 26 19 27 // Get distance from detection to centre of image 20 28 static float mergeDistance(const ppMopsDetections *detections, // Detections of interest … … 22 30 ) 23 31 { 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;32 float dx = detections->x->data.F32[index] - detections->naxis1 / 2.0; 33 float dy = detections->y->data.F32[index] - detections->naxis2 / 2.0; 26 34 return PS_SQR(dx) + PS_SQR(dy); 27 35 } 28 36 29 37 30 ppMopsDetections *ppMopsMerge(const psArray *detections)38 bool ppMopsPurgeDuplicates(const psArray *detections) 31 39 { 32 40 PS_ASSERT_ARRAY_NON_NULL(detections, NULL); 33 41 34 psTrace("ppMops.merge", 1, "Merging detections from %ld inputs\n", detections->n); 35 36 ppMopsDetections *merged = NULL; // Merged list 37 int num = 1; // Number of merged files 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 38 157 for (int i = 0; i < detections->n; i++) { 39 158 ppMopsDetections *det = detections->data[i]; // Detections of interest 40 159 if (!det) { 41 psTrace("ppMops.merge", 3, "Ignoring NULL input %d\n", i);42 160 continue; 43 161 } else if (det->num == 0) { 44 psTrace("ppMops.merge", 3, "Ignoring empty input %d\n", i); 45 continue; 46 } 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); 54 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; 70 } 71 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; 76 } 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 } 102 103 merged->seeing += det->seeing; // Taking average 104 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; 110 } 111 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 112 167 psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of interest 113 168 for (int j = 0; j < det->num; j++) { … … 119 174 psFree(coords); 120 175 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); 176 return false; 177 } 178 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 126 183 psFree(indices); 127 ppMopsDetectionsCopySingle(merged, det, j);128 184 continue; 129 185 } 130 psTrace("ppMops.merge", 5, "%ld matches for source %d from input %d\n", indices->n, j, i);131 186 132 187 // Which one do we keep? … … 134 189 long bestIndex = -1; // Index with best distance 135 190 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 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 138 203 if (distance < bestDistance) { 139 204 bestDistance = distance; … … 142 207 } 143 208 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; 151 } 152 ppMopsDetectionsCopySingle(merged, det, j); 153 } else { 154 psTrace("ppMops.merge", 6, "Old sources clobber new source\n"); 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]++; 230 } 155 231 } 156 232 psFree(indices); 157 233 } 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 165 psTrace("ppMops.merge", 2, "%ld sources in merged detections list\n", merged->num); 166 167 merged->seeing /= num; 168 169 return merged; 234 psFree(coords); 235 } 236 psFree(tree); 237 psFree(raMerged); 238 psFree(decMerged); 239 psFree(sourceMerged); 240 psFree(indexMerged); 241 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); 251 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 } 286 } 287 288 return true; 170 289 } 171 290
Note:
See TracChangeset
for help on using the changeset viewer.
