- Timestamp:
- Mar 5, 2012, 5:19:48 PM (14 years ago)
- Location:
- branches/meh_branches/ppstack_test
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
ppTranslate/src/ppMopsMerge.c (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/meh_branches/ppstack_test
- Property svn:mergeinfo changed
-
branches/meh_branches/ppstack_test/ppTranslate/src/ppMopsMerge.c
r30656 r33415 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 = (float) (detections->x->data.F32[index] - detections->naxis1->data.S32[index] / 2.0);25 float dy = (float) (detections->y->data.F32[index] - detections->naxis2->data.S32[index] / 2.0);26 return PS_SQR(dx) + PS_SQR(dy);32 float dx = detections->x->data.F32[index] - detections->naxis1 / 2.0; 33 float dy = detections->y->data.F32[index] - detections->naxis2 / 2.0; 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 = 0; // 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++) { 169 // XXX: added by Bill 170 if (det->mask->data.U8[j]) { 171 // we marked this source marked as bad when we read it, mark it as a duplicate 172 dupes->data.U8[j] = 0xFF; 173 continue; 174 } 175 114 176 coords->data.F64[0] = det->ra->data.F64[j]; 115 177 coords->data.F64[1] = det->dec->data.F64[j]; … … 119 181 psFree(coords); 120 182 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); 183 return false; 184 } 185 psTrace("ppMops.merge", 5, "%ld matches for source %d from input %d\n", indices->n, j, i); 186 psAssert(indices->n > 0, "Expect at least one match for source %d in input %d", j, i); 187 188 if (indices->n == 1 && sourceMerged->data.U16[indices->data.S64[0]] == i) { 189 // It's myself 126 190 psFree(indices); 127 ppMopsDetectionsCopySingle(merged, det, j);128 191 continue; 129 192 } 130 psTrace("ppMops.merge", 5, "%ld matches for source %d from input %d\n", indices->n, j, i);131 193 132 194 // Which one do we keep? … … 134 196 long bestIndex = -1; // Index with best distance 135 197 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 198 long mergeIndex = indices->data.S64[k]; // Index of point in merged list 199 int source = sourceMerged->data.U16[mergeIndex]; // Source image 200 if (source == i) { 201 continue; 202 } 203 long index = indexMerged->data.U32[mergeIndex]; // Index in source 204 psVector *dupes = duplicates->data[source]; // Duplicates list 205 if (dupes->data.U8[index]) { 206 continue; 207 } 208 209 float distance = mergeDistance(detections->data[source], index); // Distance to centre of image 138 210 if (distance < bestDistance) { 139 211 bestDistance = distance; … … 142 214 } 143 215 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; 216 if (bestIndex >= 0) { 217 float distance = mergeDistance(det, j); // Distance to centre of image 218 if (bestIndex >= 0 && distance < bestDistance) { 219 psTrace("ppMops.merge", 6, "New source clobbers old sources\n"); 220 for (int k = 0; k < indices->n; k++) { 221 long mergeIndex = indices->data.S64[k]; // Index of point 222 int source = sourceMerged->data.U16[mergeIndex]; // Source image 223 if (source == i) { 224 continue; 225 } 226 long index = indexMerged->data.U32[mergeIndex]; // Index in source 227 psVector *dupes = duplicates->data[source]; // Duplicates list 228 if (!dupes->data.U8[index]) { 229 dupes->data.U8[index] = 0xFF; 230 dupNum->data.U32[source]++; 231 } 232 } 233 } else { 234 psTrace("ppMops.merge", 6, "Old sources clobber new source\n"); 235 dupes->data.U8[j] = 0xFF; 236 dupNum->data.U32[i]++; 151 237 } 152 ppMopsDetectionsCopySingle(merged, det, j);153 } else {154 psTrace("ppMops.merge", 6, "Old sources clobber new source\n");155 238 } 156 239 psFree(indices); 157 240 } 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 if (num == 0) { 166 //All detections were NULL?! 167 psTrace("ppMops.merge", 3, "All %ld detections were NULL\n", detections->n); 168 return NULL; 169 } 170 psTrace("ppMops.merge", 2, "%ld sources in merged detections list\n", merged->num); 171 172 merged->seeing /= (float) num; 173 174 return merged; 241 psFree(coords); 242 } 243 psFree(tree); 244 psFree(raMerged); 245 psFree(decMerged); 246 psFree(sourceMerged); 247 psFree(indexMerged); 248 249 // Remove duplicates 250 for (int i = 0; i < detections->n; i++) { 251 ppMopsDetections *det = detections->data[i]; // Detections of interest 252 if (!det) { 253 continue; 254 } else if (det->num == 0) { 255 continue; 256 } 257 psTrace("ppMops.merge", 3, "Purging %d duplicates from input %d\n", dupNum->data.U32[i], i); 258 259 #define VECTOR_PURGE_CASE(TYPE) \ 260 case PS_TYPE_##TYPE: { \ 261 long j = 0; \ 262 for (long i = 0; i < vector->n; i++) { \ 263 if (!dupes->data.U8[i]) { \ 264 if (i == j) { \ 265 j++; \ 266 continue; \ 267 } \ 268 vector->data.TYPE[j++] = vector->data.TYPE[i]; \ 269 } \ 270 } \ 271 vector->n = j; \ 272 } \ 273 break 274 275 psVector *dupes = duplicates->data[i]; // Duplicates 276 psArray *table = psListToArray(det->table->list); // Table of data 277 long newLength = -1; 278 for (int t = 0; t < table->n; t++) { 279 psMetadataItem *item = table->data[t]; // Table item 280 psAssert(item->type == PS_DATA_VECTOR, "Table column is not a vector: %x", item->type); 281 psVector *vector = item->data.V; // Vector to purge 282 switch (vector->type.type) { 283 VECTOR_PURGE_CASE(U8); 284 VECTOR_PURGE_CASE(U16); 285 VECTOR_PURGE_CASE(U32); 286 VECTOR_PURGE_CASE(U64); 287 VECTOR_PURGE_CASE(S8); 288 VECTOR_PURGE_CASE(S16); 289 VECTOR_PURGE_CASE(S32); 290 VECTOR_PURGE_CASE(S64); 291 VECTOR_PURGE_CASE(F32); 292 VECTOR_PURGE_CASE(F64); 293 default: 294 psAbort("Unrecognised vector type: %x", vector->type.type); 295 } 296 if (newLength == -1) { 297 newLength = vector->n; 298 } else if (vector->n != newLength) { 299 psAbort("Unexpected new length found : %ld expected: %ld", 300 vector->n, newLength); 301 } 302 303 } 304 // XXX IS this safe? Perhaps save in numGood? 305 det->num = newLength; 306 psFree(table); 307 } 308 psFree(dupNum); 309 psFree(duplicates); 310 311 return true; 175 312 } 176 313
Note:
See TracChangeset
for help on using the changeset viewer.
