IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 21366


Ignore:
Timestamp:
Feb 5, 2009, 5:03:33 PM (17 years ago)
Author:
Paul Price
Message:

Merging pap_branch_20090128. Resolved a small number of conflicts. Compiles, but not tested in detail.

Location:
trunk
Files:
23 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppStack/src/ppStack.h

    r21092 r21366  
    3838    psArray *imageFits;                 // FITS file pointers for images
    3939    psArray *maskFits;                  // FITS file pointers for masks
    40     psArray *weightFits;                // FITS file pointers for weights
     40    psArray *varianceFits;                // FITS file pointers for variances
    4141} ppStackThreadData;
    4242
    … …  
    4545                                          const psArray *imageNames, // Names of images to read
    4646                                          const psArray *maskNames, // Names of masks to read
    47                                           const psArray *weightNames, // Names of weight maps to read
     47                                          const psArray *varianceNames, // Names of variance maps to read
    4848                                          const pmConfig *config // Configuration
    4949    );

  • trunk/ppStack/src/ppStackArguments.c

    r21199 r21366  
    2424            "\tIMAGE(STR):     Image filename\n"
    2525            "\tMASK(STR):      Mask filename\n"
    26             "\tWEIGHT(STR):    Weight map filename\n"
     26            "\tVARIANCE(STR):  Variance map filename\n"
    2727            "\tPSF(STR):       PSF filename\n"
    2828            "\tSOURCES(STR):   Sources filename\n"
    … …  
    151151    psMetadataAddStr(arguments, PS_LIST_TAIL, "-mask-poor", 0, "Mask value to give poor pixels", NULL);
    152152    psMetadataAddF32(arguments, PS_LIST_TAIL, "-threshold-mask", 0, "Threshold for mask deconvolution", NAN);
    153     psMetadataAddF32(arguments, PS_LIST_TAIL, "-poor-frac", 0, "Fraction of weight for poor pixels", NAN);
     153    psMetadataAddF32(arguments, PS_LIST_TAIL, "-poor-frac", 0, "Fraction of variance for poor pixels", NAN);
    154154    psMetadataAddF32(arguments, PS_LIST_TAIL, "-deconv-limit", 0, "Maximum deconvolution fraction limit", NAN);
    155155    psMetadataAddF32(arguments, PS_LIST_TAIL, "-image-rej", 0,
    … …  
    183183    psMetadataAddStr(arguments, PS_LIST_TAIL, "-temp-image", 0, "Suffix for temporary images", NULL);
    184184    psMetadataAddStr(arguments, PS_LIST_TAIL, "-temp-mask", 0, "Suffix for temporary masks", NULL);
    185     psMetadataAddStr(arguments, PS_LIST_TAIL, "-temp-weight", 0, "Suffix for temporary weight maps", NULL);
     185    psMetadataAddStr(arguments, PS_LIST_TAIL, "-temp-variance", 0, "Suffix for temporary variance maps", NULL);
    186186    psMetadataAddBool(arguments, PS_LIST_TAIL, "-temp-delete", 0,
    187187                      "Delete temporary files on completion?", false);
    … …  
    283283    VALUE_ARG_RECIPE_INT("-renorm-num", "RENORM.NUM", S32, 0);
    284284    VALUE_ARG_RECIPE_FLOAT("-renorm-width", "RENORM.WIDTH", F32);
    285     valueArgRecipeStr(arguments, recipe, "-renorm-mean", "RENORM.MEAN", recipe);
     285    valueArgRecipeStr(arguments, recipe, "-renorm-mean",  "RENORM.MEAN", recipe);
    286286    valueArgRecipeStr(arguments, recipe, "-renorm-stdev", "RENORM.STDEV", recipe);
    287287
    288     valueArgRecipeStr(arguments, recipe, "-temp-image",  "TEMP.IMAGE",  recipe);
    289     valueArgRecipeStr(arguments, recipe, "-temp-mask",   "TEMP.MASK",   recipe);
    290     valueArgRecipeStr(arguments, recipe, "-temp-weight", "TEMP.WEIGHT", recipe);
     288    valueArgRecipeStr(arguments, recipe, "-temp-image",    "TEMP.IMAGE",  recipe);
     289    valueArgRecipeStr(arguments, recipe, "-temp-mask",     "TEMP.MASK",   recipe);
     290    valueArgRecipeStr(arguments, recipe, "-temp-variance", "TEMP.VARIANCE", recipe);
    291291
    292292    if (psMetadataLookupBool(NULL, arguments, "-temp-delete") ||

  • trunk/ppStack/src/ppStackCamera.c

    r20995 r21366  
    7171bool ppStackCamera(pmConfig *config)
    7272{
    73     bool haveWeights = false;           // Do we have weight maps?
     73    bool haveVariances = false;           // Do we have variance maps?
    7474    bool havePSFs = false;              // Do we have PSFs?
    7575
    … …  
    9797        bool mdok;
    9898        psString mask = psMetadataLookupStr(&mdok, input, "MASK"); // Name of mask
    99         psString weight = psMetadataLookupStr(&mdok, input, "WEIGHT"); // Name of weight map
     99        psString variance = psMetadataLookupStr(&mdok, input, "VARIANCE"); // Name of variance map
    100100        psString psf = psMetadataLookupStr(&mdok, input, "PSF"); // Name of PSF
    101101        psString sources = psMetadataLookupStr(&mdok, input, "SOURCES"); // Name of sources
    … …  
    140140        }
    141141
    142         // Optionally add the weight file
    143         if (weight && strlen(weight) > 0) {
    144             haveWeights = true;
    145             psArray *weightFiles = psArrayAlloc(1); // Array of filenames for this FPA
    146             weightFiles->data[0] = psMemIncrRefCounter(weight);
    147             psMetadataAddArray(config->arguments, PS_LIST_TAIL, "WEIGHT.FILENAMES", PS_META_REPLACE,
    148                                "Filenames of weight files", weightFiles);
    149             psFree(weightFiles);
     142        // Optionally add the variance file
     143        if (variance && strlen(variance) > 0) {
     144            haveVariances = true;
     145            psArray *varianceFiles = psArrayAlloc(1); // Array of filenames for this FPA
     146            varianceFiles->data[0] = psMemIncrRefCounter(variance);
     147            psMetadataAddArray(config->arguments, PS_LIST_TAIL, "VARIANCE.FILENAMES", PS_META_REPLACE,
     148                               "Filenames of variance files", varianceFiles);
     149            psFree(varianceFiles);
    150150
    151151            bool status;
    152             pmFPAfile *weightFile = pmFPAfileBindFromArgs(&status, imageFile, config, "PPSTACK.INPUT.WEIGHT",
    153                                                           "WEIGHT.FILENAMES");
     152            pmFPAfile *varianceFile = pmFPAfileBindFromArgs(&status, imageFile, config,
     153                                                            "PPSTACK.INPUT.VARIANCE", "VARIANCE.FILENAMES");
    154154            if (!status) {
    155                 psError(PS_ERR_UNKNOWN, false, "Unable to define file from weight %d (%s)", i, weight);
    156                 return false;
    157             }
    158             if (weightFile->type != PM_FPA_FILE_WEIGHT) {
    159                 psError(PS_ERR_IO, true, "PPSTACK.INPUT.WEIGHT is not of type WEIGHT");
     155                psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)", i, variance);
     156                return false;
     157            }
     158            if (varianceFile->type != PM_FPA_FILE_VARIANCE) {
     159                psError(PS_ERR_IO, true, "PPSTACK.INPUT.VARIANCE is not of type VARIANCE");
    160160                return false;
    161161            }
    … …  
    217217#if 0
    218218        // Output convolved files
    219         pmFPAfile *outconvImage  = defineOutputConvolved("PPSTACK.OUTCONV", imageFile->fpa, config,
     219        pmFPAfile *outconvImage    = defineOutputConvolved("PPSTACK.OUTCONV", imageFile->fpa, config,
     220                                                           PM_FPA_FILE_IMAGE);
     221        pmFPAfile *outconvMask     = defineOutputConvolved("PPSTACK.OUTCONV.MASK", imageFile->fpa, config,
     222                                                           PM_FPA_FILE_MASK);
     223        pmFPAfile *outconvVariance = defineOutputConvolved("PPSTACK.OUTCONV.VARIANCE", imageFile->fpa, config,
     224                                                           PM_FPA_FILE_VARIANCE);
     225        if (!outconvImage || !outconvMask || !outconvVariance) {
     226            return false;
     227        }
     228
     229        // Input convolved files
     230        pmFPAfile *inconvImage    = defineInputConvolved("PPSTACK.INCONV", outconvImage, config,
    220231                                                         PM_FPA_FILE_IMAGE);
    221         pmFPAfile *outconvMask   = defineOutputConvolved("PPSTACK.OUTCONV.MASK", imageFile->fpa, config,
     232        pmFPAfile *inconvMask     = defineInputConvolved("PPSTACK.INCONV.MASK", outconvMask, config,
    222233                                                         PM_FPA_FILE_MASK);
    223         pmFPAfile *outconvWeight = defineOutputConvolved("PPSTACK.OUTCONV.WEIGHT", imageFile->fpa, config,
    224                                                          PM_FPA_FILE_WEIGHT);
    225         if (!outconvImage || !outconvMask || !outconvWeight) {
    226             return false;
    227         }
    228 
    229         // Input convolved files
    230         pmFPAfile *inconvImage  = defineInputConvolved("PPSTACK.INCONV", outconvImage, config,
    231                                                        PM_FPA_FILE_IMAGE);
    232         pmFPAfile *inconvMask   = defineInputConvolved("PPSTACK.INCONV.MASK", outconvMask, config,
    233                                                        PM_FPA_FILE_MASK);
    234         pmFPAfile *inconvWeight = defineInputConvolved("PPSTACK.INCONV.WEIGHT", outconvWeight, config,
    235                                                        PM_FPA_FILE_WEIGHT);
    236         if (!inconvImage || !inconvMask || !inconvWeight) {
     234        pmFPAfile *inconvVariance = defineInputConvolved("PPSTACK.INCONV.VARIANCE", outconvVariance, config,
     235                                                         PM_FPA_FILE_VARIANCE);
     236        if (!inconvImage || !inconvMask || !inconvVariance) {
    237237            return false;
    238238        }
    … …  
    246246        psMetadataRemoveKey(config->arguments, "MASK.FILENAMES");
    247247    }
    248     if (psMetadataLookup(config->arguments, "WEIGHT.FILENAMES")) {
    249         psMetadataRemoveKey(config->arguments, "WEIGHT.FILENAMES");
     248    if (psMetadataLookup(config->arguments, "VARIANCE.FILENAMES")) {
     249        psMetadataRemoveKey(config->arguments, "VARIANCE.FILENAMES");
    250250    }
    251251    if (psMetadataLookup(config->arguments, "PSF.FILENAMES")) {
    … …  
    295295    outMask->save = true;
    296296
    297     // Output weight
    298     if (haveWeights) {
    299         pmFPAfile *outWeight = pmFPAfileDefineOutput(config, output->fpa, "PPSTACK.OUTPUT.WEIGHT");
    300         if (!outWeight) {
    301             psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.OUTPUT.WEIGHT"));
    302             return false;
    303         }
    304         if (outWeight->type != PM_FPA_FILE_WEIGHT) {
    305             psError(PS_ERR_IO, true, "PPSTACK.OUTPUT.WEIGHT is not of type WEIGHT");
    306             return false;
    307         }
    308         outWeight->save = true;
     297    // Output variance
     298    if (haveVariances) {
     299        pmFPAfile *outVariance = pmFPAfileDefineOutput(config, output->fpa, "PPSTACK.OUTPUT.VARIANCE");
     300        if (!outVariance) {
     301            psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.OUTPUT.VARIANCE"));
     302            return false;
     303        }
     304        if (outVariance->type != PM_FPA_FILE_VARIANCE) {
     305            psError(PS_ERR_IO, true, "PPSTACK.OUTPUT.VARIANCE is not of type VARIANCE");
     306            return false;
     307        }
     308        outVariance->save = true;
    309309    }
    310310

  • trunk/ppStack/src/ppStackLoop.c

    r21258 r21366  
    2424
    2525// Files required for the convolution
    26 static char *convolveFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.WEIGHT", NULL };
     26static char *convolveFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.VARIANCE", NULL };
    2727
    2828// Output files for the combination
    29 static char *combineFiles[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT",
     29static char *combineFiles[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.VARIANCE",
    3030                                "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2", NULL };
    3131
    … …  
    204204    const char *tempImage = psMetadataLookupStr(NULL, recipe, "TEMP.IMAGE"); // Suffix for temporary images
    205205    const char *tempMask = psMetadataLookupStr(NULL, recipe, "TEMP.MASK"); // Suffix for temporary masks
    206     const char *tempWeight = psMetadataLookupStr(NULL, recipe, "TEMP.WEIGHT"); // Suffix for temp weight maps
    207     if (!tempImage || !tempMask || !tempWeight) {
     206    const char *tempVariance = psMetadataLookupStr(NULL, recipe, "TEMP.VARIANCE"); // Suffix for temp variance maps
     207    if (!tempImage || !tempMask || !tempVariance) {
    208208        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    209                 "Unable to find TEMP.IMAGE, TEMP.MASK and TEMP.WEIGHT in recipe");
     209                "Unable to find TEMP.IMAGE, TEMP.MASK and TEMP.VARIANCE in recipe");
    210210        return false;
    211211    }
    … …  
    387387    psArray *imageNames = psArrayAlloc(num);
    388388    psArray *maskNames = psArrayAlloc(num);
    389     psArray *weightNames = psArrayAlloc(num);
     389    psArray *varianceNames = psArrayAlloc(num);
    390390    for (int i = 0; i < num; i++) {
    391         psString imageName = NULL, maskName = NULL, weightName = NULL; // Names for convolved images
     391        psString imageName = NULL, maskName = NULL, varianceName = NULL; // Names for convolved images
    392392        psStringAppend(&imageName, "%s/%s.%d.%s", tempDir, tempName, i, tempImage);
    393393        psStringAppend(&maskName, "%s/%s.%d.%s", tempDir, tempName, i, tempMask);
    394         psStringAppend(&weightName, "%s/%s.%d.%s", tempDir, tempName, i, tempWeight);
    395         psTrace("ppStack", 5, "Temporary files: %s %s %s\n", imageName, maskName, weightName);
     394        psStringAppend(&varianceName, "%s/%s.%d.%s", tempDir, tempName, i, tempVariance);
     395        psTrace("ppStack", 5, "Temporary files: %s %s %s\n", imageName, maskName, varianceName);
    396396        imageNames->data[i] = imageName;
    397397        maskNames->data[i] = maskName;
    398         weightNames->data[i] = weightName;
     398        varianceNames->data[i] = varianceName;
    399399    }
    400400    // Free the outputName that we grabbed above to set up tempName.
    … …  
    417417    psList *fpaList = psListAlloc(NULL); // List of input FPAs, for concept averaging
    418418    psList *cellList = psListAlloc(NULL); // List of input cells, for concept averaging
     419    psArray *covariances = psArrayAlloc(num); // Covariance matrices
    419420    for (int i = 0; i < num; i++) {
    420421        psTrace("ppStack", 2, "Convolving input %d of %d to target PSF....\n", i, num);
    … …  
    431432            psFree(fpaList);
    432433            psFree(cellList);
     434            psFree(covariances);
    433435            return false;
    434436        }
    … …  
    450452            psFree(fpaList);
    451453            psFree(cellList);
     454            psFree(covariances);
    452455            return false;
    453456        }
    … …  
    464467            continue;
    465468        }
     469        covariances->data[i] = psMemIncrRefCounter(readout->covariance);
    466470
    467471        if (stats) {
    … …  
    486490        pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image
    487491        assert(hdu);
    488         writeImage(imageNames->data[i],  hdu->header, readout->image, config);
    489492        psMetadata *maskHeader = psMetadataCopy(NULL, hdu->header); // Copy of header, for mask
    490493        pmConfigMaskWriteHeader(config, maskHeader);
    491         writeImage(maskNames->data[i],   maskHeader, readout->mask, config);
     494        writeImage(maskNames->data[i], maskHeader, readout->mask, config);
    492495        psFree(maskHeader);
    493         writeImage(weightNames->data[i], hdu->header, readout->weight, config);
     496        writeImage(varianceNames->data[i], hdu->header, readout->variance, config);
    494497
    495498        pmCell *inCell = readout->parent; // Input cell
    … …  
    502505            psFree(fpaList);
    503506            psFree(cellList);
     507            psFree(covariances);
    504508            return false;
    505509        }
    … …  
    523527        psFree(fpaList);
    524528        psFree(cellList);
     529        psFree(covariances);
    525530        return false;
    526531    }
    … …  
    570575            psFree(matchChi2);
    571576            psFree(values);
     577            psFree(covariances);
    572578            return false;
    573579        }
    … …  
    611617        psFree(inputMask);
    612618        psFree(matchChi2);
     619        psFree(covariances);
    613620        return false;
    614621    }
    … …  
    622629    // Start threading
    623630    ppStackThreadInit();
    624     ppStackThreadData *stack = ppStackThreadDataSetup(cells, imageNames, maskNames, weightNames, config);
     631    ppStackThreadData *stack = ppStackThreadDataSetup(cells, imageNames, maskNames, varianceNames, config);
    625632
    626633    psTimerStart("PPSTACK_INITIAL");
    … …  
    645652            psFree(matchChi2);
    646653            psFree(cells);
     654            psFree(covariances);
    647655            return false;
    648656        }
    … …  
    658666            psFree(matchChi2);
    659667            psFree(cells);
     668            psFree(covariances);
    660669            return false;
    661670        }
    … …  
    676685            psFree(outRO);
    677686            psFree(cells);
     687            psFree(covariances);
    678688            return false;
    679689        }
    … …  
    695705                psFree(view);
    696706                psFree(outRO);
     707                psFree(covariances);
    697708                return false;
    698709            }
    … …  
    721732                psFree(view);
    722733                psFree(outRO);
     734                psFree(covariances);
    723735                return false;
    724736            }
    … …  
    735747            psFree(view);
    736748            psFree(outRO);
     749            psFree(covariances);
    737750            return false;
    738751        }
    … …  
    803816                psFree(inspect);
    804817                psFree(rejected);
     818                psFree(covariances);
    805819                return false;
    806820            }
    … …  
    818832            psFree(inspect);
    819833            psFree(rejected);
     834            psFree(covariances);
    820835            return false;
    821836        }
    … …  
    912927            psFree(view);
    913928            psFree(outRO);
     929            psFree(covariances);
    914930            return false;
    915931        }
    … …  
    941957                psFree(view);
    942958                psFree(outRO);
     959                psFree(covariances);
    943960                return false;
    944961            }
    … …  
    962979                psFree(view);
    963980                psFree(outRO);
     981                psFree(covariances);
    964982                return false;
    965983            }
    … …  
    974992            psFree(view);
    975993            psFree(outRO);
     994            psFree(covariances);
    976995            return false;
    977996        }
    … …  
    979998
    980999    memDump("final");
     1000
     1001    // Sum covariance matrices
     1002    for (int i = 0; i < num; i++) {
     1003        if (inputMask->data.U8[i]) {
     1004            psFree(covariances->data[i]);
     1005            covariances->data[i] = NULL;
     1006        }
     1007    }
     1008    outRO->covariance = psImageCovarianceSum(covariances);
     1009    psFree(covariances);
    9811010
    9821011    if (stats) {
    … …  
    10401069            psString imageResolved = pmConfigConvertFilename(imageNames->data[i], config, false, false);
    10411070            psString maskResolved = pmConfigConvertFilename(maskNames->data[i], config, false, false);
    1042             psString weightResolved = pmConfigConvertFilename(weightNames->data[i], config, false, false);
     1071            psString varianceResolved = pmConfigConvertFilename(varianceNames->data[i], config, false, false);
    10431072            if (unlink(imageResolved) == -1 || unlink(maskResolved) == -1 ||
    1044                 unlink(weightResolved) == -1) {
     1073                unlink(varianceResolved) == -1) {
    10451074                psWarning("Unable to delete temporary files for image %d", i);
    10461075            }
    10471076            psFree(imageResolved);
    10481077            psFree(maskResolved);
    1049             psFree(weightResolved);
     1078            psFree(varianceResolved);
    10501079        }
    10511080    }
    10521081    psFree(imageNames);
    10531082    psFree(maskNames);
    1054     psFree(weightNames);
     1083    psFree(varianceNames);
    10551084
    10561085    psFree(inputMask);
    … …  
    10921121        float renormWidth = psMetadataLookupS32(&mdok, recipe, "RENORM.WIDTH"); // Width of Gaussian phot
    10931122        psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
    1094         if (!pmReadoutWeightRenormPhot(outRO, maskValue, renormNum, renormWidth,
    1095                                        renormMean, renormStdev, NULL)) {
     1123        if (!pmReadoutVarianceRenormPhot(outRO, maskValue, renormNum, renormWidth,
     1124                                         renormMean, renormStdev, NULL)) {
    10961125            psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances.");
    10971126            psFree(outRO);

  • trunk/ppStack/src/ppStackMatch.c

    r21199 r21366  
    1616                     PM_SOURCE_MODE_CR_LIMIT) // Mask to apply to input sources
    1717#define FAINT_SOURCE_FRAC 1.0e-4         // Set minimum flux to this fraction of faintest source flux
     18#define COVAR_FRAC 0.01                 // Truncation fraction for covariance matrix
    1819
    1920//#define TESTING                         // Enable debugging output
    … …  
    233234        }
    234235
    235         // Read image, mask, weight
     236        // Read image, mask, variance
    236237        const char *tempImage = psMetadataLookupStr(NULL, recipe, "TEMP.IMAGE"); // Suffix for image
    237238        const char *tempMask = psMetadataLookupStr(NULL, recipe, "TEMP.MASK"); // Suffix for mask
    238         const char *tempWeight = psMetadataLookupStr(NULL, recipe, "TEMP.WEIGHT"); // Suffix for weight map
    239         psString imageName = NULL, maskName = NULL, weightName = NULL; // Names for convolved images
     239        const char *tempVariance = psMetadataLookupStr(NULL, recipe, "TEMP.VARIANCE"); // Suffix for variance map
     240        psString imageName = NULL, maskName = NULL, varianceName = NULL; // Names for convolved images
    240241        psStringAppend(&imageName, "%s.%d.%s", outName, numInput, tempImage);
    241242        psStringAppend(&maskName, "%s.%d.%s", outName, numInput, tempMask);
    242         psStringAppend(&weightName, "%s.%d.%s", outName, numInput, tempWeight);
     243        psStringAppend(&varianceName, "%s.%d.%s", outName, numInput, tempVariance);
    243244
    244245        if (!readImage(&readout->image, imageName, config) || !readImage(&readout->mask, maskName, config) ||
    245             !readImage(&readout->weight, weightName, config)) {
     246            !readImage(&readout->variance, varianceName, config)) {
    246247            psError(PS_ERR_IO, false, "Unable to read previously produced image.");
    247248            psFree(imageName);
    248249            psFree(maskName);
    249             psFree(weightName);
     250            psFree(varianceName);
    250251            return false;
    251252        }
    252253        psFree(imageName);
    253254        psFree(maskName);
    254         psFree(weightName);
     255        psFree(varianceName);
    255256
    256257        psRegion *region = psMetadataLookupPtr(NULL, output->analysis,
    … …  
    374375                                                                                         "RENORM.STDEV"));
    375376
    376                 if (!pmReadoutWeightRenormPixels(readout, maskBad, renormMean, renormStdev, rng)) {
     377                if (!pmReadoutVarianceRenormPixels(readout, maskBad, renormMean, renormStdev, rng)) {
    377378                    psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances.");
    378379                    psFree(output);
    … …  
    444445            psFree(readout->image);
    445446            psFree(readout->mask);
    446             psFree(readout->weight);
     447            psFree(readout->variance);
     448            psFree(readout->covariance);
    447449            readout->image  = psMemIncrRefCounter(output->image);
    448450            readout->mask   = psMemIncrRefCounter(output->mask);
    449             readout->weight = psMemIncrRefCounter(output->weight);
     451            readout->variance = psMemIncrRefCounter(output->variance);
     452            readout->covariance = psImageCovarianceTruncate(output->covariance, COVAR_FRAC);
    450453        } else {
    451454            // Fake the convolution
    … …  
    559562                                                                                 "RENORM.STDEV"));
    560563
    561         if (!pmReadoutWeightRenormPixels(readout, maskBad, renormMean, renormStdev, rng)) {
     564        if (!pmReadoutVarianceRenormPixels(readout, maskBad, renormMean, renormStdev, rng)) {
    562565            psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances.");
    563566            psFree(output);
    … …  
    583586    psFree(bg);
    584587
    585 
    586588#if 0
    587589#define RADIUS 10                       // Radius of photometry

  • trunk/ppStack/src/ppStackThread.c

    r20711 r21366  
    3838        psFitsClose(stack->imageFits->data[i]);
    3939        psFitsClose(stack->maskFits->data[i]);
    40         psFitsClose(stack->weightFits->data[i]);
    41         stack->imageFits->data[i] = stack->maskFits->data[i] = stack->weightFits->data[i] = NULL;
     40        psFitsClose(stack->varianceFits->data[i]);
     41        stack->imageFits->data[i] = stack->maskFits->data[i] = stack->varianceFits->data[i] = NULL;
    4242    }
    4343    psFree(stack->imageFits);
    4444    psFree(stack->maskFits);
    45     psFree(stack->weightFits);
     45    psFree(stack->varianceFits);
    4646    return;
    4747}
    4848
    4949ppStackThreadData *ppStackThreadDataSetup(const psArray *cells, const psArray *imageNames,
    50                                           const psArray *maskNames, const psArray *weightNames,
     50                                          const psArray *maskNames, const psArray *varianceNames,
    5151                                          const pmConfig *config)
    5252{
    … …  
    5454    PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, imageNames, NULL);
    5555    PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, maskNames, NULL);
    56     PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, weightNames, NULL);
     56    PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, varianceNames, NULL);
    5757
    5858    ppStackThreadData *stack = psAlloc(sizeof(ppStackThreadData)); // Thread data, to return
    … …  
    6464    stack->imageFits  = psArrayAlloc(numInputs);
    6565    stack->maskFits   = psArrayAlloc(numInputs);
    66     stack->weightFits = psArrayAlloc(numInputs);
     66    stack->varianceFits = psArrayAlloc(numInputs);
    6767    for (int i = 0; i < numInputs; i++) {
    6868        if (!cells->data[i]) {
    … …  
    7373        psString imageResolved = pmConfigConvertFilename(imageNames->data[i], config, false, false);
    7474        psString maskResolved = pmConfigConvertFilename(maskNames->data[i], config, false, false);
    75         psString weightResolved = pmConfigConvertFilename(weightNames->data[i], config, false, false);
     75        psString varianceResolved = pmConfigConvertFilename(varianceNames->data[i], config, false, false);
    7676        stack->imageFits->data[i] = psFitsOpen(imageResolved, "r");
    7777        stack->maskFits->data[i] = psFitsOpen(maskResolved, "r");
    78         stack->weightFits->data[i] = psFitsOpen(weightResolved, "r");
     78        stack->varianceFits->data[i] = psFitsOpen(varianceResolved, "r");
    7979        psFree(imageResolved);
    8080        psFree(maskResolved);
    81         psFree(weightResolved);
    82         if (!stack->imageFits->data[i] || !stack->maskFits->data[i] || !stack->weightFits->data[i]) {
     81        psFree(varianceResolved);
     82        if (!stack->imageFits->data[i] || !stack->maskFits->data[i] || !stack->varianceFits->data[i]) {
    8383            psError(PS_ERR_UNKNOWN, false, "Unable to open convolved files %s, %s, %s",
    84                     (char*)imageNames->data[i], (char*)maskNames->data[i], (char*)weightNames->data[i]);
     84                    (char*)imageNames->data[i], (char*)maskNames->data[i], (char*)varianceNames->data[i]);
    8585            return false;
    8686        }
    … …  
    156156                // override the recorded last scan
    157157                ro->thisImageScan  = thread->firstScan;
    158                 ro->thisWeightScan = thread->firstScan;
     158                ro->thisVarianceScan = thread->firstScan;
    159159                ro->thisMaskScan   = thread->firstScan;
    160160                ro->lastImageScan  = thread->lastScan;
    161161                ro->lastMaskScan   = thread->lastScan;
    162                 ro->lastWeightScan = thread->lastScan;
     162                ro->lastVarianceScan = thread->lastScan;
    163163                ro->forceScan      = true;
    164164
    165165                psFits *imageFits  = stack->imageFits->data[i]; // FITS file for image
    166166                psFits *maskFits   = stack->maskFits->data[i]; // FITS file for mask
    167                 psFits *weightFits = stack->weightFits->data[i]; // FITS file for weight
     167                psFits *varianceFits = stack->varianceFits->data[i]; // FITS file for variance
    168168
    169169
    … …  
    189189                }
    190190
    191                 if (pmReadoutMoreWeight(ro, weightFits, 0, rows, config)) {
     191                if (pmReadoutMoreVariance(ro, varianceFits, 0, rows, config)) {
    192192                    keepReading = true;
    193                     if (!pmReadoutReadChunkWeight(ro, weightFits, 0, rows, overlap, config)) {
    194                         psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPSTACK.INPUT.WEIGHT %d",
     193                    if (!pmReadoutReadChunkVariance(ro, varianceFits, 0, rows, overlap, config)) {
     194                        psError(PS_ERR_IO, false,
     195                                "Unable to read chunk %d for file PPSTACK.INPUT.VARIANCE %d",
    195196                                numChunk, i);
    196197                        *status = false;

  • trunk/psphot/src/psphot.h

    r21255 r21366  
    9090void            psphotModelClassInit (void);
    9191bool            psphotGrowthCurve (pmReadout *readout, pmPSF *psf, bool ignore, psImageMaskType maskVal);
    92 bool            psphotSetMaskAndWeight (pmConfig *config, pmReadout *readout, psMetadata *recipe);
     92bool            psphotSetMaskAndVariance (pmConfig *config, pmReadout *readout, psMetadata *recipe);
    9393void            psphotSourceFreePixels (psArray *sources);
    9494

  • trunk/psphot/src/psphotAnnuli.c

    r21183 r21366  
    1111
    1212  psVector *radius = source->extpars->profile->radius;
    13   psVector *weight = source->extpars->profile->weight;
     13  psVector *variance = source->extpars->profile->variance;
    1414  psVector *flux = source->extpars->profile->flux;
    1515
    … …  
    2222  psVector *fluxValues = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32);
    2323  psVector *fluxSquare = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32);
    24   psVector *fluxWeight = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32);
     24  psVector *fluxVariance = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32);
    2525  psVector *pixelCount = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32);
    2626  psVectorInit (fluxValues, 0.0);
    2727  psVectorInit (fluxSquare, 0.0);
    28   psVectorInit (fluxWeight, 0.0);
     28  psVectorInit (fluxVariance, 0.0);
    2929  psVectorInit (pixelCount, 0.0);
    3030
    3131  // XXX this code assumes the radii are in pixels.  convert from arcsec with plate scale
    3232  // XXX assume the annulii above are not overlapping?  much faster...
    33   // XXX this might be must faster in the reverse order: loop over annulii and use disection to 
     33  // XXX this might be must faster in the reverse order: loop over annulii and use disection to
    3434  // skip to the start of the annulus.
    3535  for (int i = 0; i < flux->n; i++) {
    … …  
    3939      fluxValues->data.F32[j] += flux->data.F32[i];
    4040      fluxSquare->data.F32[j] += PS_SQR(flux->data.F32[i]);
    41       fluxWeight->data.F32[j] += weight->data.F32[i];
     41      fluxVariance->data.F32[j] += variance->data.F32[i];
    4242      pixelCount->data.F32[j] += 1.0;
    4343    }
    … …  
    4949    fluxSquare->data.F32[j] -= PS_SQR(fluxValues->data.F32[j]);
    5050  }
    51  
     51
    5252  source->extpars->annuli = pmSourceAnnuliAlloc ();
    5353  source->extpars->annuli->flux    = fluxValues;
    54   source->extpars->annuli->fluxErr = fluxWeight;
     54  source->extpars->annuli->fluxErr = fluxVariance;
    5555  source->extpars->annuli->fluxVar = fluxSquare;
    5656

  • trunk/psphot/src/psphotBlendFit.c

    r21359 r21366  
    2121    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
    2222    if (!status) {
    23         nThreads = 0;
     23        nThreads = 0;
    2424    }
    2525
    … …  
    4444    psphotInitRadiusPSF (recipe, psf->type);
    4545
    46     // starts the timer, sets up the array of fitSets 
     46    // starts the timer, sets up the array of fitSets
    4747    psphotFitInit (nThreads);
    4848
    … …  
    5656    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
    5757
    58     fprintf (stderr, "starting with %ld sources\n", sources->n);
    59 
    6058    for (int i = 0; i < cellGroups->n; i++) {
    6159
    62         psArray *cells = cellGroups->data[i];
    63 
    64         for (int j = 0; j < cells->n; j++) {
    65 
    66             if (nThreads) {
    67                 // allocate a job -- if threads are not defined, this just runs the job
    68                 psThreadJob *job = psThreadJobAlloc ("PSPHOT_BLEND_FIT");
    69                 psArray *newSources = psArrayAllocEmpty(16);
    70 
    71                 psArrayAdd(job->args, 1, readout);
    72                 psArrayAdd(job->args, 1, recipe);
    73                 psArrayAdd(job->args, 1, cells->data[j]); // sources
    74                 psArrayAdd(job->args, 1, psf);
    75                 psArrayAdd(job->args, 1, newSources); // return for new sources
    76                 psFree (newSources);
    77 
    78                 PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfit
    79                 PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npsf
    80                 PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
    81                 PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
    82 
    83                 if (!psThreadJobAddPending(job)) {
    84                     psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    85                     psFree (job);
    86                     return NULL;
    87                 }
    88                 psFree(job);
    89             } else {
    90                 int nfit = 0;
    91                 int npsf = 0;
    92                 int next = 0;
    93                 int nfail = 0;
    94                 psArray *newSources = psArrayAllocEmpty(16);
    95 
    96                 if (!psphotBlendFit_Unthreaded (&nfit, &npsf, &next, &nfail, readout, recipe, cells->data[j], psf, newSources)) {
    97                     psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    98                     return NULL;
    99                 }
    100                 Nfit += nfit;
    101                 Npsf += npsf;
    102                 Next += next;
    103                 Nfail += nfail;
    104 
    105                 // add these back onto sources
    106                 for (int k = 0; k < newSources->n; k++) {
    107                     psArrayAdd (sources, 16, newSources->data[k]);
    108                 }
    109                 psFree (newSources);
    110             }
    111         }
    112 
    113         if (nThreads) {
    114             // wait for the threads to finish and manage results
    115             if (!psThreadPoolWait (false)) {
    116                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    117                 return NULL;
    118             }
    119 
    120             // we have only supplied one type of job, so we can assume the types here
    121             psThreadJob *job = NULL;
    122             while ((job = psThreadJobGetDone()) != NULL) {
    123                 if (job->args->n < 1) {
    124                     fprintf (stderr, "error with job\n");
    125                 } else {
    126                     psScalar *scalar = NULL;
    127                     scalar = job->args->data[5];
    128                     Nfit += scalar->data.S32;
    129                     scalar = job->args->data[6];
    130                     Npsf += scalar->data.S32;
    131                     scalar = job->args->data[7];
    132                     Next += scalar->data.S32;
    133                     scalar = job->args->data[8];
    134                     Nfail += scalar->data.S32;
    135 
    136                     // add these back onto sources
    137                     psArray *newSources = job->args->data[4];
    138                     for (int j = 0; j < newSources->n; j++) {
    139                         psArrayAdd (sources, 16, newSources->data[j]);
    140                     }
    141                 }
    142                 psFree(job);
    143             }
    144         }
     60        psArray *cells = cellGroups->data[i];
     61
     62        for (int j = 0; j < cells->n; j++) {
     63
     64            if (nThreads) {
     65                // allocate a job -- if threads are not defined, this just runs the job
     66                psThreadJob *job = psThreadJobAlloc ("PSPHOT_BLEND_FIT");
     67                psArray *newSources = psArrayAllocEmpty(16);
     68
     69                psArrayAdd(job->args, 1, readout);
     70                psArrayAdd(job->args, 1, recipe);
     71                psArrayAdd(job->args, 1, cells->data[j]); // sources
     72                psArrayAdd(job->args, 1, psf);
     73                psArrayAdd(job->args, 1, newSources); // return for new sources
     74                psFree (newSources);
     75
     76                PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfit
     77                PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npsf
     78                PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
     79                PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
     80
     81                if (!psThreadJobAddPending(job)) {
     82                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     83                    psFree (job);
     84                    return NULL;
     85                }
     86                psFree(job);
     87            } else {
     88                int nfit = 0;
     89                int npsf = 0;
     90                int next = 0;
     91                int nfail = 0;
     92                psArray *newSources = psArrayAllocEmpty(16);
     93
     94                if (!psphotBlendFit_Unthreaded (&nfit, &npsf, &next, &nfail, readout, recipe, cells->data[j], psf, newSources)) {
     95                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     96                    return NULL;
     97                }
     98                Nfit += nfit;
     99                Npsf += npsf;
     100                Next += next;
     101                Nfail += nfail;
     102
     103                // add these back onto sources
     104                for (int k = 0; k < newSources->n; k++) {
     105                    psArrayAdd (sources, 16, newSources->data[k]);
     106                }
     107                psFree (newSources);
     108            }
     109        }
     110
     111        if (nThreads) {
     112            // wait for the threads to finish and manage results
     113            if (!psThreadPoolWait (false)) {
     114                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     115                return NULL;
     116            }
     117
     118            // we have only supplied one type of job, so we can assume the types here
     119            psThreadJob *job = NULL;
     120            while ((job = psThreadJobGetDone()) != NULL) {
     121                if (job->args->n < 1) {
     122                    fprintf (stderr, "error with job\n");
     123                } else {
     124                    psScalar *scalar = NULL;
     125                    scalar = job->args->data[5];
     126                    Nfit += scalar->data.S32;
     127                    scalar = job->args->data[6];
     128                    Npsf += scalar->data.S32;
     129                    scalar = job->args->data[7];
     130                    Next += scalar->data.S32;
     131                    scalar = job->args->data[8];
     132                    Nfail += scalar->data.S32;
     133
     134                    // add these back onto sources
     135                    psArray *newSources = job->args->data[4];
     136                    for (int j = 0; j < newSources->n; j++) {
     137                        psArrayAdd (sources, 16, newSources->data[j]);
     138                    }
     139                }
     140                psFree(job);
     141            }
     142        }
    145143    }
    146144    psFree (cellGroups);
    … …  
    207205        if (source->peak->SN < FIT_SN_LIM) continue;
    208206
    209         // exclude sources outside optional analysis region
     207        // exclude sources outside optional analysis region
    210208        if (source->peak->xf < AnalysisRegion.x0) continue;
    211209        if (source->peak->yf < AnalysisRegion.y0) continue;
    … …  
    216214        if (source->modelPSF == NULL) continue;
    217215
    218         // skip sources which are insignificant flux? 
    219         // XXX this is somewhat ad-hoc
     216        // skip sources which are insignificant flux?
     217        // XXX this is somewhat ad-hoc
    220218        if (source->modelPSF->params->data.F32[1] < 0.1) {
    221219            psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n",
    … …  
    234232        // try fitting PSFs or extended sources depending on source->mode
    235233        // these functions subtract the resulting fitted source
    236         if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
    237             if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) {
    238                 source->type = PM_SOURCE_TYPE_EXTENDED;
    239                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
    240                 Next ++;
    241                 continue;
    242             }
    243         } else {
    244             if (psphotFitBlend (readout, source, psf, maskVal, markVal)) {
    245                 source->type = PM_SOURCE_TYPE_STAR;
    246                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
    247                 Npsf ++;
    248                 continue;
    249             }
    250         }
     234        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
     235            if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) {
     236                source->type = PM_SOURCE_TYPE_EXTENDED;
     237                psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
     238                Next ++;
     239                continue;
     240            }
     241        } else {
     242            if (psphotFitBlend (readout, source, psf, maskVal, markVal)) {
     243                source->type = PM_SOURCE_TYPE_STAR;
     244                psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
     245                Npsf ++;
     246                continue;
     247            }
     248        }
    251249
    252250        psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->peak->xf, source->peak->yf);
    … …  
    317315        if (source->peak->SN < FIT_SN_LIM) continue;
    318316
    319         // exclude sources outside optional analysis region
     317        // exclude sources outside optional analysis region
    320318        if (source->peak->xf < AnalysisRegion.x0) continue;
    321319        if (source->peak->yf < AnalysisRegion.y0) continue;
    … …  
    326324        if (source->modelPSF == NULL) continue;
    327325
    328         // skip sources which are insignificant flux? 
    329         // XXX this is somewhat ad-hoc
     326        // skip sources which are insignificant flux?
     327        // XXX this is somewhat ad-hoc
    330328        if (source->modelPSF->params->data.F32[1] < 0.1) {
    331329            psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n",
    … …  
    344342        // try fitting PSFs or extended sources depending on source->mode
    345343        // these functions subtract the resulting fitted source
    346         if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
    347             if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) {
    348                 source->type = PM_SOURCE_TYPE_EXTENDED;
    349                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
    350                 Next ++;
    351                 continue;
    352             }
    353         } else {
    354             if (psphotFitBlend (readout, source, psf, maskVal, markVal)) {
    355                 source->type = PM_SOURCE_TYPE_STAR;
    356                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
    357                 Npsf ++;
    358                 continue;
    359             }
    360         }
     344        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
     345            if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) {
     346                source->type = PM_SOURCE_TYPE_EXTENDED;
     347                psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
     348                Next ++;
     349                continue;
     350            }
     351        } else {
     352            if (psphotFitBlend (readout, source, psf, maskVal, markVal)) {
     353                source->type = PM_SOURCE_TYPE_STAR;
     354                psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
     355                Npsf ++;
     356                continue;
     357            }
     358        }
    361359
    362360        psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->peak->xf, source->peak->yf);

  • trunk/psphot/src/psphotFindFootprints.c

    r21183 r21366  
    6868    }
    6969
    70     psphotCullPeaks(readout->image, readout->weight, recipe, detections->footprints);
     70    psphotCullPeaks(readout->image, readout->variance, recipe, detections->footprints);
    7171    detections->peaks = pmFootprintArrayToPeaks(detections->footprints);
    7272    psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks, %ld total footprints: %f sec\n", detections->peaks->n, detections->footprints->n, psTimerMark ("psphot.footprints"));

  • trunk/psphot/src/psphotFitSourcesLinear.c

    r21183 r21366  
    6868
    6969        // if (source->type == PM_SOURCE_TYPE_STAR &&
    70         if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
     70        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
    7171
    7272        if (final) {
    … …  
    7676        }
    7777
    78         // generate model for sources without, or skip if we can't
    79         if (!source->modelFlux) {
    80             if (!pmSourceCacheModel (source, maskVal)) continue;
    81         }
     78        // generate model for sources without, or skip if we can't
     79        if (!source->modelFlux) {
     80            if (!pmSourceCacheModel (source, maskVal)) continue;
     81        }
    8282
    8383        // save the original coords
    … …  
    9696
    9797    // vectors to store stats for each object
    98     // psVector *weight = psVectorAlloc (fitSources->n, PS_TYPE_F32);
     98    // psVector *variance = psVectorAlloc (fitSources->n, PS_TYPE_F32);
    9999    psVector *errors = psVectorAlloc (fitSources->n, PS_TYPE_F32);
    100100
    … …  
    128128        psSparseVectorElement (sparse, i, f);
    129129
    130         // add the per-source weights (border region)
     130        // add the per-source variances (border region)
    131131        switch (SKY_FIT_ORDER) {
    132132          case 1:
    … …  
    208208        pmSource *source = fitSources->data[i];
    209209        pmModel *model = pmSourceGetModel (NULL, source);
    210         pmSourceChisq (model, source->pixels, source->maskObj, source->weight, maskVal);
     210        pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal);
    211211    }
    212212    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
    … …  
    253253    // accumulate the image statistics from the masked regions
    254254    psF32 **image  = readout->image->data.F32;
    255     psF32 **weight = readout->weight->data.F32;
     255    psF32 **variance = readout->variance->data.F32;
    256256    psImageMaskType  **mask   = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA;
    257257
    … …  
    268268                wt = 1.0;
    269269            } else {
    270                 wt = weight[j][i];
     270                wt = variance[j][i];
    271271            }
    272272            f = image[j][i];

  • trunk/psphot/src/psphotMakeResiduals.c

    r21183 r21366  
    5858    // for each input source:
    5959    // - construct a residual image, renormalized
    60     // - construct a renormalized weight image
     60    // - construct a renormalized variance image
    6161    // - construct a new mask image
    6262
    6363    // construct the output residual table (Nx*DX,Ny*DY)
    6464    // for each output pixel:
    65     // - construct a histogram of the values & weights (interpolate to the common pixel coordinate)
     65    // - construct a histogram of the values & variances (interpolate to the common pixel coordinate)
    6666    // - measure the robust median & sigma
    6767    // - reject (mask) input pixels which are outliers
    6868    // - re-measure the robust median & sigma
    69     // - set output pixel, weight, and mask
     69    // - set output pixel, variance, and mask
    7070
    7171    // these mask values do not correspond to the recipe values: they
    … …  
    7676    // psVectorStats
    7777
    78     const psImageMaskType badMask     = 0x01;   // mask bits
    79     const psImageMaskType poorMask    = 0x02;   // from psImageInterpolate
    80     const psImageMaskType clippedMask = 0x04;   // mask bit set for clipped values
     78    const psImageMaskType badMask     = 0x01;   // mask bits
     79    const psImageMaskType poorMask    = 0x02;   // from psImageInterpolate
     80    const psImageMaskType clippedMask = 0x04;   // mask bit set for clipped values
    8181    const psVectorMaskType fmaskVal = badMask | poorMask | clippedMask;
    8282
    … …  
    101101
    102102        psImage *image  = psImageCopy (NULL, source->pixels,   PS_TYPE_F32);
    103         psImage *weight = psImageCopy (NULL, source->weight,   PS_TYPE_F32);
     103        psImage *variance = psImageCopy (NULL, source->variance,   PS_TYPE_F32);
    104104        psImage *mask   = psImageCopy (NULL, source->maskView, PS_TYPE_IMAGE_MASK);
    105105        pmModelSub (image, mask, model, PM_MODEL_OP_FUNC, maskVal);
    106106
    107         // re-normalize image and weight
     107        // re-normalize image and variance
    108108        float Io = model->params->data.F32[PM_PAR_I0];
    109109        psBinaryOp (image, image, "/", psScalarAlloc(Io, PS_TYPE_F32));
    110         psBinaryOp (weight, weight, "/", psScalarAlloc(Io*Io, PS_TYPE_F32));
    111 
    112         // we interpolate the image and weight - include the mask or not?
    113         // XXX why not the mask?
    114         // psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, weight, mask, maskVal, 0.0, 0.0, badMask, poorMask, 0.0, 0);
    115         psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, weight, NULL, 0xff, 0.0, 0.0, badMask, poorMask, 0.0, 0);
     110        psBinaryOp (variance, variance, "/", psScalarAlloc(Io*Io, PS_TYPE_F32));
     111
     112        // we interpolate the image and variance - include the mask or not?
     113        // XXX why not the mask?
     114        // psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, variance, mask, maskVal, 0.0, 0.0, badMask, poorMask, 0.0, 0);
     115        psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, variance, NULL, 0xff, 0.0, 0.0, badMask, poorMask, 0.0, 0);
    116116        psArrayAdd (input, 100, interp);
    117117
    … …  
    128128        psFree (mask);
    129129        psFree (image);
    130         psFree (weight);
     130        psFree (variance);
    131131        psFree (interp);
    132132    }
    … …  
    176176                    fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = badMask;
    177177                } else {
    178                   fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = mflux; // XXX is mflux IMAGE or VECTOR type?
    179                 }
    180                 fluxes->data.F32[i] = flux;
    181                 dfluxes->data.F32[i] = dflux;
     178                  fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = mflux; // XXX is mflux IMAGE or VECTOR type?
     179                }
     180                fluxes->data.F32[i] = flux;
     181                dfluxes->data.F32[i] = dflux;
    182182                if (isnan(flux)) {
    183183                    fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = badMask;
    … …  
    204204
    205205            // mark input pixels which are more than N sigma from the median
    206             int nKeep = 0;
     206            int nKeep = 0;
    207207            for (int i = 0; i < fluxes->n; i++) {
    208208                float delta = fluxes->data.F32[i] - fluxClip->robustMedian;
    … …  
    210210                float swing = fabs(delta) / sigma;
    211211
    212                 // mask pixels which are out of range
     212                // mask pixels which are out of range
    213213                if (swing > nSigma) {
    214214                    fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = clippedMask;
    215215                }
    216                 if (!fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) nKeep++;
     216                if (!fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) nKeep++;
    217217            }
    218218
    … …  
    225225                resid->Ro->data.F32[oy][ox] = psStatsGetValue(fluxStats, statOption);
    226226                resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0;
    227                 //resid->weight->data.F32[oy][ox] = fluxStats->sampleStdev;
     227                //resid->variance->data.F32[oy][ox] = fluxStats->sampleStdev;
    228228
    229229                if (fabs(resid->Ro->data.F32[oy][ox]) < pixelSN*fluxStats->sampleStdev/sqrt(nKeep)) {
    … …  
    231231                }
    232232
    233                 // fprintf (stderr, "res: %2d %2d : %6.4f  %6.4f  %6.4f   %3d  %1d\n", ox, oy, resid->Ro->data.F32[oy][ox], fluxStats->sampleStdev, fluxStats->sampleStdev/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]);
     233                // fprintf (stderr, "res: %2d %2d : %6.4f  %6.4f  %6.4f   %3d  %1d\n", ox, oy, resid->Ro->data.F32[oy][ox], fluxStats->sampleStdev, fluxStats->sampleStdev/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]);
    234234
    235235            } else {
    … …  
    268268
    269269                float dRo = sqrt(A->data.F32[0][0]);
    270                 // fprintf (stderr, "res: %2d %2d : %6.4f  %6.4f  %6.4f   %3d  %1d\n",
    271                 // ox, oy, resid->Ro->data.F32[oy][ox], dRo, dRo/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]);
     270                // fprintf (stderr, "res: %2d %2d : %6.4f  %6.4f  %6.4f   %3d  %1d\n",
     271                // ox, oy, resid->Ro->data.F32[oy][ox], dRo, dRo/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]);
    272272
    273273                if (fabs(resid->Ro->data.F32[oy][ox]) < pixelSN*dRo/sqrt(nKeep)) {
    274274                  resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
    275275                }
    276                 //resid->weight->data.F32[oy][ox] = XXX;
     276                //resid->variance->data.F32[oy][ox] = XXX;
    277277            }
    278278        }
    … …  
    301301      psphotSaveImage (NULL, resid->Rx,     "resid.rx.fits");
    302302      psphotSaveImage (NULL, resid->Ry,     "resid.ry.fits");
    303       psphotSaveImage (NULL, resid->weight, "resid.wt.fits");
     303      psphotSaveImage (NULL, resid->variance, "resid.wt.fits");
    304304      psphotSaveImage (NULL, resid->mask,   "resid.mk.fits");
    305305    }

  • trunk/psphot/src/psphotMaskReadout.c

    r21183 r21366  
    11# include "psphotInternal.h"
    22
    3 // generate mask and weight if not defined, additional mask for restricted subregion
    4 bool psphotSetMaskAndWeight (pmConfig *config, pmReadout *readout, psMetadata *recipe) {
     3// generate mask and variance if not defined, additional mask for restricted subregion
     4bool psphotSetMaskAndVariance (pmConfig *config, pmReadout *readout, psMetadata *recipe) {
    55
    66    bool status;
    … …  
    1515    psMetadataAddImageMask (recipe, PS_LIST_TAIL, "MASK.BAD", PS_META_REPLACE, "user-defined mask", maskBad);
    1616
    17     // generate mask & weight images if they don't already exit
     17    // generate mask & variance images if they don't already exit
    1818    if (!readout->mask) {
    1919        if (!pmReadoutGenerateMask(readout, maskSat, maskBad)) {
    … …  
    2222        }
    2323    }
    24     if (!readout->weight) {
    25         if (!pmReadoutGenerateWeight(readout, true)) {
    26             psError (PSPHOT_ERR_CONFIG, false, "trouble creating weight");
     24    if (!readout->variance) {
     25        if (!pmReadoutGenerateVariance(readout, true)) {
     26            psError (PSPHOT_ERR_CONFIG, false, "trouble creating variance");
    2727            return false;
    2828        }
    … …  
    4949        psphotSaveImage (NULL, readout->image,  "image.fits");
    5050        psphotSaveImage (NULL, readout->mask,   "mask.fits");
    51         psphotSaveImage (NULL, readout->weight, "weight.fits");
     51        psphotSaveImage (NULL, readout->variance, "variance.fits");
    5252    }
    5353

  • trunk/psphot/src/psphotModelWithPSF.c

    r21183 r21366  
    2020        paramMask = constraint->paramMask;
    2121        if (paramMask != NULL) {
    22           PS_ASSERT_VECTOR_TYPE(paramMask, PS_TYPE_VECTOR_MASK, false);
     22          PS_ASSERT_VECTOR_TYPE(paramMask, PS_TYPE_VECTOR_MASK, false);
    2323            PS_ASSERT_VECTORS_SIZE_EQUAL(params, paramMask, false);
    2424        }
    … …  
    4242    // Alpha & Beta only contain elements to represent the unmasked parameters
    4343    if (!psMinLM_AllocAB (&Alpha, &Beta, params, paramMask)) {
    44         psAbort ("programming error: no unmasked parameters to be fit\n");
    45     }
    46    
     44        psAbort ("programming error: no unmasked parameters to be fit\n");
     45    }
     46
    4747    // allocate internal arrays (current vs Guess)
    4848    psImage *alpha   = psImageAlloc(Alpha->numCols, Alpha->numRows, PS_TYPE_F32);
    … …  
    127127            lambda *= 0.25;
    128128
    129             // save the new convolved model image
    130             psFree (source->modelFlux);
    131             source->modelFlux = pmPCMDataSaveImage(pcm);
     129            // save the new convolved model image
     130            psFree (source->modelFlux);
     131            source->modelFlux = pmPCMDataSaveImage(pcm);
    132132        } else {
    133133            lambda *= 10.0;
    … …  
    142142            psTrace ("psphot", 5, "failure to calculate covariance matrix\n");
    143143        }
    144         // set covar values which are not masked
    145         psImageInit (covar, 0.0);
    146         for (int j = 0, J = 0; j < params->n; j++) {
    147             if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[j])) {
    148                 covar->data.F32[j][j] = 1.0;
    149                 continue;
    150             }
    151             for (int k = 0, K = 0; k < params->n; k++) {
    152                 if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[k])) continue;
    153                 covar->data.F32[j][k] = Alpha->data.F32[J][K];
    154                 K++;
    155             }
    156             J++;
    157         }
     144        // set covar values which are not masked
     145        psImageInit (covar, 0.0);
     146        for (int j = 0, J = 0; j < params->n; j++) {
     147            if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[j])) {
     148                covar->data.F32[j][j] = 1.0;
     149                continue;
     150            }
     151            for (int k = 0, K = 0; k < params->n; k++) {
     152                if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[k])) continue;
     153                covar->data.F32[j][k] = Alpha->data.F32[J][K];
     154                K++;
     155            }
     156            J++;
     157        }
    158158    }
    159159
    … …  
    192192    PS_ASSERT_PTR_NON_NULL(source, NAN);
    193193    PS_ASSERT_IMAGE_NON_NULL(source->pixels, NAN);
    194     PS_ASSERT_IMAGE_NON_NULL(source->weight, NAN);
     194    PS_ASSERT_IMAGE_NON_NULL(source->variance, NAN);
    195195    PS_ASSERT_IMAGE_NON_NULL(source->maskObj, NAN);
    196196
    … …  
    210210    psImageInit (pcm->model, 0.0);
    211211    for (int n = 0; n < params->n; n++) {
    212         if (!pcm->dmodels->data[n]) continue;
    213         psImageInit (pcm->dmodels->data[n], 0.0);
     212        if (!pcm->dmodels->data[n]) continue;
     213        psImageInit (pcm->dmodels->data[n], 0.0);
    214214    }
    215215
    … …  
    218218        for (psS32 j = 0; j < source->pixels->numCols; j++) {
    219219
    220             // XXX can we skip some of the data points where the model
    221             // is not going to be fitted??
     220            // XXX can we skip some of the data points where the model
     221            // is not going to be fitted??
    222222
    223223            // skip masked points
    224             // XXX probably should not skipped masked points:
    225             // XXX skip if convolution of unmasked pixels will not see this pixel
     224            // XXX probably should not skipped masked points:
     225            // XXX skip if convolution of unmasked pixels will not see this pixel
    226226            // if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[i][j]) {
    227             // continue;
    228             // }
    229 
    230             // skip zero-weight points
    231             // XXX why is this not masked?
    232             // if (source->weight->data.F32[i][j] == 0) {
    233             // continue;
    234             // }
     227            // continue;
     228            // }
     229
     230            // skip zero-variance points
     231            // XXX why is this not masked?
     232            // if (source->variance->data.F32[i][j] == 0) {
     233            // continue;
     234            // }
    235235            // skip nan value points
    236             // XXX why is this not masked?
     236            // XXX why is this not masked?
    237237            // if (!isfinite(source->pixels->data.F32[i][j])) {
    238             // continue;
    239             // }
     238            // continue;
     239            // }
    240240
    241241            // Convert i/j to image space:
    … …  
    243243            coord->data.F32[1] = (psF32) (i + source->pixels->row0);
    244244
    245             pcm->model->data.F32[i][j] = func (deriv, params, coord);
    246 
    247             for (int n = 0; n < params->n; n++) {
    248                 if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; }
    249                 psImage *dmodel = pcm->dmodels->data[n];
    250                 dmodel->data.F32[i][j] = deriv->data.F32[n];
    251             }
     245            pcm->model->data.F32[i][j] = func (deriv, params, coord);
     246
     247            for (int n = 0; n < params->n; n++) {
     248                if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; }
     249                psImage *dmodel = pcm->dmodels->data[n];
     250                dmodel->data.F32[i][j] = deriv->data.F32[n];
     251            }
    252252        }
    253253    }
    … …  
    258258    psImageConvolveDirect (pcm->modelConv, pcm->model, psf);
    259259    for (int n = 0; n < pcm->dmodels->n; n++) {
    260         if (pcm->dmodels->data[n] == NULL) continue;
    261         psImage *dmodel = pcm->dmodels->data[n];
    262         psImage *dmodelConv = pcm->dmodelsConv->data[n];
    263         psImageConvolveDirect (dmodelConv, dmodel, psf);
     260        if (pcm->dmodels->data[n] == NULL) continue;
     261        psImage *dmodel = pcm->dmodels->data[n];
     262        psImage *dmodelConv = pcm->dmodelsConv->data[n];
     263        psImageConvolveDirect (dmodelConv, dmodel, psf);
    264264    }
    265265
    266266    // XXX TEST : SAVE IMAGES
    267 # if (SAVE_IMAGES) 
     267# if (SAVE_IMAGES)
    268268    psphotSaveImage (NULL, psf->image, "psf.fits");
    269269    psphotSaveImage (NULL, pcm->model, "model.fits");
    … …  
    271271    psphotSaveImage (NULL, source->pixels, "obj.fits");
    272272    psphotSaveImage (NULL, source->maskObj, "mask.fits");
    273     psphotSaveImage (NULL, source->weight, "weight.fits");
     273    psphotSaveImage (NULL, source->variance, "variance.fits");
    274274# endif
    275275
    276     // 2 *** accumulate alpha & beta 
     276    // 2 *** accumulate alpha & beta
    277277
    278278    // zero alpha and beta for summing below
    … …  
    283283    for (psS32 i = 0; i < source->pixels->numRows; i++) {
    284284        for (psS32 j = 0; j < source->pixels->numCols; j++) {
    285             // XXX are we doing the right thing with the mask?
     285            // XXX are we doing the right thing with the mask?
    286286            // skip masked points
    287287            if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[i][j]) {
    288288                continue;
    289289            }
    290             // skip zero-weight points
    291             if (source->weight->data.F32[i][j] == 0) {
     290            // skip zero-variance points
     291            if (source->variance->data.F32[i][j] == 0) {
    292292                continue;
    293293            }
    … …  
    297297            }
    298298
    299             float ymodel  = pcm->modelConv->data.F32[i][j];
    300             float yweight = 1.0 / source->weight->data.F32[i][j];
    301             float delta = ymodel - source->pixels->data.F32[i][j];
    302 
    303             chisq += PS_SQR(delta) * yweight;
    304 
    305             if (isnan(delta)) psAbort("nan in delta");
    306             if (isnan(chisq)) psAbort("nan in chisq");
    307 
    308             // alpha & beta only contain unmasked elements
    309             for (int n1 = 0, N1 = 0; n1 < params->n; n1++) {
    310                 if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n1])) continue;
    311                 psImage *dmodel = pcm->dmodelsConv->data[n1];
    312                 float weight = dmodel->data.F32[i][j] * yweight;
    313                 for (int n2 = 0, N2 = 0; n2 <= n1; n2++) {
    314                     if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n2])) continue;
    315                     dmodel = pcm->dmodelsConv->data[n2];
    316                     alpha->data.F32[N1][N2] += weight * dmodel->data.F32[i][j];
    317                     N2++;
    318                 }
    319                 beta->data.F32[N1] += weight * delta;
    320                 N1++;
    321             }
    322         }
     299            float ymodel  = pcm->modelConv->data.F32[i][j];
     300            float yweight = 1.0 / source->variance->data.F32[i][j];
     301            float delta = ymodel - source->pixels->data.F32[i][j];
     302
     303            chisq += PS_SQR(delta) * yweight;
     304
     305            if (isnan(delta)) psAbort("nan in delta");
     306            if (isnan(chisq)) psAbort("nan in chisq");
     307
     308            // alpha & beta only contain unmasked elements
     309            for (int n1 = 0, N1 = 0; n1 < params->n; n1++) {
     310                if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n1])) continue;
     311                psImage *dmodel = pcm->dmodelsConv->data[n1];
     312                float weight = dmodel->data.F32[i][j] * yweight;
     313                for (int n2 = 0, N2 = 0; n2 <= n1; n2++) {
     314                    if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n2])) continue;
     315                    dmodel = pcm->dmodelsConv->data[n2];
     316                    alpha->data.F32[N1][N2] += weight * dmodel->data.F32[i][j];
     317                    N2++;
     318                }
     319                beta->data.F32[N1] += weight * delta;
     320                N1++;
     321            }
     322        }
    323323    }
    324324
    … …  
    356356    pcm->dmodels = psArrayAlloc (params->n);
    357357    for (psS32 n = 0; n < params->n; n++) {
    358         pcm->dmodels->data[n] = NULL;
    359         if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; }
    360         pcm->dmodels->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
     358        pcm->dmodels->data[n] = NULL;
     359        if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; }
     360        pcm->dmodels->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
    361361    }
    362362
    … …  
    365365    pcm->dmodelsConv = psArrayAlloc (params->n);
    366366    for (psS32 n = 0; n < params->n; n++) {
    367         pcm->dmodelsConv->data[n] = NULL;
    368         if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; }
    369         pcm->dmodelsConv->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
     367        pcm->dmodelsConv->data[n] = NULL;
     368        if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; }
     369        pcm->dmodelsConv->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
    370370    }
    371371
    … …  
    376376
    377377    psImage *model = psImageCopy (NULL, pcm->modelConv, PS_TYPE_F32);
    378    
     378
    379379    return model;
    380380}
    … …  
    383383 *
    384384 * we have a function func(param; value)
    385  
     385
    386386 * basic LMM:
    387  
     387
    388388 - fill in the data (x, y)
    389  
     389
    390390 chisq = SetABX (alpha, beta, params, paramMask, x, y, dy, func)
    391  
     391
    392392 while () {
    393393 GuessABP (Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda)
    … …  
    396396 convergence tests...
    397397 }
    398  
    399  
     398
     399
    400400
    401401 ** GuessABP:

  • trunk/psphot/src/psphotOutput.c

    r21183 r21366  
    1818
    1919pmReadout *psphotSelectBackgroundStdev (pmConfig *config,
    20                                         const pmFPAview *view) {
     20                                        const pmFPAview *view) {
    2121
    2222    bool status;
    … …  
    7777                continue;
    7878            }
    79             // skip zero-weight points
    80             if (source->weight->data.F32[i][j] == 0) {
     79            // skip zero-variance points
     80            if (source->variance->data.F32[i][j] == 0) {
    8181                continue;
    8282            }
    … …  
    8686                     (i + source->pixels->row0),
    8787                     source->pixels->data.F32[i][j],
    88                      1.0 / source->weight->data.F32[i][j],
     88                     1.0 / source->variance->data.F32[i][j],
    8989                     source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[i][j]);
    9090        }
    … …  
    124124        if (model == NULL)
    125125            continue;
    126         if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
    127             nEXT ++;
    128         }
    129         if (source->mode & PM_SOURCE_MODE_CR_LIMIT) {
    130             nCR ++;
    131         }
     126        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
     127            nEXT ++;
     128        }
     129        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) {
     130            nCR ++;
     131        }
    132132        nSrc ++;
    133133    }
    … …  
    242242
    243243    for (int i = 0; i < try->sources->n; i++) {
    244         // masked for: bad model fit, outlier in parameters
    245         if (try->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL)
    246             continue;
    247 
    248         pmSource *source = try->sources->data[i];
    249         float x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
    250         float y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
    251 
    252         // set the mask and subtract the PSF model
    253         // XXX should we be using maskObj? should we be unsetting the mask?
    254         // use pmModelSub because modelFlux has not been generated
    255         assert (source->maskObj);
    256         psImageKeepCircle (source->maskObj, x, y, radius, "OR", markVal);
    257         pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL, maskVal);
    258         psImageKeepCircle (source->maskObj, x, y, radius, "AND", PS_NOT_IMAGE_MASK(markVal));
     244        // masked for: bad model fit, outlier in parameters
     245        if (try->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL)
     246            continue;
     247
     248        pmSource *source = try->sources->data[i];
     249        float x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
     250        float y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
     251
     252        // set the mask and subtract the PSF model
     253        // XXX should we be using maskObj? should we be unsetting the mask?
     254        // use pmModelSub because modelFlux has not been generated
     255        assert (source->maskObj);
     256        psImageKeepCircle (source->maskObj, x, y, radius, "OR", markVal);
     257        pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL, maskVal);
     258        psImageKeepCircle (source->maskObj, x, y, radius, "AND", PS_NOT_IMAGE_MASK(markVal));
    259259    }
    260260
    261261    FILE *f = fopen ("shapes.dat", "w");
    262262    for (int i = 0; i < try->sources->n; i++) {
    263         psF32 inPar[10];  // must be psF32 to pmPSF_FitToModel
    264 
    265         // masked for: bad model fit, outlier in parameters
    266         if (try->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) continue;
    267 
    268         pmSource *source = try->sources->data[i];
    269         psF32 *outPar = source->modelEXT->params->data.F32;
    270 
    271         psEllipseShape shape;
    272 
    273         shape.sx  = outPar[PM_PAR_SXX] / M_SQRT2;
    274         shape.sy  = outPar[PM_PAR_SYY] / M_SQRT2;
    275         shape.sxy = outPar[PM_PAR_SXY];
    276 
    277         psEllipsePol pol = pmPSF_ModelToFit (outPar);
    278         inPar[PM_PAR_E0] = pol.e0;
    279         inPar[PM_PAR_E1] = pol.e1;
    280         inPar[PM_PAR_E2] = pol.e2;
    281         pmPSF_FitToModel (inPar, 0.1);
    282 
    283         psEllipseAxes axes1 = psEllipseShapeToAxes (shape, 20.0);
    284         psEllipseAxes axes2 = psEllipsePolToAxes(pol, 0.1);
    285 
    286         psEllipsePol pol2 = psEllipseAxesToPol (axes1);
    287 
    288         fprintf (f, "%3d  %7.2f %7.2f  %7.4f %7.4f %7.4f  --  %7.4f %7.4f %7.4f  :  %7.4f %7.4f %7.4f  --  %7.4f %7.4f %7.4f : %7.4f %7.4f %6.1f : %7.4f %7.4f %6.1f\n",
    289                 i, outPar[PM_PAR_XPOS], outPar[PM_PAR_YPOS],
    290                 outPar[PM_PAR_SXX], outPar[PM_PAR_SXY], outPar[PM_PAR_SYY],
    291                 pol.e0, pol.e1, pol.e2,
    292                 pol2.e0, pol2.e1, pol2.e2,
    293                 inPar[PM_PAR_SXX], inPar[PM_PAR_SXY], inPar[PM_PAR_SYY],
    294                 axes1.major, axes1.minor, axes1.theta*PM_DEG_RAD,
    295                 axes2.major, axes2.minor, axes2.theta*PM_DEG_RAD
    296             );
     263        psF32 inPar[10];  // must be psF32 to pmPSF_FitToModel
     264
     265        // masked for: bad model fit, outlier in parameters
     266        if (try->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) continue;
     267
     268        pmSource *source = try->sources->data[i];
     269        psF32 *outPar = source->modelEXT->params->data.F32;
     270
     271        psEllipseShape shape;
     272
     273        shape.sx  = outPar[PM_PAR_SXX] / M_SQRT2;
     274        shape.sy  = outPar[PM_PAR_SYY] / M_SQRT2;
     275        shape.sxy = outPar[PM_PAR_SXY];
     276
     277        psEllipsePol pol = pmPSF_ModelToFit (outPar);
     278        inPar[PM_PAR_E0] = pol.e0;
     279        inPar[PM_PAR_E1] = pol.e1;
     280        inPar[PM_PAR_E2] = pol.e2;
     281        pmPSF_FitToModel (inPar, 0.1);
     282
     283        psEllipseAxes axes1 = psEllipseShapeToAxes (shape, 20.0);
     284        psEllipseAxes axes2 = psEllipsePolToAxes(pol, 0.1);
     285
     286        psEllipsePol pol2 = psEllipseAxesToPol (axes1);
     287
     288        fprintf (f, "%3d  %7.2f %7.2f  %7.4f %7.4f %7.4f  --  %7.4f %7.4f %7.4f  :  %7.4f %7.4f %7.4f  --  %7.4f %7.4f %7.4f : %7.4f %7.4f %6.1f : %7.4f %7.4f %6.1f\n",
     289                i, outPar[PM_PAR_XPOS], outPar[PM_PAR_YPOS],
     290                outPar[PM_PAR_SXX], outPar[PM_PAR_SXY], outPar[PM_PAR_SYY],
     291                pol.e0, pol.e1, pol.e2,
     292                pol2.e0, pol2.e1, pol2.e2,
     293                inPar[PM_PAR_SXX], inPar[PM_PAR_SXY], inPar[PM_PAR_SYY],
     294                axes1.major, axes1.minor, axes1.theta*PM_DEG_RAD,
     295                axes2.major, axes2.minor, axes2.theta*PM_DEG_RAD
     296            );
    297297    }
    298298    fclose (f);

  • trunk/psphot/src/psphotRadialProfile.c

    r21183 r21366  
    1111    flux->data.F32[A] = flux->data.F32[B]; \
    1212    flux->data.F32[B] = tmp; \
    13     tmp = weight->data.F32[A]; \
    14     weight->data.F32[A] = weight->data.F32[B]; \
    15     weight->data.F32[B] = tmp; \
     13    tmp = variance->data.F32[A]; \
     14    variance->data.F32[A] = variance->data.F32[B]; \
     15    variance->data.F32[B] = tmp; \
    1616  } \
    1717}
    … …  
    3131    source->extpars->profile->radius = psVectorAllocEmpty (nPts, PS_TYPE_F32);
    3232    source->extpars->profile->flux   = psVectorAllocEmpty (nPts, PS_TYPE_F32);
    33     source->extpars->profile->weight = psVectorAllocEmpty (nPts, PS_TYPE_F32);
     33    source->extpars->profile->variance = psVectorAllocEmpty (nPts, PS_TYPE_F32);
    3434
    3535    psVector *radius = source->extpars->profile->radius;
    3636    psVector *flux   = source->extpars->profile->flux;
    37     psVector *weight = source->extpars->profile->weight;
     37    psVector *variance = source->extpars->profile->variance;
    3838
    3939    // XXX use the extended source model here for Xo, Yo?
    … …  
    4141
    4242    int n = 0;
    43    
     43
    4444    float Xo = 0.0;
    4545    float Yo = 0.0;
    … …  
    5757            radius->data.F32[n] = hypot (ix - Xo, iy - Yo) ;
    5858            flux->data.F32[n]   = source->pixels->data.F32[iy][ix];
    59             weight->data.F32[n] = source->weight->data.F32[iy][ix];
     59            variance->data.F32[n] = source->variance->data.F32[iy][ix];
    6060            n++;
    6161        }
    6262    }
    6363    radius->n = n;
    64     weight->n = n;
     64    variance->n = n;
    6565    flux->n = n;
    6666

  • trunk/psphot/src/psphotReadout.c

    r21250 r21366  
    4040
    4141    // Generate the mask and weight images, including the user-defined analysis region of interest
    42     psphotSetMaskAndWeight (config, readout, recipe);
     42    psphotSetMaskAndVariance (config, readout, recipe);
    4343    if (!strcasecmp (breakPt, "NOTHING")) {
    4444        return psphotReadoutCleanup(config, readout, recipe, NULL, NULL, NULL);

  • trunk/psphot/src/psphotReadoutFindPSF.c

    r21249 r21366  
    1818    PS_ASSERT_PTR_NON_NULL (readout, false);
    1919
    20     // Generate the mask and weight images, including the user-defined analysis region of interest
    21     psphotSetMaskAndWeight (config, readout, recipe);
     20    // Generate the mask and variance images, including the user-defined analysis region of interest
     21    psphotSetMaskAndVariance (config, readout, recipe);
    2222
    2323    // display the image, weight, mask (ch 1,2,3)

  • trunk/psphot/src/psphotReadoutKnownSources.c

    r21248 r21366  
    1919
    2020    // Generate the mask and weight images, including the user-defined analysis region of interest
    21     psphotSetMaskAndWeight (config, readout, recipe);
     21    psphotSetMaskAndVariance (config, readout, recipe);
    2222
    2323    // display the image, weight, mask (ch 1,2,3)
    … …  
    4141    // use the peak measured in the moments analysis:
    4242    for (int i = 0; i < sources->n; i++) {
    43         pmSource *source = sources->data[i];
    44         source->peak->flux = source->moments->Peak;
     43        pmSource *source = sources->data[i];
     44        source->peak->flux = source->moments->Peak;
    4545    }
    46    
     46
    4747    // classify sources based on moments, brightness (psf is not known)
    4848    if (!psphotRoughClass (readout, sources, recipe, false)) {
    … …  
    5353    pmPSF *psf = psphotChoosePSF (readout, sources, recipe);
    5454    if (!psf) {
    55         psLogMsg ("psphot", 3, "failure to construct a psf model");
    56         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
     55        psLogMsg ("psphot", 3, "failure to construct a psf model");
     56        return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
    5757    }
    5858    psphotVisualShowPSFModel (readout, psf);

  • trunk/psphot/src/psphotReadoutMinimal.c

    r21248 r21366  
    55// used by ppSub.
    66
    7 // NOTE: ppSub needs to perform extended source analysis for comets and trails. 
     7// NOTE: ppSub needs to perform extended source analysis for comets and trails.
    88
    99bool psphotReadoutMinimal(pmConfig *config, const pmFPAview *view) {
    … …  
    2727
    2828    // Generate the mask and weight images, including the user-defined analysis region of interest
    29     psphotSetMaskAndWeight (config, readout, recipe);
     29    psphotSetMaskAndVariance (config, readout, recipe);
    3030
    3131    // display the image, weight, mask (ch 1,2,3)
    … …  
    7878
    7979// XXX eventually, add the extended source fits here
    80 # if (0) 
     80# if (0)
    8181    // measure source size for the remaining sources
    8282    psphotSourceSize (config, readout, sources, recipe, 0);

  • trunk/psphot/src/psphotSignificanceImage.c

    r21183 r21366  
    11# include "psphotInternal.h"
    22
    3 // In this function, we smooth the image and weight, then generate the significance image :
     3// In this function, we smooth the image and variance, then generate the significance image :
    44// (S/N)^2.  If FWMH_X,Y have been recorded, use them, otherwise use PEAKS_SMOOTH_SIGMA for the
    55// smoothing kernel.
    … …  
    1010    bool guess = false;
    1111
    12     // smooth the image and weight map
     12    // smooth the image and variance map
    1313    psTimerStart ("psphot.smooth");
    1414    bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading in psImageConvolve
    … …  
    4848    psLogMsg("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark("psphot.smooth"));
    4949
    50     // Smooth the weight, applying the mask as we go.  The variance is smoothed by the PSF^2,
     50    // Smooth the variance, applying the mask as we go.  The variance is smoothed by the PSF^2,
    5151    // renomalized to maintain the input level of the variance.  We achieve this by smoothing
    5252    // with a Gaussian with sigma = SIGMA_SMTH/sqrt(2) with unity normalization.  Note that
    … …  
    5555    // measurements based on apertures comparable to or larger than the smoothing kernel, the
    5656    // effective per-pixel variance is maintained.
    57     psImage *smooth_wt = psImageCopy(NULL, readout->weight, PS_TYPE_F32);
     57    psImage *smooth_wt = psImageCopy(NULL, readout->variance, PS_TYPE_F32);
    5858    psImageSmoothMask_Threaded(smooth_wt, smooth_wt, readout->mask, maskVal, SIGMA_SMTH * M_SQRT1_2,
    5959                      NSIGMA_SMTH, minGauss);
    60     psLogMsg("psphot", PS_LOG_MINUTIA, "smooth weight: %f sec\n", psTimerMark("psphot.smooth"));
     60    psLogMsg("psphot", PS_LOG_MINUTIA, "smooth variance: %f sec\n", psTimerMark("psphot.smooth"));
    6161
    6262    psImage *mask = readout->mask;
    … …  
    100100
    101101    float factor = guess ? 4.0 * M_PI * PS_SQR(SIGMA_SMTH) : 4.0 * M_PI * PS_SQR(SIGMA_SMTH);
     102
     103    // Correct the correction factor for the covariance produced by the (potentially multiple) smoothing
     104    psKernel *kernel = psImageSmoothKernel(SIGMA_SMTH, NSIGMA_SMTH); // Kernel used for smoothing
     105    psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix
     106    psFree(kernel);
     107    factor /= psImageCovarianceFactor(covar);
     108    psFree(covar);
     109
    102110    // record the effective area and significance scaling factor
    103111    float effArea = 8.0 * M_PI * PS_SQR(SIGMA_SMTH);

  • trunk/psphot/src/psphotSourceSize.c

    r21183 r21366  
    22# include <gsl/gsl_sf_gamma.h>
    33
    4 static float psphotModelContour(const psImage *image, const psImage *weight, const psImage *mask,
     4static float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask,
    55                                psImageMaskType maskVal, const pmModel *model, float Ro);
    66
    … …  
    6262
    6363        psF32 **resid  = source->pixels->data.F32;
    64         psF32 **weight = source->weight->data.F32;
     64        psF32 **variance = source->variance->data.F32;
    6565        psImageMaskType **mask    = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
    6666
    6767        // check for extendedness: measure the delta flux significance at the 1 sigma contour
    68         source->extNsigma = psphotModelContour(source->pixels, source->weight, source->maskObj, maskVal,
     68        source->extNsigma = psphotModelContour(source->pixels, source->variance, source->maskObj, maskVal,
    6969                                               source->modelPSF, 1.0);
    7070
    … …  
    103103        // Compare the central pixel with those on either side, for the four possible lines through it.
    104104
    105         // Soften weights (add systematic error)
    106         float softening = soft * PS_SQR(source->peak->flux); // Softening for weights
     105        // Soften variances (add systematic error)
     106        float softening = soft * PS_SQR(source->peak->flux); // Softening for variances
    107107
    108108        // Across the middle: y = 0
    109109        float cX = 2*resid[yPeak][xPeak]   - resid[yPeak+0][xPeak-1]  - resid[yPeak+0][xPeak+1];
    110         float dcX = 4*weight[yPeak][xPeak] + weight[yPeak+0][xPeak-1] + weight[yPeak+0][xPeak+1];
     110        float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1];
    111111        float nX = cX / sqrtf(dcX + softening);
    112112
    113113        // Up the centre: x = 0
    114114        float cY = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak+0]  - resid[yPeak+1][xPeak+0];
    115         float dcY = 4*weight[yPeak][xPeak] + weight[yPeak-1][xPeak+0] + weight[yPeak+1][xPeak+0];
     115        float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0];
    116116        float nY = cY / sqrtf(dcY + softening);
    117117
    118118        // Diagonal: x = y
    119119        float cL = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak-1]  - resid[yPeak+1][xPeak+1];
    120         float dcL = 4*weight[yPeak][xPeak] + weight[yPeak-1][xPeak-1] + weight[yPeak+1][xPeak+1];
     120        float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1];
    121121        float nL = cL / sqrtf(dcL + softening);
    122122
    123123        // Diagonal: x = - y
    124124        float cR = 2*resid[yPeak][xPeak]   - resid[yPeak+1][xPeak-1]  - resid[yPeak-1][xPeak+1];
    125         float dcR = 4*weight[yPeak][xPeak] + weight[yPeak+1][xPeak-1] + weight[yPeak-1][xPeak+1];
     125        float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1];
    126126        float nR = cR / sqrtf(dcR + softening);
    127127
    … …  
    160160        // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
    161161        if (source->crNsigma > CR_NSIGMA_LIMIT) {
    162             // XXX still testing... : psphotMaskCosmicRay_New (readout->mask, source, maskVal, crMask);
    163             psphotMaskCosmicRay_Old (source, maskVal, crMask);
     162            // XXX still testing... : psphotMaskCosmicRay_New (readout->mask, source, maskVal, crMask);
     163            psphotMaskCosmicRay_Old (source, maskVal, crMask);
    164164        }
    165165    }
    … …  
    190190// deviation in sigmas.  This is measured on the residual image - should we ignore negative
    191191// deviations?
    192 static float psphotModelContour(const psImage *image, const psImage *weight, const psImage *mask,
     192static float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask,
    193193                                psImageMaskType maskVal, const pmModel *model, float Ro)
    194194{
    … …  
    240240        if (yPixM >= 0 && yPixM < image->numRows &&
    241241            !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) {
    242             float dSigma = image->data.F32[yPixM][xPix] / sqrtf(weight->data.F32[yPixM][xPix]);
     242            float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]);
    243243            nSigma += dSigma;
    244244            nPts++;
    … …  
    251251        if (yPixP >= 0 && yPixP < image->numRows &&
    252252            !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) {
    253             float dSigma = image->data.F32[yPixP][xPix] / sqrtf(weight->data.F32[yPixP][xPix]);
     253            float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]);
    254254            nSigma += dSigma;
    255255            nPts++;
    … …  
    274274    pmFootprint *footprint = peak->footprint;
    275275    if (!footprint) {
    276         // if we have not footprint, use the old code to mask by isophot
    277         psphotMaskCosmicRay_Old (source, maskVal, crMask);
    278         return true;
     276        // if we have not footprint, use the old code to mask by isophot
     277        psphotMaskCosmicRay_Old (source, maskVal, crMask);
     278        return true;
    279279    }
    280280
    281281    if (!footprint->spans) {
    282         // if we have not footprint, use the old code to mask by isophot
    283         psphotMaskCosmicRay_Old (source, maskVal, crMask);
    284         return true;
     282        // if we have not footprint, use the old code to mask by isophot
     283        psphotMaskCosmicRay_Old (source, maskVal, crMask);
     284        return true;
    285285    }
    286286
    287287    // mask all of the pixels covered by the spans of the footprint
    288288    for (int j = 1; j < footprint->spans->n; j++) {
    289         pmSpan *span1 = footprint->spans->data[j];
    290 
    291         int iy = span1->y;
    292         int xs = span1->x0;
    293         int xe = span1->x1;
    294        
    295         for (int ix = xs; ix < xe; ix++) {
    296             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
    297         }
     289        pmSpan *span1 = footprint->spans->data[j];
     290
     291        int iy = span1->y;
     292        int xs = span1->x0;
     293        int xe = span1->x1;
     294
     295        for (int ix = xs; ix < xe; ix++) {
     296            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
     297        }
    298298    }
    299299    return true;
    … …  
    308308    psImage *mask   = source->maskView;
    309309    psImage *pixels = source->pixels;
    310     psImage *weight = source->weight;
     310    psImage *variance = source->variance;
    311311
    312312    // XXX This should be a recipe variable
    … …  
    318318    // mark the pixels in this row to the left, then the right
    319319    for (int ix = xo; ix >= 0; ix--) {
    320         float SN = pixels->data.F32[yo][ix] / sqrt(weight->data.F32[yo][ix]);
    321         if (SN > SN_LIMIT) {
    322             mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
    323         }
     320        float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
     321        if (SN > SN_LIMIT) {
     322            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
     323        }
    324324    }
    325325    for (int ix = xo + 1; ix < pixels->numCols; ix++) {
    326         float SN = pixels->data.F32[yo][ix] / sqrt(weight->data.F32[yo][ix]);
    327         if (SN > SN_LIMIT) {
    328             mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
    329         }
     326        float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
     327        if (SN > SN_LIMIT) {
     328            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
     329        }
    330330    }
    331331
    … …  
    333333    // first go up:
    334334    for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) {
    335         // mark the pixels in this row to the left, then the right
    336         for (int ix = 0; ix < pixels->numCols; ix++) {
    337             float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]);
    338             if (SN < SN_LIMIT) continue;
    339 
    340             bool valid = false;
    341             valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask);
    342             valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0;
    343             valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0;
    344 
    345             if (!valid) continue;
    346             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
    347         }
     335        // mark the pixels in this row to the left, then the right
     336        for (int ix = 0; ix < pixels->numCols; ix++) {
     337            float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
     338            if (SN < SN_LIMIT) continue;
     339
     340            bool valid = false;
     341            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask);
     342            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0;
     343            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0;
     344
     345            if (!valid) continue;
     346            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
     347        }
    348348    }
    349349    // next go down:
    350350    for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) {
    351         // mark the pixels in this row to the left, then the right
    352         for (int ix = 0; ix < pixels->numCols; ix++) {
    353             float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]);
    354             if (SN < SN_LIMIT) continue;
    355 
    356             bool valid = false;
    357             valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask);
    358             valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0;
    359             valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0;
    360 
    361             if (!valid) continue;
    362             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
    363         }
     351        // mark the pixels in this row to the left, then the right
     352        for (int ix = 0; ix < pixels->numCols; ix++) {
     353            float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
     354            if (SN < SN_LIMIT) continue;
     355
     356            bool valid = false;
     357            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask);
     358            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0;
     359            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0;
     360
     361            if (!valid) continue;
     362            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
     363        }
    364364    }
    365365    return true;

  • trunk/psphot/src/psphotVisual.c

    r21183 r21366  
    22
    33// this function displays representative images as the psphot analysis progresses:
    4 // 0 : image, 1 : weight
    5 // 0 : backsub, 1 : weight, 2 : backgnd
    6 // 0 : backsub, 1 : weight, 2 : signif
     4// 0 : image, 1 : variance
     5// 0 : backsub, 1 : variance, 2 : backgnd
     6// 0 : backsub, 1 : variance, 2 : signif
    77// (overlay peaks on images)
    88// (overlay footprints on images)
    99// (overlay moments on images)
    10 // (overlay rough class on images) 
     10// (overlay rough class on images)
    1111// 0 : backsub, 1 : psfpos, 2: psfsub
    1212// 0 : backsub, 1 : lin_resid, 2: psfsub
    … …  
    4040    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
    4141    if (!psImageBackground(stats, NULL, inImage, NULL, 0, rng)) {
    42         fprintf (stderr, "failed to get background values\n");
    43         return false;
     42        fprintf (stderr, "failed to get background values\n");
     43        return false;
    4444    }
    4545
    … …  
    4949    ALLOCATE (image.data2d, float *, image.Ny);
    5050    for (int iy = 0; iy < image.Ny; iy++) {
    51         ALLOCATE (image.data2d[iy], float, image.Nx);
    52         for (int ix = 0; ix < image.Nx; ix++) {
    53             image.data2d[iy][ix] = inImage->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix];
    54         }
     51        ALLOCATE (image.data2d[iy], float, image.Nx);
     52        for (int ix = 0; ix < image.Nx; ix++) {
     53            image.data2d[iy][ix] = inImage->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix];
     54        }
    5555    }
    5656
    … …  
    6060    data.range = 32;
    6161    data.logflux = 0;
    62  
     62
    6363    KiiSetChannel (kapaFD, channel);
    6464    KiiNewPicture2D (kapaFD, &image, &data, &coords);
    6565
    6666    for (int iy = 0; iy < image.Ny; iy++) {
    67         free (image.data2d[iy]);
     67        free (image.data2d[iy]);
    6868    }
    6969    free (image.data2d);
    … …  
    8686    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
    8787    if (!psImageBackground(stats, NULL, inImage, NULL, 0, rng)) {
    88         fprintf (stderr, "failed to get background values\n");
    89         return false;
     88        fprintf (stderr, "failed to get background values\n");
     89        return false;
    9090    }
    9191
    … …  
    9999    data.range = 5*stats->robustStdev;
    100100    data.logflux = 0;
    101  
     101
    102102    KiiSetChannel (kapaFD, channel);
    103103    KiiNewPicture2D (kapaFD, &image, &data, &coords);
    … …  
    126126    data.range = max - min;
    127127    data.logflux = 0;
    128  
     128
    129129    KiiSetChannel (kapaFD, channel);
    130130    KiiNewPicture2D (kapaFD, &image, &data, &coords);
    … …  
    139139    if (kapa == -1) {
    140140        kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
    141         if (kapa == -1) {
    142             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    143             isVisual = false;
    144             return false;
    145         }
    146     } 
     141        if (kapa == -1) {
     142            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     143            isVisual = false;
     144            return false;
     145        }
     146    }
    147147
    148148    // psphotVisualShowMask (kapa, readout->mask, "mask", 2);
    149     psphotVisualScaleImage (kapa, readout->weight, "weight", 1);
     149    psphotVisualScaleImage (kapa, readout->variance, "variance", 1);
    150150    psphotVisualScaleImage (kapa, readout->image, "image", 0);
    151151
    … …  
    155155    fprintf (stdout, "[c]ontinue? ");
    156156    if (!fgets(key, 8, stdin)) {
    157         psWarning("Unable to read option");
     157        psWarning("Unable to read option");
    158158    }
    159159    return true;
    … …  
    168168    if (kapa == -1) {
    169169        kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
    170         if (kapa == -1) {
    171             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    172             isVisual = false;
    173             return false;
    174         }
    175     } 
     170        if (kapa == -1) {
     171            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     172            isVisual = false;
     173            return false;
     174        }
     175    }
    176176
    177177    bool status = false;
    … …  
    179179
    180180    if (file->mode == PM_FPA_MODE_INTERNAL) {
    181         backgnd = file->readout;
     181        backgnd = file->readout;
    182182    } else {
    183         backgnd = pmFPAviewThisReadout (view, file->fpa);
     183        backgnd = pmFPAviewThisReadout (view, file->fpa);
    184184    }
    185185
    … …  
    192192    fprintf (stdout, "[c]ontinue? ");
    193193    if (!fgets(key, 8, stdin)) {
    194         psWarning("Unable to read option");
     194        psWarning("Unable to read option");
    195195    }
    196196    return true;
    … …  
    203203    if (kapa == -1) {
    204204        kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
    205         if (kapa == -1) {
    206             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    207             isVisual = false;
    208             return false;
    209         }
    210     } 
     205        if (kapa == -1) {
     206            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     207            isVisual = false;
     208            return false;
     209        }
     210    }
    211211
    212212    // XXX test: image->data.F32[10][10] = 10000;
    … …  
    218218    fprintf (stdout, "[c]ontinue? ");
    219219    if (!fgets(key, 8, stdin)) {
    220         psWarning("Unable to read option");
     220        psWarning("Unable to read option");
    221221    }
    222222    return true;
    … …  
    227227    int Noverlay;
    228228    KiiOverlay *overlay;
    229  
     229
    230230    if (!isVisual) return true;
    231231
    232232    if (kapa == -1) {
    233         fprintf (stderr, "kapa not opened, skipping\n");
    234         return false;
     233        fprintf (stderr, "kapa not opened, skipping\n");
     234        return false;
    235235    }
    236236
    … …  
    244244    for (int i = 0; i < peaks->n; i++) {
    245245
    246         pmPeak *peak = peaks->data[i];
    247         if (peak == NULL) continue;
    248 
    249         overlay[Noverlay].type = KII_OVERLAY_BOX;
    250         overlay[Noverlay].x = peak->xf;
    251         overlay[Noverlay].y = peak->yf;
    252         overlay[Noverlay].dx = 2.0;
    253         overlay[Noverlay].dy = 2.0;
    254         overlay[Noverlay].angle = 0.0;
    255         overlay[Noverlay].text = NULL;
    256         Noverlay ++;
     246        pmPeak *peak = peaks->data[i];
     247        if (peak == NULL) continue;
     248
     249        overlay[Noverlay].type = KII_OVERLAY_BOX;
     250        overlay[Noverlay].x = peak->xf;
     251        overlay[Noverlay].y = peak->yf;
     252        overlay[Noverlay].dx = 2.0;
     253        overlay[Noverlay].dy = 2.0;
     254        overlay[Noverlay].angle = 0.0;
     255        overlay[Noverlay].text = NULL;
     256        Noverlay ++;
    257257
    258258# if (0)
    259         overlay[Noverlay].type = KII_OVERLAY_BOX;
    260         overlay[Noverlay].x = peak->x;
    261         overlay[Noverlay].y = peak->y;
    262         overlay[Noverlay].dx = 1.0;
    263         overlay[Noverlay].dy = 1.0;
    264         overlay[Noverlay].angle = 0.0;
    265         overlay[Noverlay].text = NULL;
    266         Noverlay ++;
    267 
    268         overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
    269         overlay[Noverlay].x = peak->xf;
    270         overlay[Noverlay].y = peak->yf;
    271         overlay[Noverlay].dx = 2.0;
    272         overlay[Noverlay].dy = 2.0;
    273         overlay[Noverlay].angle = 0.0;
    274         overlay[Noverlay].text = NULL;
    275         Noverlay ++;
     259        overlay[Noverlay].type = KII_OVERLAY_BOX;
     260        overlay[Noverlay].x = peak->x;
     261        overlay[Noverlay].y = peak->y;
     262        overlay[Noverlay].dx = 1.0;
     263        overlay[Noverlay].dy = 1.0;
     264        overlay[Noverlay].angle = 0.0;
     265        overlay[Noverlay].text = NULL;
     266        Noverlay ++;
     267
     268        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     269        overlay[Noverlay].x = peak->xf;
     270        overlay[Noverlay].y = peak->yf;
     271        overlay[Noverlay].dx = 2.0;
     272        overlay[Noverlay].dy = 2.0;
     273        overlay[Noverlay].angle = 0.0;
     274        overlay[Noverlay].text = NULL;
     275        Noverlay ++;
    276276# endif
    277277    }
    … …  
    296296    fprintf (stdout, "[c]ontinue? ");
    297297    if (!fgets(key, 8, stdin)) {
    298         psWarning("Unable to read option");
     298        psWarning("Unable to read option");
    299299    }
    300300    return true;
    … …  
    305305    int Noverlay;
    306306    KiiOverlay *overlay;
    307  
     307
    308308    if (!isVisual) return true;
    309309
    310310    if (kapa == -1) {
    311         fprintf (stderr, "kapa not opened, skipping\n");
    312         return false;
     311        fprintf (stderr, "kapa not opened, skipping\n");
     312        return false;
    313313    }
    314314
    … …  
    323323    for (int i = 0; i < footprints->n; i++) {
    324324
    325         pmSpan *span = NULL;
    326 
    327         pmFootprint *footprint = footprints->data[i];
    328         if (footprint == NULL) continue;
    329         if (footprint->spans == NULL) continue;
    330         if (footprint->spans->n < 1) continue;
    331 
    332         // draw the top
    333         span = footprint->spans->data[0];
    334         overlay[Noverlay].type = KII_OVERLAY_LINE;
    335         overlay[Noverlay].x = span->x0;
    336         overlay[Noverlay].y = span->y;
    337         overlay[Noverlay].dx = span->x1 - span->x0;
    338         overlay[Noverlay].dy = 0;
    339         overlay[Noverlay].angle = 0.0;
    340         overlay[Noverlay].text = NULL;
    341         Noverlay ++;
    342         CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    343 
    344         int ys = span->y;
    345         int x0s = span->x0;
    346         int x1s = span->x1;
    347 
    348         // draw the outer span edges
    349         for (int j = 1; j < footprint->spans->n; j++) {
    350             pmSpan *span1 = footprint->spans->data[j];
    351 
    352             int ye = span1->y;
    353             int x0e = span1->x0;
    354             int x1e = span1->x1;
    355 
    356             // we cannot have two discontinuous spans on the top or bottom, right? (no, probably not right)
    357             // find all of the spans in this row and generate x0e, x01:
    358             for (int k = j + 1; k < footprint->spans->n; k++) {
    359                 pmSpan *span2 = footprint->spans->data[k];
    360                 if (span2->y > span1->y) break;
    361                 x0e = PS_MIN (x0e, span2->x0);
    362                 x1e = PS_MAX (x1e, span2->x1);
    363                 j++;
    364             }
    365 
    366             overlay[Noverlay].type = KII_OVERLAY_LINE;
    367             overlay[Noverlay].x = x0s;
    368             overlay[Noverlay].y = ys;
    369             overlay[Noverlay].dx = x0e - x0s;
    370             overlay[Noverlay].dy = ye - ys;
    371             overlay[Noverlay].angle = 0.0;
    372             overlay[Noverlay].text = NULL;
    373             Noverlay ++;
    374             CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    375 
    376             overlay[Noverlay].type = KII_OVERLAY_LINE;
    377             overlay[Noverlay].x = x1s;
    378             overlay[Noverlay].y = ys;
    379             overlay[Noverlay].dx = x1e - x1s;
    380             overlay[Noverlay].dy = ye - ys;
    381             overlay[Noverlay].angle = 0.0;
    382             overlay[Noverlay].text = NULL;
    383             Noverlay ++;
    384             CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    385 
    386             ys = ye;
    387             x0s = x0e;
    388             x1s = x1e;
    389         }
    390 
    391         // draw the bottom
    392         span = footprint->spans->data[footprint->spans->n - 1];
    393         overlay[Noverlay].type = KII_OVERLAY_LINE;
    394         overlay[Noverlay].x = span->x0;
    395         overlay[Noverlay].y = span->y;
    396         overlay[Noverlay].dx = span->x1 - span->x0;
    397         overlay[Noverlay].dy = 0;
    398         overlay[Noverlay].angle = 0.0;
    399         overlay[Noverlay].text = NULL;
    400         Noverlay ++;
    401         CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     325        pmSpan *span = NULL;
     326
     327        pmFootprint *footprint = footprints->data[i];
     328        if (footprint == NULL) continue;
     329        if (footprint->spans == NULL) continue;
     330        if (footprint->spans->n < 1) continue;
     331
     332        // draw the top
     333        span = footprint->spans->data[0];
     334        overlay[Noverlay].type = KII_OVERLAY_LINE;
     335        overlay[Noverlay].x = span->x0;
     336        overlay[Noverlay].y = span->y;
     337        overlay[Noverlay].dx = span->x1 - span->x0;
     338        overlay[Noverlay].dy = 0;
     339        overlay[Noverlay].angle = 0.0;
     340        overlay[Noverlay].text = NULL;
     341        Noverlay ++;
     342        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     343
     344        int ys = span->y;
     345        int x0s = span->x0;
     346        int x1s = span->x1;
     347
     348        // draw the outer span edges
     349        for (int j = 1; j < footprint->spans->n; j++) {
     350            pmSpan *span1 = footprint->spans->data[j];
     351
     352            int ye = span1->y;
     353            int x0e = span1->x0;
     354            int x1e = span1->x1;
     355
     356            // we cannot have two discontinuous spans on the top or bottom, right? (no, probably not right)
     357            // find all of the spans in this row and generate x0e, x01:
     358            for (int k = j + 1; k < footprint->spans->n; k++) {
     359                pmSpan *span2 = footprint->spans->data[k];
     360                if (span2->y > span1->y) break;
     361                x0e = PS_MIN (x0e, span2->x0);
     362                x1e = PS_MAX (x1e, span2->x1);
     363                j++;
     364            }
     365
     366            overlay[Noverlay].type = KII_OVERLAY_LINE;
     367            overlay[Noverlay].x = x0s;
     368            overlay[Noverlay].y = ys;
     369            overlay[Noverlay].dx = x0e - x0s;
     370            overlay[Noverlay].dy = ye - ys;
     371            overlay[Noverlay].angle = 0.0;
     372            overlay[Noverlay].text = NULL;
     373            Noverlay ++;
     374            CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     375
     376            overlay[Noverlay].type = KII_OVERLAY_LINE;
     377            overlay[Noverlay].x = x1s;
     378            overlay[Noverlay].y = ys;
     379            overlay[Noverlay].dx = x1e - x1s;
     380            overlay[Noverlay].dy = ye - ys;
     381            overlay[Noverlay].angle = 0.0;
     382            overlay[Noverlay].text = NULL;
     383            Noverlay ++;
     384            CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     385
     386            ys = ye;
     387            x0s = x0e;
     388            x1s = x1e;
     389        }
     390
     391        // draw the bottom
     392        span = footprint->spans->data[footprint->spans->n - 1];
     393        overlay[Noverlay].type = KII_OVERLAY_LINE;
     394        overlay[Noverlay].x = span->x0;
     395        overlay[Noverlay].y = span->y;
     396        overlay[Noverlay].dx = span->x1 - span->x0;
     397        overlay[Noverlay].dy = 0;
     398        overlay[Noverlay].angle = 0.0;
     399        overlay[Noverlay].text = NULL;
     400        Noverlay ++;
     401        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    402402    }
    403403
    … …  
    410410    fprintf (stdout, "[c]ontinue? ");
    411411    if (!fgets(key, 8, stdin)) {
    412         psWarning("Unable to read option");
     412        psWarning("Unable to read option");
    413413    }
    414414    return true;
    … …  
    419419    int Noverlay;
    420420    KiiOverlay *overlay;
    421  
     421
    422422    psEllipseMoments emoments;
    423423    psEllipseAxes axes;
    … …  
    426426
    427427    if (kapa == -1) {
    428         fprintf (stderr, "kapa not opened, skipping\n");
    429         return false;
     428        fprintf (stderr, "kapa not opened, skipping\n");
     429        return false;
    430430    }
    431431
    … …  
    436436    for (int i = 0; i < sources->n; i++) {
    437437
    438         pmSource *source = sources->data[i];
    439         if (source == NULL) continue;
    440 
    441         pmMoments *moments = source->moments;
    442         if (moments == NULL) continue;
    443 
    444         overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
    445         overlay[Noverlay].x = moments->Mx;
    446         overlay[Noverlay].y = moments->My;
    447 
    448         emoments.x2 = moments->Mxx;
    449         emoments.xy = moments->Mxy;
    450         emoments.y2 = moments->Myy;
    451 
    452         axes = psEllipseMomentsToAxes (emoments, 20.0);
    453 
    454         overlay[Noverlay].dx = 2.0*axes.major;
    455         overlay[Noverlay].dy = 2.0*axes.minor;
    456         overlay[Noverlay].angle = -axes.theta * PS_DEG_RAD;  // XXXXXXXX the axes angle is negative to display of object on kapa
    457         overlay[Noverlay].text = NULL;
    458         Noverlay ++;
     438        pmSource *source = sources->data[i];
     439        if (source == NULL) continue;
     440
     441        pmMoments *moments = source->moments;
     442        if (moments == NULL) continue;
     443
     444        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     445        overlay[Noverlay].x = moments->Mx;
     446        overlay[Noverlay].y = moments->My;
     447
     448        emoments.x2 = moments->Mxx;
     449        emoments.xy = moments->Mxy;
     450        emoments.y2 = moments->Myy;
     451
     452        axes = psEllipseMomentsToAxes (emoments, 20.0);
     453
     454        overlay[Noverlay].dx = 2.0*axes.major;
     455        overlay[Noverlay].dy = 2.0*axes.minor;
     456        overlay[Noverlay].angle = -axes.theta * PS_DEG_RAD;  // XXXXXXXX the axes angle is negative to display of object on kapa
     457        overlay[Noverlay].text = NULL;
     458        Noverlay ++;
    459459    }
    460460
    … …  
    467467    fprintf (stdout, "[c]ontinue? ");
    468468    if (!fgets(key, 8, stdin)) {
    469         psWarning("Unable to read option");
     469        psWarning("Unable to read option");
    470470    }
    471471
    … …  
    482482    if (kapa3 == -1) {
    483483        kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots");
    484         if (kapa3 == -1) {
    485             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    486             isVisual = false;
    487             return false;
    488         }
    489     } 
     484        if (kapa3 == -1) {
     485            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     486            isVisual = false;
     487            return false;
     488        }
     489    }
    490490
    491491    KapaClearPlots (kapa3);
    … …  
    494494    // there are N regions: use the first (guaranteed to exist) to get the overall limits
    495495    psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, "PSF.CLUMP.REGION.000");
    496    
     496
    497497    float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
    498498    float psfY  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
    … …  
    558558
    559559    // XXX draw N circles to outline the clumps
    560     { 
    561         // draw a circle centered on psfX,Y with size of the psf limit
    562         psVector *xLimit  = psVectorAlloc (120, PS_TYPE_F32);
    563         psVector *yLimit  = psVectorAlloc (120, PS_TYPE_F32);
    564 
    565         int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
    566         float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
    567 
    568         graphdata.color = KapaColorByName ("blue");
    569         graphdata.style = 0;
    570 
    571         for (int n = 0; n < nRegions; n++) {
    572 
    573             char regionName[64];
    574             snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
    575             psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
    576    
    577             float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
    578             float psfY  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
    579             float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
    580             float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
    581             float Rx = psfdX * PSF_CLUMP_NSIGMA;
    582             float Ry = psfdY * PSF_CLUMP_NSIGMA;
    583 
    584             for (int i = 0; i < xLimit->n; i++) {
    585                 xLimit->data.F32[i] = Rx*cos(i*2.0*M_PI/120.0) + psfX;
    586                 yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY;
    587             }
    588             KapaPrepPlot (kapa3, xLimit->n, &graphdata);
    589             KapaPlotVector (kapa3, xLimit->n, xLimit->data.F32, "x");
    590             KapaPlotVector (kapa3, yLimit->n, yLimit->data.F32, "y");
    591         }
    592         psFree (xLimit);
    593         psFree (yLimit);
     560    {
     561        // draw a circle centered on psfX,Y with size of the psf limit
     562        psVector *xLimit  = psVectorAlloc (120, PS_TYPE_F32);
     563        psVector *yLimit  = psVectorAlloc (120, PS_TYPE_F32);
     564
     565        int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
     566        float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
     567
     568        graphdata.color = KapaColorByName ("blue");
     569        graphdata.style = 0;
     570
     571        for (int n = 0; n < nRegions; n++) {
     572
     573            char regionName[64];
     574            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
     575            psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
     576
     577            float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
     578            float psfY  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
     579            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
     580            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
     581            float Rx = psfdX * PSF_CLUMP_NSIGMA;
     582            float Ry = psfdY * PSF_CLUMP_NSIGMA;
     583
     584            for (int i = 0; i < xLimit->n; i++) {
     585                xLimit->data.F32[i] = Rx*cos(i*2.0*M_PI/120.0) + psfX;
     586                yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY;
     587            }
     588            KapaPrepPlot (kapa3, xLimit->n, &graphdata);
     589            KapaPlotVector (kapa3, xLimit->n, xLimit->data.F32, "x");
     590            KapaPlotVector (kapa3, yLimit->n, yLimit->data.F32, "y");
     591        }
     592        psFree (xLimit);
     593        psFree (yLimit);
    594594    }
    595595
    596596# if (0)
    597     // *** make a histogram of the source counts in the x and y directions 
     597    // *** make a histogram of the source counts in the x and y directions
    598598    psHistogram *nX = psHistogramAlloc (graphdata.xmin, graphdata.xmax, 50.0);
    599599    psHistogram *nY = psHistogramAlloc (graphdata.ymin, graphdata.ymax, 50.0);
    … …  
    605605    psVector *vY = psVectorAlloc (nY->nums->n, PS_TYPE_F32);
    606606    for (int i = 0; i < nX->nums->n; i++) {
    607         dX->data.F32[i] = nX->nums->data.S32[i];
    608         vX->data.F32[i] = 0.5*(nX->bounds->data.F32[i] + nX->bounds->data.F32[i+1]);
     607        dX->data.F32[i] = nX->nums->data.S32[i];
     608        vX->data.F32[i] = 0.5*(nX->bounds->data.F32[i] + nX->bounds->data.F32[i+1]);
    609609    }
    610610    for (int i = 0; i < nY->nums->n; i++) {
    611         dY->data.F32[i] = nY->nums->data.S32[i];
    612         vY->data.F32[i] = 0.5*(nY->bounds->data.F32[i] + nY->bounds->data.F32[i+1]);
     611        dY->data.F32[i] = nY->nums->data.S32[i];
     612        vY->data.F32[i] = 0.5*(nY->bounds->data.F32[i] + nY->bounds->data.F32[i+1]);
    613613    }
    614614
    … …  
    640640    fprintf (stdout, "[c]ontinue? ");
    641641    if (!fgets(key, 8, stdin)) {
    642         psWarning("Unable to read option");
     642        psWarning("Unable to read option");
    643643    }
    644644    return true;
    … …  
    650650    int Noverlay;
    651651    KiiOverlay *overlay;
    652  
     652
    653653    psEllipseMoments emoments;
    654654    psEllipseAxes axes;
    … …  
    660660    for (int i = 0; i < sources->n; i++) {
    661661
    662         pmSource *source = sources->data[i];
    663         if (source == NULL) continue;
    664 
    665         if (source->type != type) continue;
    666         if (mode && !(source->mode & mode)) continue;
    667 
    668         pmMoments *moments = source->moments;
    669         if (moments == NULL) continue;
    670 
    671         overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
    672         overlay[Noverlay].x = moments->Mx;
    673         overlay[Noverlay].y = moments->My;
    674 
    675         emoments.x2 = moments->Mxx;
    676         emoments.y2 = moments->Myy;
    677         emoments.xy = moments->Mxy;
    678 
    679         axes = psEllipseMomentsToAxes (emoments, 20.0);
    680 
    681         overlay[Noverlay].dx = 2.0*axes.major;
    682         overlay[Noverlay].dy = 2.0*axes.minor;
    683         overlay[Noverlay].angle = -axes.theta * PS_DEG_RAD;
    684         overlay[Noverlay].text = NULL;
    685         Noverlay ++;
     662        pmSource *source = sources->data[i];
     663        if (source == NULL) continue;
     664
     665        if (source->type != type) continue;
     666        if (mode && !(source->mode & mode)) continue;
     667
     668        pmMoments *moments = source->moments;
     669        if (moments == NULL) continue;
     670
     671        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     672        overlay[Noverlay].x = moments->Mx;
     673        overlay[Noverlay].y = moments->My;
     674
     675        emoments.x2 = moments->Mxx;
     676        emoments.y2 = moments->Myy;
     677        emoments.xy = moments->Mxy;
     678
     679        axes = psEllipseMomentsToAxes (emoments, 20.0);
     680
     681        overlay[Noverlay].dx = 2.0*axes.major;
     682        overlay[Noverlay].dy = 2.0*axes.minor;
     683        overlay[Noverlay].angle = -axes.theta * PS_DEG_RAD;
     684        overlay[Noverlay].text = NULL;
     685        Noverlay ++;
    686686    }
    687687
    … …  
    697697
    698698    if (kapa == -1) {
    699         fprintf (stderr, "kapa not opened, skipping\n");
    700         return false;
     699        fprintf (stderr, "kapa not opened, skipping\n");
     700        return false;
    701701    }
    702702
    … …  
    716716    fprintf (stdout, "[c]ontinue? ");
    717717    if (!fgets(key, 8, stdin)) {
    718         psWarning("Unable to read option");
     718        psWarning("Unable to read option");
    719719    }
    720720    return true;
    … …  
    727727    if (kapa2 == -1) {
    728728        kapa2 = KapaOpenNamedSocket ("kapa", "psphot:psfstars");
    729         if (kapa2 == -1) {
    730             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    731             isVisual = false;
    732             return false;
    733         }
    734     } 
     729        if (kapa2 == -1) {
     730            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     731            isVisual = false;
     732            return false;
     733        }
     734    }
    735735
    736736    int DX = 64;
    … …  
    739739    psImage *psfMosaic = psImageAlloc (5*DX, 5*DY, PS_TYPE_F32);
    740740    psImageInit (psfMosaic, 0.0);
    741    
     741
    742742    psImage *funMosaic = psImageAlloc (5*DX, 5*DY, PS_TYPE_F32);
    743743    psImageInit (funMosaic, 0.0);
    … …  
    750750    // generate a fake model at each of the 3x3 image grid positions
    751751    for (int x = -2; x <= +2; x ++) {
    752         for (int y = -2; y <= +2; y ++) {
    753             // use the center of the center pixel of the image
    754             float xc = (int)((0.5 + 0.225*x)*readout->image->numCols) + readout->image->col0 + 0.5;
    755             float yc = (int)((0.5 + 0.225*y)*readout->image->numRows) + readout->image->row0 + 0.5;
    756 
    757             // assign the x and y coords to the image center
    758             // create an object with center intensity of 1000
    759             modelRef->params->data.F32[PM_PAR_SKY] = 0;
    760             modelRef->params->data.F32[PM_PAR_I0] = 1000;
    761             modelRef->params->data.F32[PM_PAR_XPOS] = xc;
    762             modelRef->params->data.F32[PM_PAR_YPOS] = yc;
    763    
    764             // create modelPSF from this model
    765             pmModel *model = pmModelFromPSF (modelRef, psf);
    766             if (!model) continue;
    767 
    768             // place the reference object in the image center
    769             // no need to mask the source here
    770             // XXX should we measure this for the analytical model only or the full model?
    771             pmModelAddWithOffset (psfMosaic, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_CENTER, 0, x*DX, y*DY);
    772             pmModelAddWithOffset (funMosaic, NULL, model, PM_MODEL_OP_FUNC | PM_MODEL_OP_CENTER, 0, x*DX, y*DY);
    773             pmModelAddWithOffset (resMosaic, NULL, model, PM_MODEL_OP_RES0 | PM_MODEL_OP_RES1 | PM_MODEL_OP_CENTER, 0, x*DX, y*DY);
    774             psFree (model);
    775         }
     752        for (int y = -2; y <= +2; y ++) {
     753            // use the center of the center pixel of the image
     754            float xc = (int)((0.5 + 0.225*x)*readout->image->numCols) + readout->image->col0 + 0.5;
     755            float yc = (int)((0.5 + 0.225*y)*readout->image->numRows) + readout->image->row0 + 0.5;
     756
     757            // assign the x and y coords to the image center
     758            // create an object with center intensity of 1000
     759            modelRef->params->data.F32[PM_PAR_SKY] = 0;
     760            modelRef->params->data.F32[PM_PAR_I0] = 1000;
     761            modelRef->params->data.F32[PM_PAR_XPOS] = xc;
     762            modelRef->params->data.F32[PM_PAR_YPOS] = yc;
     763
     764            // create modelPSF from this model
     765            pmModel *model = pmModelFromPSF (modelRef, psf);
     766            if (!model) continue;
     767
     768            // place the reference object in the image center
     769            // no need to mask the source here
     770            // XXX should we measure this for the analytical model only or the full model?
     771            pmModelAddWithOffset (psfMosaic, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_CENTER, 0, x*DX, y*DY);
     772            pmModelAddWithOffset (funMosaic, NULL, model, PM_MODEL_OP_FUNC | PM_MODEL_OP_CENTER, 0, x*DX, y*DY);
     773            pmModelAddWithOffset (resMosaic, NULL, model, PM_MODEL_OP_RES0 | PM_MODEL_OP_RES1 | PM_MODEL_OP_CENTER, 0, x*DX, y*DY);
     774            psFree (model);
     775        }
    776776    }
    777777
    … …  
    792792    fprintf (stdout, "[c]ontinue? ");
    793793    if (!fgets(key, 8, stdin)) {
    794         psWarning("Unable to read option");
     794        psWarning("Unable to read option");
    795795    }
    796796    return true;
    … …  
    805805    if (kapa2 == -1) {
    806806        kapa2 = KapaOpenNamedSocket ("kapa", "psphot:psfstars");
    807         if (kapa2 == -1) {
    808             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    809             isVisual = false;
    810             return false;
    811         }
    812     } 
     807        if (kapa2 == -1) {
     808            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     809            isVisual = false;
     810            return false;
     811        }
     812    }
    813813
    814814    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    827827
    828828    // counters to track the size of the image and area used in a row
    829     int dX = 0;                         // starting corner of next box
    830     int dY = 0;                         // height of row so far
    831     int NX = 20*DX;                     // full width of output image
    832     int NY = 0;                         // total height of output image
     829    int dX = 0;                         // starting corner of next box
     830    int dY = 0;                         // height of row so far
     831    int NX = 20*DX;                     // full width of output image
     832    int NY = 0;                         // total height of output image
    833833
    834834    // first, examine the PSF stars:
    … …  
    838838        pmSource *source = sources->data[i];
    839839
    840         bool keep = false;
     840        bool keep = false;
    841841        keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR);
    842         if (!keep) continue;
    843 
    844         // how does this subimage get placed into the output image?
    845         // DX = source->pixels->numCols
    846         // DY = source->pixels->numRows
    847 
    848         if (dX + DX > NX) {
    849             // too wide for the rest of this row
    850             if (dX == 0) {
    851                 // alone on this row
    852                 NY += DY;
    853                 dX = 0;
    854                 dY = 0;
    855             } else {
    856                 // start the next row
    857                 NY += dY;
    858                 dX = DX;
    859                 dY = DY;
    860             }
    861         } else {
    862             // extend this row
    863             dX += DX;
    864             dY = PS_MAX (dY, DY);
    865         }
     842        if (!keep) continue;
     843
     844        // how does this subimage get placed into the output image?
     845        // DX = source->pixels->numCols
     846        // DY = source->pixels->numRows
     847
     848        if (dX + DX > NX) {
     849            // too wide for the rest of this row
     850            if (dX == 0) {
     851                // alone on this row
     852                NY += DY;
     853                dX = 0;
     854                dY = 0;
     855            } else {
     856                // start the next row
     857                NY += dY;
     858                dX = DX;
     859                dY = DY;
     860            }
     861        } else {
     862            // extend this row
     863            dX += DX;
     864            dY = PS_MAX (dY, DY);
     865        }
    866866    }
    867867    NY += DY;
    … …  
    873873    psImageInit (outsub, 0.0);
    874874
    875     int Xo = 0;                         // starting corner of next box
    876     int Yo = 0;                         // starting corner of next box
    877     dY = 0;                             // height of row so far
     875    int Xo = 0;                         // starting corner of next box
     876    int Yo = 0;                         // starting corner of next box
     877    dY = 0;                             // height of row so far
    878878
    879879    int nPSF = 0;
    … …  
    885885        pmSource *source = sources->data[i];
    886886
    887         bool keep = false;
     887        bool keep = false;
    888888        if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
    889             nPSF ++;
    890             keep = true;
    891         }
    892         if (!keep) continue;
    893 
    894         if (Xo + DX > NX) {
    895             // too wide for the rest of this row
    896             if (Xo == 0) {
    897                 // place source alone on this row
    898                 psphotAddWithTest (source, true, maskVal); // replace source if subtracted
    899                 psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
    900 
    901                 psphotSubWithTest (source, false, maskVal); // remove source (force)
    902                 psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
    903                 psphotSetState (source, false, maskVal); // reset source Add/Sub state to recorded
    904 
    905                 Yo += DY;
    906                 Xo = 0;
    907                 dY = 0;
    908             } else {
    909                 // start the next row
    910                 Yo += dY;
    911                 Xo = 0;
    912 
    913                 psphotAddWithTest (source, true, maskVal); // replace source if subtracted
    914                 psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
    915 
    916                 psphotSubWithTest (source, false, maskVal); // remove source (force)
    917                 psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
    918                 psphotSetState (source, false, maskVal); // replace source (has been subtracted)
    919 
    920                 Xo = DX;
    921                 dY = DY;
    922             }
    923         } else {
    924             // extend this row
    925             psphotAddWithTest (source, true, maskVal); // replace source if subtracted
    926             psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
    927 
    928             psphotSubWithTest (source, false, maskVal); // remove source (force)
    929             psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
    930             psphotSetState (source, false, maskVal); // replace source (has been subtracted)
    931 
    932             Xo += DX;
    933             dY = PS_MAX (dY, DY);
    934         }
     889            nPSF ++;
     890            keep = true;
     891        }
     892        if (!keep) continue;
     893
     894        if (Xo + DX > NX) {
     895            // too wide for the rest of this row
     896            if (Xo == 0) {
     897                // place source alone on this row
     898                psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     899                psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
     900
     901                psphotSubWithTest (source, false, maskVal); // remove source (force)
     902                psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
     903                psphotSetState (source, false, maskVal); // reset source Add/Sub state to recorded
     904
     905                Yo += DY;
     906                Xo = 0;
     907                dY = 0;
     908            } else {
     909                // start the next row
     910                Yo += dY;
     911                Xo = 0;
     912
     913                psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     914                psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
     915
     916                psphotSubWithTest (source, false, maskVal); // remove source (force)
     917                psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
     918                psphotSetState (source, false, maskVal); // replace source (has been subtracted)
     919
     920                Xo = DX;
     921                dY = DY;
     922            }
     923        } else {
     924            // extend this row
     925            psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     926            psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true);
     927
     928            psphotSubWithTest (source, false, maskVal); // remove source (force)
     929            psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true);
     930            psphotSetState (source, false, maskVal); // replace source (has been subtracted)
     931
     932            Xo += DX;
     933            dY = PS_MAX (dY, DY);
     934        }
    935935    }
    936936
    … …  
    943943    fprintf (stdout, "[c]ontinue? ");
    944944    if (!fgets(key, 8, stdin)) {
    945         psWarning("Unable to read option");
     945        psWarning("Unable to read option");
    946946    }
    947947
    … …  
    965965    if (kapa2 == -1) {
    966966        kapa2 = KapaOpenNamedSocket ("kapa", "psphot:images");
    967         if (kapa2 == -1) {
    968             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    969             isVisual = false;
    970             return false;
    971         }
    972     } 
     967        if (kapa2 == -1) {
     968            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     969            isVisual = false;
     970            return false;
     971        }
     972    }
    973973
    974974    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    987987
    988988    // counters to track the size of the image and area used in a row
    989     int dX = 0;                         // starting corner of next box
    990     int dY = 0;                         // height of row so far
    991     int NX = 10*DX;                     // full width of output image
    992     int NY = 0;                         // total height of output image
     989    int dX = 0;                         // starting corner of next box
     990    int dY = 0;                         // height of row so far
     991    int NX = 10*DX;                     // full width of output image
     992    int NY = 0;                         // total height of output image
    993993
    994994    // first, examine the PSF and SAT stars:
    … …  
    998998        pmSource *source = sources->data[i];
    999999
    1000         bool keep = false;
     1000        bool keep = false;
    10011001        keep |= (source->mode & PM_SOURCE_MODE_SATSTAR);
    1002         if (!keep) continue;
    1003 
    1004         // how does this subimage get placed into the output image?
    1005         // DX = source->pixels->numCols
    1006         // DY = source->pixels->numRows
    1007 
    1008         if (dX + DX > NX) {
    1009             // too wide for the rest of this row
    1010             if (dX == 0) {
    1011                 // alone on this row
    1012                 NY += DY;
    1013                 dX = 0;
    1014                 dY = 0;
    1015             } else {
    1016                 // start the next row
    1017                 NY += dY;
    1018                 dX = DX;
    1019                 dY = DY;
    1020             }
    1021         } else {
    1022             // extend this row
    1023             dX += DX;
    1024             dY = PS_MAX (dY, DY);
    1025         }
     1002        if (!keep) continue;
     1003
     1004        // how does this subimage get placed into the output image?
     1005        // DX = source->pixels->numCols
     1006        // DY = source->pixels->numRows
     1007
     1008        if (dX + DX > NX) {
     1009            // too wide for the rest of this row
     1010            if (dX == 0) {
     1011                // alone on this row
     1012                NY += DY;
     1013                dX = 0;
     1014                dY = 0;
     1015            } else {
     1016                // start the next row
     1017                NY += dY;
     1018                dX = DX;
     1019                dY = DY;
     1020            }
     1021        } else {
     1022            // extend this row
     1023            dX += DX;
     1024            dY = PS_MAX (dY, DY);
     1025        }
    10261026    }
    10271027    NY += DY;
    … …  
    10311031    psImageInit (outsat, 0.0);
    10321032
    1033     int Xo = 0;                         // starting corner of next box
    1034     int Yo = 0;                         // starting corner of next box
    1035     dY = 0;                             // height of row so far
     1033    int Xo = 0;                         // starting corner of next box
     1034    int Yo = 0;                         // starting corner of next box
     1035    dY = 0;                             // height of row so far
    10361036
    10371037    int nSAT = 0;
    … …  
    10431043        pmSource *source = sources->data[i];
    10441044
    1045         bool keep = false;
     1045        bool keep = false;
    10461046        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
    1047             nSAT ++;
    1048             keep = true;
    1049         }           
    1050         if (!keep) continue;
    1051 
    1052         if (Xo + DX > NX) {
    1053             // too wide for the rest of this row
    1054             if (Xo == 0) {
    1055                 // place source alone on this row
    1056                 psphotAddWithTest (source, true, maskVal); // replace source if subtracted
    1057                 psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
    1058                 psphotSetState (source, true, maskVal); // reset source Add/Sub state to recorded
    1059 
    1060                 Yo += DY;
    1061                 Xo = 0;
    1062                 dY = 0;
    1063             } else {
    1064                 // start the next row
    1065                 Yo += dY;
    1066                 Xo = 0;
    1067                 psphotAddWithTest (source, true, maskVal); // replace source if subtracted
    1068                 psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
    1069                 psphotSetState (source, true, maskVal); // replace source (has been subtracted)
    1070 
    1071                 Xo = DX;
    1072                 dY = DY;
    1073             }
    1074         } else {
    1075             // extend this row
    1076             psphotAddWithTest (source, true, maskVal); // replace source if subtracted
    1077             psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
    1078             psphotSetState (source, true, maskVal); // replace source (has been subtracted)
    1079 
    1080             Xo += DX;
    1081             dY = PS_MAX (dY, DY);
    1082         }
     1047            nSAT ++;
     1048            keep = true;
     1049        }
     1050        if (!keep) continue;
     1051
     1052        if (Xo + DX > NX) {
     1053            // too wide for the rest of this row
     1054            if (Xo == 0) {
     1055                // place source alone on this row
     1056                psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     1057                psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
     1058                psphotSetState (source, true, maskVal); // reset source Add/Sub state to recorded
     1059
     1060                Yo += DY;
     1061                Xo = 0;
     1062                dY = 0;
     1063            } else {
     1064                // start the next row
     1065                Yo += dY;
     1066                Xo = 0;
     1067                psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     1068                psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
     1069                psphotSetState (source, true, maskVal); // replace source (has been subtracted)
     1070
     1071                Xo = DX;
     1072                dY = DY;
     1073            }
     1074        } else {
     1075            // extend this row
     1076            psphotAddWithTest (source, true, maskVal); // replace source if subtracted
     1077            psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false);
     1078            psphotSetState (source, true, maskVal); // replace source (has been subtracted)
     1079
     1080            Xo += DX;
     1081            dY = PS_MAX (dY, DY);
     1082        }
    10831083    }
    10841084
    … …  
    10901090    fprintf (stdout, "[c]ontinue? ");
    10911091    if (!fgets(key, 8, stdin)) {
    1092         psWarning("Unable to read option");
     1092        psWarning("Unable to read option");
    10931093    }
    10941094
    … …  
    11171117    float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS] - source->pixels->row0;
    11181118    for (int iy = 0; iy < source->pixels->numRows; iy++) {
    1119         for (int ix = 0; ix < source->pixels->numCols; ix++) {
    1120             if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
    1121                 // rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
    1122                 rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ;
    1123                 Rb->data.F32[nb] = log10(rb->data.F32[nb]);
    1124                 fb->data.F32[nb] = log10(source->pixels->data.F32[iy][ix]);
    1125                 nb++;
    1126             } else {
    1127                 // rg->data.F32[ng] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
    1128                 rg->data.F32[ng] = hypot (ix - Xo, iy - Yo) ;
    1129                 Rg->data.F32[ng] = log10(rg->data.F32[ng]);
    1130                 fg->data.F32[ng] = log10(source->pixels->data.F32[iy][ix]);
    1131                 ng++;
    1132             }
    1133         }
    1134     }
    1135  
     1119        for (int ix = 0; ix < source->pixels->numCols; ix++) {
     1120            if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
     1121                // rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
     1122                rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ;
     1123                Rb->data.F32[nb] = log10(rb->data.F32[nb]);
     1124                fb->data.F32[nb] = log10(source->pixels->data.F32[iy][ix]);
     1125                nb++;
     1126            } else {
     1127                // rg->data.F32[ng] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
     1128                rg->data.F32[ng] = hypot (ix - Xo, iy - Yo) ;
     1129                Rg->data.F32[ng] = log10(rg->data.F32[ng]);
     1130                fg->data.F32[ng] = log10(source->pixels->data.F32[iy][ix]);
     1131                ng++;
     1132            }
     1133        }
     1134    }
     1135
    11361136    // reset source Add/Sub state to recorded
    1137     psphotSetState (source, state, maskVal); 
     1137    psphotSetState (source, state, maskVal);
    11381138
    11391139    KapaInitGraph (&graphdata);
    … …  
    11481148    graphdata.ymax = +5.05;
    11491149    KapaSetLimits (myKapa, &graphdata);
    1150  
     1150
    11511151    KapaSetFont (myKapa, "helvetica", 14);
    11521152    KapaBox (myKapa, &graphdata);
    11531153    KapaSendLabel (myKapa, "radius (pixels)", KAPA_LABEL_XM);
    11541154    KapaSendLabel (myKapa, "log flux (counts)", KAPA_LABEL_YM);
    1155                
     1155
    11561156    graphdata.color = KapaColorByName ("black");
    11571157    graphdata.ptype = 2;
    … …  
    11691169    KapaPlotVector (myKapa, nb, rb->data.F32, "x");
    11701170    KapaPlotVector (myKapa, nb, fb->data.F32, "y");
    1171  
     1171
    11721172    // ** loglog **
    11731173    KapaSelectSection (myKapa, "loglog");
    … …  
    11791179    graphdata.ymax = +5.05;
    11801180    KapaSetLimits (myKapa, &graphdata);
    1181  
     1181
    11821182    KapaSetFont (myKapa, "helvetica", 14);
    11831183    KapaBox (myKapa, &graphdata);
    11841184    KapaSendLabel (myKapa, "log radius (pixels)", KAPA_LABEL_XM);
    11851185    KapaSendLabel (myKapa, "log flux (counts)", KAPA_LABEL_YM);
    1186                
     1186
    11871187    graphdata.color = KapaColorByName ("black");
    11881188    graphdata.ptype = 2;
    … …  
    12121212bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources) {
    12131213
    1214     KapaSection section;  // put the positive profile in one and the residuals in another? 
     1214    KapaSection section;  // put the positive profile in one and the residuals in another?
    12151215
    12161216    if (!isVisual) return true;
    … …  
    12181218    if (kapa3 == -1) {
    12191219        kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots");
    1220         if (kapa3 == -1) {
    1221             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    1222             isVisual = false;
    1223             return false;
    1224         }
    1225     } 
     1220        if (kapa3 == -1) {
     1221            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     1222            isVisual = false;
     1223            return false;
     1224        }
     1225    }
    12261226
    12271227    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    12571257        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
    12581258
    1259         psphotVisualPlotRadialProfile (kapa3, source, maskVal);
    1260 
    1261         // pause and wait for user input:
    1262         // continue, save (provide name), ??
    1263         char key[10];
    1264         fprintf (stdout, "[e]rase and continue? [o]verplot and continue? [s]kip rest of stars? : ");
    1265         if (!fgets(key, 8, stdin)) {
    1266             psWarning("Unable to read option");
    1267         }
    1268         if (key[0] == 'e') {
    1269             KapaClearPlots (kapa3);
    1270         }
    1271         if (key[0] == 's') {
    1272             break;
    1273         }
    1274     }       
     1259        psphotVisualPlotRadialProfile (kapa3, source, maskVal);
     1260
     1261        // pause and wait for user input:
     1262        // continue, save (provide name), ??
     1263        char key[10];
     1264        fprintf (stdout, "[e]rase and continue? [o]verplot and continue? [s]kip rest of stars? : ");
     1265        if (!fgets(key, 8, stdin)) {
     1266            psWarning("Unable to read option");
     1267        }
     1268        if (key[0] == 'e') {
     1269            KapaClearPlots (kapa3);
     1270        }
     1271        if (key[0] == 's') {
     1272            break;
     1273        }
     1274    }
    12751275
    12761276    return true;
    … …  
    12821282    int NoverlayO, NOVERLAYO;
    12831283    KiiOverlay *overlayE, *overlayO;
    1284  
     1284
    12851285    psEllipseMoments emoments;
    12861286    psEllipseAxes axes;
    … …  
    12891289
    12901290    if (kapa == -1) {
    1291         fprintf (stderr, "kapa not opened, skipping\n");
    1292         return false;
     1291        fprintf (stderr, "kapa not opened, skipping\n");
     1292        return false;
    12931293    }
    12941294
    … …  
    13041304    for (int i = 0; i < sources->n; i++) {
    13051305
    1306         float Xo, Yo, Rmaj, Rmin, cs, sn;
    1307 
    1308         pmSource *source = sources->data[i];
    1309         if (source == NULL) continue;
    1310 
    1311         pmMoments *moments = source->moments;
    1312         if (0) {
    1313             emoments.x2 = moments->Mxx;
    1314             emoments.y2 = moments->Myy;
    1315             emoments.xy = moments->Mxy;
    1316             Xo = moments->Mx;
    1317             Yo = moments->My;
    1318 
    1319             axes = psEllipseMomentsToAxes (emoments, 20.0);
    1320             Rmaj = 2.0*axes.major;
    1321             Rmin = 2.0*axes.minor;
    1322             cs = cos(axes.theta);
    1323             sn = sin(axes.theta);
    1324         } else {
    1325             Rmaj = Rmin = 5.0;
    1326             cs = 1.0;
    1327             sn = 0.0;
    1328             Xo = source->peak->xf;
    1329             Yo = source->peak->yf;
    1330         }
    1331 
    1332         unsigned short int flagMask = 0x01;
    1333         for (int j = 0; j < 8; j++) {
    1334             if (source->mode & flagMask) {
    1335                 overlayE[NoverlayE].type = KII_OVERLAY_LINE;
    1336                 overlayE[NoverlayE].x = Xo;
    1337                 overlayE[NoverlayE].y = Yo;
    1338 
    1339                 float phi = j*M_PI/4.0;
    1340                 overlayE[NoverlayE].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn;
    1341                 overlayE[NoverlayE].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs;
    1342                 overlayE[NoverlayE].angle = 0;
    1343                 overlayE[NoverlayE].text = NULL;
    1344                 NoverlayE ++;
    1345                 CHECK_REALLOCATE (overlayE, KiiOverlay, NOVERLAYE, NoverlayE, 100);
    1346             }
    1347             flagMask <<= 1;
    1348 
    1349             if (source->mode & flagMask) {
    1350                 overlayO[NoverlayO].type = KII_OVERLAY_LINE;
    1351                 overlayO[NoverlayO].x = Xo + 1;
    1352                 overlayO[NoverlayO].y = Yo;
    1353 
    1354                 float phi = j*M_PI/4.0;
    1355                 overlayO[NoverlayO].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn;
    1356                 overlayO[NoverlayO].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs;
    1357                 overlayO[NoverlayO].angle = 0;
    1358                 overlayO[NoverlayO].text = NULL;
    1359                 NoverlayO ++;
    1360                 CHECK_REALLOCATE (overlayO, KiiOverlay, NOVERLAYO, NoverlayO, 100);
    1361             }
    1362             flagMask <<= 1;
    1363         }
     1306        float Xo, Yo, Rmaj, Rmin, cs, sn;
     1307
     1308        pmSource *source = sources->data[i];
     1309        if (source == NULL) continue;
     1310
     1311        pmMoments *moments = source->moments;
     1312        if (0) {
     1313            emoments.x2 = moments->Mxx;
     1314            emoments.y2 = moments->Myy;
     1315            emoments.xy = moments->Mxy;
     1316            Xo = moments->Mx;
     1317            Yo = moments->My;
     1318
     1319            axes = psEllipseMomentsToAxes (emoments, 20.0);
     1320            Rmaj = 2.0*axes.major;
     1321            Rmin = 2.0*axes.minor;
     1322            cs = cos(axes.theta);
     1323            sn = sin(axes.theta);
     1324        } else {
     1325            Rmaj = Rmin = 5.0;
     1326            cs = 1.0;
     1327            sn = 0.0;
     1328            Xo = source->peak->xf;
     1329            Yo = source->peak->yf;
     1330        }
     1331
     1332        unsigned short int flagMask = 0x01;
     1333        for (int j = 0; j < 8; j++) {
     1334            if (source->mode & flagMask) {
     1335                overlayE[NoverlayE].type = KII_OVERLAY_LINE;
     1336                overlayE[NoverlayE].x = Xo;
     1337                overlayE[NoverlayE].y = Yo;
     1338
     1339                float phi = j*M_PI/4.0;
     1340                overlayE[NoverlayE].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn;
     1341                overlayE[NoverlayE].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs;
     1342                overlayE[NoverlayE].angle = 0;
     1343                overlayE[NoverlayE].text = NULL;
     1344                NoverlayE ++;
     1345                CHECK_REALLOCATE (overlayE, KiiOverlay, NOVERLAYE, NoverlayE, 100);
     1346            }
     1347            flagMask <<= 1;
     1348
     1349            if (source->mode & flagMask) {
     1350                overlayO[NoverlayO].type = KII_OVERLAY_LINE;
     1351                overlayO[NoverlayO].x = Xo + 1;
     1352                overlayO[NoverlayO].y = Yo;
     1353
     1354                float phi = j*M_PI/4.0;
     1355                overlayO[NoverlayO].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn;
     1356                overlayO[NoverlayO].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs;
     1357                overlayO[NoverlayO].angle = 0;
     1358                overlayO[NoverlayO].text = NULL;
     1359                NoverlayO ++;
     1360                CHECK_REALLOCATE (overlayO, KiiOverlay, NOVERLAYO, NoverlayO, 100);
     1361            }
     1362            flagMask <<= 1;
     1363        }
    13641364    }
    13651365
    … …  
    13761376    fprintf (stdout, "[c]ontinue? ");
    13771377    if (!fgets(key, 8, stdin)) {
    1378         psWarning("Unable to read option");
     1378        psWarning("Unable to read option");
    13791379    }
    13801380
    … …  
    13901390
    13911391    if (kapa == -1) {
    1392         fprintf (stderr, "kapa not opened, skipping\n");
    1393         return false;
     1392        fprintf (stderr, "kapa not opened, skipping\n");
     1393        return false;
    13941394    }
    13951395
    … …  
    14021402    for (int i = 0; i < sources->n; i++) {
    14031403
    1404         pmSource *source = sources->data[i];
    1405         if (source == NULL) continue;
    1406 
    1407         if (!(source->mode & PM_SOURCE_MODE_CR_LIMIT)) continue;
    1408 
    1409         overlay[Noverlay].type = KII_OVERLAY_BOX;
    1410         overlay[Noverlay].x = source->peak->xf;
    1411         overlay[Noverlay].y = source->peak->yf;
    1412 
    1413         overlay[Noverlay].dx = 4;
    1414         overlay[Noverlay].dy = 4;
    1415         overlay[Noverlay].angle = 0;
    1416         overlay[Noverlay].text = NULL;
    1417         Noverlay ++;
    1418         CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     1404        pmSource *source = sources->data[i];
     1405        if (source == NULL) continue;
     1406
     1407        if (!(source->mode & PM_SOURCE_MODE_CR_LIMIT)) continue;
     1408
     1409        overlay[Noverlay].type = KII_OVERLAY_BOX;
     1410        overlay[Noverlay].x = source->peak->xf;
     1411        overlay[Noverlay].y = source->peak->yf;
     1412
     1413        overlay[Noverlay].dx = 4;
     1414        overlay[Noverlay].dy = 4;
     1415        overlay[Noverlay].angle = 0;
     1416        overlay[Noverlay].text = NULL;
     1417        Noverlay ++;
     1418        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    14191419    }
    14201420    KiiLoadOverlay (kapa, overlay, Noverlay, "red");
    … …  
    14241424    for (int i = 0; i < sources->n; i++) {
    14251425
    1426         pmSource *source = sources->data[i];
    1427         if (source == NULL) continue;
    1428 
    1429         // mark EXTs with yellow circles
    1430         if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
    1431 
    1432         overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
    1433         overlay[Noverlay].x = source->peak->xf;
    1434         overlay[Noverlay].y = source->peak->yf;
    1435 
    1436         overlay[Noverlay].dx = 10;
    1437         overlay[Noverlay].dy = 10;
    1438         overlay[Noverlay].angle = 0;
    1439         overlay[Noverlay].text = NULL;
    1440         Noverlay ++;
    1441         CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
     1426        pmSource *source = sources->data[i];
     1427        if (source == NULL) continue;
     1428
     1429        // mark EXTs with yellow circles
     1430        if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
     1431
     1432        overlay[Noverlay].type = KII_OVERLAY_CIRCLE;
     1433        overlay[Noverlay].x = source->peak->xf;
     1434        overlay[Noverlay].y = source->peak->yf;
     1435
     1436        overlay[Noverlay].dx = 10;
     1437        overlay[Noverlay].dy = 10;
     1438        overlay[Noverlay].angle = 0;
     1439        overlay[Noverlay].text = NULL;
     1440        Noverlay ++;
     1441        CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100);
    14421442    }
    14431443
    … …  
    14531453    fprintf (stdout, "[c]ontinue? ");
    14541454    if (!fgets(key, 8, stdin)) {
    1455         psWarning("Unable to read option");
     1455        psWarning("Unable to read option");
    14561456    }
    14571457
    … …  
    14681468    if (kapa3 == -1) {
    14691469        kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots");
    1470         if (kapa3 == -1) {
    1471             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    1472             isVisual = false;
    1473             return false;
    1474         }
    1475     } 
     1470        if (kapa3 == -1) {
     1471            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     1472            isVisual = false;
     1473            return false;
     1474        }
     1475    }
    14761476
    14771477    KapaClearPlots (kapa3);
    … …  
    15001500    for (int i = 0; i < sources->n; i++) {
    15011501        pmSource *source = sources->data[i];
    1502         if (!source) continue;
     1502        if (!source) continue;
    15031503        if (source->type != PM_SOURCE_TYPE_STAR) continue;
    1504         if (!isfinite (source->crNsigma)) continue;
     1504        if (!isfinite (source->crNsigma)) continue;
    15051505
    15061506        x->data.F32[n] = -2.5*log10(source->peak->flux);
    … …  
    15551555    for (int i = 0; i < sources->n; i++) {
    15561556        pmSource *source = sources->data[i];
    1557         if (!source) continue;
     1557        if (!source) continue;
    15581558        if (source->type != PM_SOURCE_TYPE_STAR) continue;
    1559         if (!isfinite (source->extNsigma)) continue;
     1559        if (!isfinite (source->extNsigma)) continue;
    15601560
    15611561        x->data.F32[n] = -2.5*log10(source->peak->flux);
    … …  
    15971597    fprintf (stdout, "[c]ontinue? ");
    15981598    if (!fgets(key, 8, stdin)) {
    1599         psWarning("Unable to read option");
     1599        psWarning("Unable to read option");
    16001600    }
    16011601    return true;
    … …  
    16081608    if (kapa == -1) {
    16091609        kapa = KapaOpenNamedSocket ("kapa", "psphot:images");
    1610         if (kapa == -1) {
    1611             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    1612             isVisual = false;
    1613             return false;
    1614         }
    1615     } 
     1610        if (kapa == -1) {
     1611            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     1612            isVisual = false;
     1613            return false;
     1614        }
     1615    }
    16161616
    16171617    psphotVisualScaleImage (kapa, readout->image, "resid", 1);
    … …  
    16221622    fprintf (stdout, "[c]ontinue? ");
    16231623    if (!fgets(key, 8, stdin)) {
    1624         psWarning("Unable to read option");
     1624        psWarning("Unable to read option");
    16251625    }
    16261626    return true;
    … …  
    16351635    if (kapa3 == -1) {
    16361636        kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots");
    1637         if (kapa3 == -1) {
    1638             fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    1639             isVisual = false;
    1640             return false;
    1641         }
    1642     } 
     1637        if (kapa3 == -1) {
     1638            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     1639            isVisual = false;
     1640            return false;
     1641        }
     1642    }
    16431643
    16441644    KapaClearPlots (kapa3);
    … …  
    16571657    for (int i = 0; i < sources->n; i++) {
    16581658        pmSource *source = sources->data[i];
    1659         if (!source) continue;
     1659        if (!source) continue;
    16601660        if (source->type != PM_SOURCE_TYPE_STAR) continue;
    1661         if (!isfinite (source->apMag)) continue;
    1662         if (!isfinite (source->psfMag)) continue;
     1661        if (!isfinite (source->apMag)) continue;
     1662        if (!isfinite (source->psfMag)) continue;
    16631663
    16641664        x->data.F32[n] = source->psfMag;
    … …  
    17051705    fprintf (stdout, "[c]ontinue? ");
    17061706    if (!fgets(key, 8, stdin)) {
    1707         psWarning("Unable to read option");
     1707        psWarning("Unable to read option");
    17081708    }
    17091709    return true;
Note: See TracChangeset for help on using the changeset viewer.