Changeset 5743 for trunk/stac/src/stac.c
- Timestamp:
- Dec 7, 2005, 2:57:14 PM (21 years ago)
- File:
-
- 1 edited
-
trunk/stac/src/stac.c (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/stac/src/stac.c
r3673 r5743 18 18 // Set trace levels 19 19 psTraceSetLevel(".",0); 20 psTraceSetLevel("stac",10); 20 21 psTraceSetLevel("stac.checkMemory",10); 21 22 psTraceSetLevel("stac.config",10); … … 38 39 39 40 // Read input files 40 psArray *inputs = stacReadImages(config->inputs); 41 psArray *headers = NULL; // Array of headers 42 psArray *inputs = stacReadImages(&headers, config->inputs); 43 44 // Generate masks 45 psArray *masks = psArrayAlloc(inputs->n); 46 for (int i = 0; i < inputs->n; i++) { 47 psImage *image = inputs->data[i]; // Image for which to get mask 48 psImage *mask = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8); 49 for (int y = 0; y < image->numRows; y++) { 50 for (int x = 0; x < image->numCols; x++) { 51 if (image->data.F32[y][x] <= config->bad) { 52 mask->data.U8[y][x] = 1; 53 } else { 54 mask->data.U8[y][x] = 0; 55 } 56 } 57 } 58 masks->data[i] = mask; 59 } 41 60 42 61 // Read maps … … 47 66 // Get size, if not input 48 67 if (config->outnx == 0 || config->outny == 0) { 49 stacSize(&config->outnx, &config->outny, &config->xMapDiff, &config->yMapDiff, inputs, maps); 68 stacSize(&config->outnx, &config->outny, &config->xMapDiff, &config->yMapDiff, inputs, maps); 69 } 70 71 // Fix up the headers 72 for (int i = 0; i < headers->n; i++) { 73 psMetadata *header = headers->data[i]; 74 psMetadataAddS32(header, PS_LIST_TAIL, "NAXIS1", PS_META_REPLACE, "Number of pixels in x", 75 config->outnx); 76 psMetadataAddS32(header, PS_LIST_TAIL, "NAXIS2", PS_META_REPLACE, "Number of pixels in y", 77 (int)config->outny); 78 psMetadataAddS32(header, PS_LIST_TAIL, "BITPIX", PS_META_REPLACE, "Bits per pixel", -32); 79 #if 0 80 bool mdok = true; // Result of MD lookup 81 float crpix1 = psMetadataLookupF32(&mdok, header, "CRPIX1"); 82 if (mdok && !isnan(crpix1)) { 83 psMetadataAddF32(header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, 84 "WCS Coordinate reference pixel", crpix1 - 85 #endif 50 86 } 51 87 … … 58 94 // Write error images out to check 59 95 for (int i = 0; i < inputs->n; i++) { 60 char errName[MAXCHAR];// Filename of error image61 62 psFits *errorFile = psFitsAlloc(errName);63 if (!psFitsWriteImage(errorFile, NULL, errors->data[i], 0, NULL)) {64 65 66 67 psFree(errorFile);96 char errName[MAXCHAR]; // Filename of error image 97 sprintf(errName,"%s.err",config->inputs->data[i]); 98 psFits *errorFile = psFitsOpen(errName, "w"); 99 if (!psFitsWriteImage(errorFile, NULL, errors->data[i], 0)) { 100 psErrorStackPrint(stderr, "Unable to write image: %s\n", errName); 101 } 102 psTrace("stac", 1, "Error image written to %s\n", errName); 103 psFitsClose(errorFile); 68 104 } 69 105 #endif … … 72 108 psArray *transformed = NULL; 73 109 psArray *transformedErrors = NULL; 74 (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, NULL, NULL, NULL, NULL,75 110 (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, masks, NULL, NULL, NULL, 111 config->outnx, config->outny); 76 112 psTrace("stac.time",1,"Transformation completed at %f seconds\n", getTime() - startTime); 77 113 … … 79 115 // Write transformed images out to check 80 116 for (int i = 0; i < inputs->n; i++) { 81 char shiftName[MAXCHAR]; // Filename of shift image 82 char errName[MAXCHAR]; // Filename of error image 83 sprintf(shiftName,"%s.shift1",config->inputs->data[i]); 84 sprintf(errName,"%s.shifterr1",config->inputs->data[i]); 85 psFits *shiftFile = psFitsAlloc(shiftName); 86 psFits *errFile = psFitsAlloc(errName); 87 psImage *trans = transformed->data[i]; // Transformed image 88 psImage *transErr = transformedErrors->data[i]; // Transformed error image 89 if (!psFitsWriteImage(shiftFile, NULL, trans, 0, NULL)) { 90 psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName); 91 } 92 psTrace("stac", 1, "Shifted image written to %s\n", shiftName); 93 if (!psFitsWriteImage(errFile, NULL, transErr, 0, NULL)) { 94 psErrorStackPrint(stderr, "Unable to write image: %s\n", errName); 95 } 96 psTrace("stac", 1, "Shifted error image written to %s\n", errName); 97 psFree(shiftFile); 98 psFree(errFile); 99 } 100 #endif 101 102 103 104 psVector *scales = NULL; // Relative scales between images 105 psVector *offsets = NULL; // Offsets between images 117 char shiftName[MAXCHAR]; // Filename of shift image 118 char errName[MAXCHAR]; // Filename of error image 119 sprintf(shiftName,"%s.shift1",config->inputs->data[i]); 120 sprintf(errName,"%s.shifterr1",config->inputs->data[i]); 121 psFits *shiftFile = psFitsOpen(shiftName, "w"); 122 psFits *errFile = psFitsOpen(errName, "w"); 123 psImage *trans = transformed->data[i]; // Transformed image 124 psImage *transErr = transformedErrors->data[i]; // Transformed error image 125 if (!psFitsWriteImage(shiftFile, NULL, trans, 0)) { 126 psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName); 127 } 128 psTrace("stac", 1, "Shifted image written to %s\n", shiftName); 129 if (!psFitsWriteImage(errFile, NULL, transErr, 0)) { 130 psErrorStackPrint(stderr, "Unable to write image: %s\n", errName); 131 } 132 psTrace("stac", 1, "Shifted error image written to %s\n", errName); 133 psFitsClose(shiftFile); 134 psFitsClose(errFile); 135 } 136 #endif 137 138 psVector *scales = NULL; // Relative scales between images 139 psVector *offsets = NULL; // Offsets between images 106 140 (void)stacScales(&scales, &offsets, transformed, config->starFile, config->starMapFile, config->xMapDiff, 107 141 config->yMapDiff, config->aper); 108 142 // Rescale the images 109 143 (void)stacRescale(transformed, transformedErrors, NULL, scales, offsets); 110 144 // Set the saturation and bad values 111 145 psVector *saturated = psVectorAlloc(transformed->n, PS_TYPE_F32); // Saturation limits 112 psVector *bad = psVectorAlloc(transformed->n, PS_TYPE_F32); 146 psVector *bad = psVectorAlloc(transformed->n, PS_TYPE_F32); // Bad limits 113 147 for (int i = 0; i < transformed->n; i++) { 114 saturated->data.F32[i] = (config->saturated - offsets->data.F32[i]) / scales->data.F32[i]; 115 bad->data.F32[i] = (config->bad - offsets->data.F32[i]) / scales->data.F32[i]; 148 saturated->data.F32[i] = (config->saturated - offsets->data.F32[i]) / scales->data.F32[i]; 149 bad->data.F32[i] = (config->bad - offsets->data.F32[i]) / scales->data.F32[i]; 150 } 151 152 // Save shifted images 153 if (config->saveShifts) { 154 psImage *image = NULL; // Copy of image, to remove NAN for output 155 for (int i = 0; i < transformed->n; i++) { 156 char shiftName[MAXCHAR]; // Filename of shifted image 157 sprintf(shiftName, "%s.shift", config->inputs->data[i]); 158 psFits *shiftFile = psFitsOpen(shiftName, "w"); 159 image = psImageCopy(NULL, transformed->data[i], PS_TYPE_F32); 160 (void)psImageClipNaN(image, 0.0); 161 if (!psFitsWriteImage(shiftFile, headers->data[i], image, 0)) { 162 psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName); 163 } 164 psTrace("stac", 1, "Shifted image %d written to %s\n", i, shiftName); 165 psFitsClose(shiftFile); 166 psFree(image); 167 } 116 168 } 117 169 … … 120 172 psImage *combined = NULL; 121 173 (void)stacCombine(&combined, &rejected, transformed, transformedErrors, config->nReject, NULL, saturated, 122 174 bad, config->reject); 123 175 124 176 psTrace("stac.time",1,"First combination completed at %f seconds\n", getTime() - startTime); … … 128 180 // Write rejection images out to check 129 181 for (int i = 0; i < rejected->n; i++) { 130 char rejName[MAXCHAR];// Filename of rejection image131 132 133 psFits *rejFile = psFitsAlloc(rejName);134 if (!psFitsWriteImage(rejFile, NULL, rejected->data[i], 0, NULL)) {135 136 137 138 psFree(rejFile);182 char rejName[MAXCHAR]; // Filename of rejection image 183 sprintf(rejName, "%s.shiftrej", config->inputs->data[i]); 184 185 psFits *rejFile = psFitsOpen(rejName, "w"); 186 if (!psFitsWriteImage(rejFile, NULL, rejected->data[i], 0)) { 187 psErrorStackPrint(stderr, "Unable to write image: %s\n", rejName); 188 } 189 psTrace("stac", 1, "Rejection image written to %s\n", rejName); 190 psFitsClose(rejFile); 139 191 } 140 192 141 193 // Write out pre-combined image 142 char preName[MAXCHAR]; 194 char preName[MAXCHAR]; // Filename of precombined image 143 195 sprintf(preName, "%s.pre", config->output); 144 psFits *preFile = psFits Alloc(preName);145 if (!psFitsWriteImage(preFile, NULL, combined, 0 , NULL)) {146 196 psFits *preFile = psFitsOpen(preName, "w"); 197 if (!psFitsWriteImage(preFile, NULL, combined, 0)) { 198 psErrorStackPrint(stderr, "Unable to write image: %s\n", preName); 147 199 } 148 200 psTrace("stac", 1, "Pre-combined image written to %s\n", preName); 149 psF ree(preFile);201 psFitsClose(preFile); 150 202 #endif 151 203 … … 153 205 psArray *regions = psArrayAlloc(inputs->n); // Array of images denoting regions of interest 154 206 for (int i = 0; i < inputs->n; i++) { 155 156 157 158 #ifdef TESTING 159 char regionName[MAXCHAR];// Filename of region image160 161 psFits *regionFile = psFitsAlloc(regionName);162 if (!psFitsWriteImage(regionFile, NULL, regions->data[i], 0, NULL)) {163 164 165 166 psFree(regionFile);207 regions->data[i] = stacAreaOfInterest(rejected->data[i], inverseMaps->data[i], 208 ((psImage*)(inputs->data[i]))->numCols, 209 ((psImage*)(inputs->data[i]))->numRows); 210 #ifdef TESTING 211 char regionName[MAXCHAR]; // Filename of region image 212 sprintf(regionName,"%s.region",config->inputs->data[i]); 213 psFits *regionFile = psFitsOpen(regionName, "w"); 214 if (!psFitsWriteImage(regionFile, NULL, regions->data[i], 0)) { 215 psErrorStackPrint(stderr, "Unable to write image: %s\n", regionName); 216 } 217 psTrace("stac", 1, "Region image written to %s\n", regionName); 218 psFitsClose(regionFile); 167 219 #endif 168 220 } … … 170 222 // Transform rejected pixels to source frame 171 223 psArray *rejectedSource = stacRejection(inputs, rejected, regions, maps, inverseMaps, config->frac, 172 224 config->grad, config->inputs); 173 225 174 226 // Get regions of interest in the output frame 175 227 psImage *combineRegion = psImageAlloc(config->outnx, config->outny, PS_TYPE_U8); 176 228 for (int y = 0; y < config->outny; y++) { 177 for (int x = 0; x < config->outnx; x++) { 178 combineRegion->data.U8[y][x] = 0; 179 } 180 } 181 for (int i = 0; i < inputs->n; i++) { 182 psImage *region = stacAreaOfInterest(rejectedSource->data[i], maps->data[i], config->outnx, 183 config->outny); 184 psBinaryOp(combineRegion, combineRegion, "+", region); 185 psFree(region); 229 for (int x = 0; x < config->outnx; x++) { 230 combineRegion->data.U8[y][x] = 0; 231 } 232 } 233 for (int i = 0; i < inputs->n; i++) { 234 psImage *region = stacAreaOfInterest(rejectedSource->data[i], maps->data[i], config->outnx, 235 config->outny); 236 psBinaryOp(combineRegion, combineRegion, "+", region); 237 psFree(region); 238 239 // Do OR of masks 240 psImage *mask = masks->data[i]; // Mask of input image (bad columns etc) 241 psImage *cr = rejectedSource->data[i]; // Mask of CRs 242 for (int y = 0; y < mask->numRows; y++) { 243 for (int x = 0; x < mask->numCols; x++) { 244 if (cr->data.U8[y][x] > 0) { 245 mask->data.U8[y][x] = 1; 246 } 247 } 248 } 186 249 } 187 250 188 251 // Redo transformation with the masks and scales/offsets 189 (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, rejectedSource,190 252 (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, masks, 253 combineRegion, scales, offsets, config->outnx, config->outny); 191 254 192 255 // Combine the newly-transformed CR-free images, no rejection … … 194 257 rejected = NULL; 195 258 (void)stacCombine(&combined, &rejected, transformed, transformedErrors, 0, combineRegion, saturated, 196 259 bad, config->reject); 197 260 198 261 // Write output image 199 psFits *outFile = psFitsAlloc(config->output); 200 if (!psFitsWriteImage(outFile, NULL, combined, 0, NULL)) { 201 psErrorStackPrint(stderr, "Unable to write image: %s\n", config->output); 262 psFits *outFile = psFitsOpen(config->output, "w"); 263 int numPix = psImageClipNaN(combined, 0.0); 264 if (numPix > 0) { 265 psTrace("stac", 3, "Clipping %d NaN pixels to zero.\n", numPix); 266 } 267 if (!psFitsWriteImage(outFile, headers->data[0], combined, 0)) { 268 psErrorStackPrint(stderr, "Unable to write image: %s\n", config->output); 202 269 } 203 270 psTrace("stac", 1, "Combined image written to %s\n", config->output); 204 psFree(outFile); 205 206 // Save shifted images 207 if (config->saveShifts) { 208 for (int i = 0; i < transformed->n; i++) { 209 char shiftName[MAXCHAR]; // Filename of shifted image 210 sprintf(shiftName, "%s.shift", config->inputs->data[i]); 211 psFits *shiftFile = psFitsAlloc(shiftName); 212 if (!psFitsWriteImage(shiftFile, NULL, transformed->data[i], 0, NULL)) { 213 psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName); 214 } 215 psTrace("stac", 1, "Shifted image %d written to %s\n", i, shiftName); 216 psFree(shiftFile); 217 } 218 } 271 psFitsClose(outFile); 219 272 220 273 // Free everything I've used … … 228 281 psFree(inverseMaps); 229 282 psFree(errors); 283 psFree(masks); 230 284 psFree(transformedErrors); 231 285 psFree(transformed); 232 286 psFree(rejectedSource); 233 287 psFree(combined); 288 psFree(saturated); 289 psFree(bad); 234 290 235 291 psTrace("stac.time",1,"Final combination completed at %f seconds\n", getTime() - startTime);
Note:
See TracChangeset
for help on using the changeset viewer.
