IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 19, 2013, 4:45:57 PM (13 years ago)
Author:
eugene
Message:

carry over changes from eam_branch/20130307

Location:
branches/eam_branches/ipp-20130419
Files:
34 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20130419

  • branches/eam_branches/ipp-20130419/psModules

  • branches/eam_branches/ipp-20130419/psModules/src/astrom/pmAstrometryUtils.c

    r24096 r35421  
    2323#include "pmAstrometryUtils.h"
    2424
     25// this is used by the test output block
     26static int Nout = 0;
     27
    2528// given a 2D transformation -- L(x,y),M(x,y) -- find the coordinates x,y
    2629// for which L,M = 0,0. tol is the allowed error on x,y.
    … …  
    3336    double Xo = (trans->y->coeff[0][0]*trans->x->coeff[0][1] - trans->x->coeff[0][0]*trans->y->coeff[0][1])/R;
    3437    double Yo = (trans->x->coeff[0][0]*trans->y->coeff[1][0] - trans->y->coeff[0][0]*trans->x->coeff[1][0])/R;
     38
     39    if (0) {
     40      // this is a test output block, not used in normal operations
     41      char filename[64];
     42      snprintf (filename, 64, "trans.%03d.md", Nout);
     43      FILE *f = fopen (filename, "w");
     44
     45      psMetadata *md = psMetadataAlloc();
     46
     47      psPolynomial2DtoMetadata (md, trans->x, "X");
     48      psPolynomial2DtoMetadata (md, trans->y, "Y");
     49      psMetadataConfigPrint (f, md);
     50      psFree (md);
     51      fclose (f);
     52      Nout ++;
     53    }
    3554
    3655    // iterate to actual solution: requires small non-linear terms
    … …  
    5271        for (int i = 0; (dPos > tol) && (i < 20); i++) {
    5372            // NOTE: order for Alpha is: [y][x]
     73            // original Alpha->data.F32[1][0] = psPolynomial2DEval (XdY, Xo, Yo);
     74            // original Alpha->data.F32[0][1] = psPolynomial2DEval (YdX, Xo, Yo);
    5475            Alpha->data.F32[0][0] = psPolynomial2DEval (XdX, Xo, Yo);
    55             Alpha->data.F32[1][0] = psPolynomial2DEval (XdY, Xo, Yo);
    56             Alpha->data.F32[0][1] = psPolynomial2DEval (YdX, Xo, Yo);
     76            Alpha->data.F32[0][1] = psPolynomial2DEval (XdY, Xo, Yo);
     77            Alpha->data.F32[1][0] = psPolynomial2DEval (YdX, Xo, Yo);
    5778            Alpha->data.F32[1][1] = psPolynomial2DEval (YdY, Xo, Yo);
    5879
    59             Beta->data.F32[0] = psPolynomial2DEval (trans->x, Xo, Yo);
    60             Beta->data.F32[1] = psPolynomial2DEval (trans->y, Xo, Yo);
     80            Beta->data.F32[0] = -1.0 * psPolynomial2DEval (trans->x, Xo, Yo);
     81            Beta->data.F32[1] = -1.0 * psPolynomial2DEval (trans->y, Xo, Yo);
    6182
    6283            if (!psMatrixGJSolve (Alpha, Beta)) {
    … …  
    7192            }
    7293
    73             Xo -= Beta->data.F32[0];
    74             Yo -= Beta->data.F32[1];
     94            Xo += Beta->data.F32[0];
     95            Yo += Beta->data.F32[1];
    7596            dPos = hypot(Beta->data.F32[0], Beta->data.F32[1]);
    76 
     97            psTrace ("psastro", 5, "i: %d, Xo,Yo : %f %f  dX, dY: %f %f dPos: %f\n", i, Xo, Yo, Beta->data.F32[0], Beta->data.F32[1], dPos);
    7798        }
    7899        psFree (Alpha);

  • branches/eam_branches/ipp-20130419/psModules/src/camera/pmFPAMaskWeight.c

    r35038 r35421  
    301301    psRegion *trimsec = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC"); // Trim section
    302302    if (!mdok || psRegionIsNaN(*trimsec)) {
    303         psError(PS_ERR_IO, true, "CELL.TRIMSEC is not set --- unable to set variance.\n");
    304         return false;
     303      // if trimsec is not defined, use the full image
     304      trimsec = psRegionAlloc(0,0,0,0);
     305      // psError(PS_ERR_IO, true, "CELL.TRIMSEC is not set --- unable to set variance.\n");
     306      // return false;
    305307    }
    306308

  • branches/eam_branches/ipp-20130419/psModules/src/camera/pmFPAfileDefine.c

    r34085 r35421  
    446446}
    447447
     448// given a filename, convert to UNIX namespace and read the PHU
     449psMetadata *readPHUfromFilename (char *filename, pmConfig *config) {
     450
     451    // Need to generate an FPA
     452    psString realName = pmConfigConvertFilename(filename, config, false, false);
     453    if (!realName) {
     454        psError(psErrorCodeLast(), false, "Failed to convert file name %s", filename);
     455        return NULL;
     456    }
     457
     458    // load the header of the first image
     459    // EXTWORD (fits->extword) is not relevant to the PHU
     460    psFits *fits = psFitsOpen(realName, "r"); // FITS file
     461    if (!fits) {
     462        psError(psErrorCodeLast(), false, "Failed to open file %s", realName);
     463        psFree(realName);
     464        return NULL;
     465    }
     466
     467    psMetadata *phu = psFitsReadHeader (NULL, fits); // Primary header
     468    if (!phu) {
     469        psError(psErrorCodeLast(), false, "Failed to read file header %s", realName);
     470        psFree(realName);
     471        return NULL;
     472    }
     473
     474    if (!psFitsClose(fits)) {
     475        psError(psErrorCodeLast(), false, "Failed to close file %s", realName);
     476        psFree(realName);
     477        psFree(phu);
     478        return NULL;
     479    }
     480
     481    return phu;
     482}
     483
     484// this this function wants to return:
     485// pmFPA, PHU, fileLevel, outConfig
     486// camera, cameraName, formatName
     487typedef struct {
     488    pmFPA *fpa;
     489    psMetadata *phu;
     490    psMetadata *format;
     491    pmFPALevel fileLevel;
     492    psString cameraName;
     493    psString formatName;
     494} pmFPAfromFilenameOutput;
     495
     496// for the given filename, read PHU and determine camera format; build an FPA for the file
     497bool pmFPAfromFilename (pmFPAfromFilenameOutput *output, pmConfig **outConfig, pmConfig *sysConfig, char *filename){
     498
     499    // Need to generate an FPA
     500    psMetadata *phu = readPHUfromFilename (filename, sysConfig);
     501    if (!phu) {
     502        psError(psErrorCodeLast(), false, "Failed to read PHU for %s", filename);
     503        return false;
     504    }
     505
     506    // if we expect the loaded FPA to differ in configuration from the current system configuration
     507    // generate an output config for this FPA
     508    pmConfig *config = NULL;
     509    if (outConfig) {
     510        config = pmConfigAlloc();
     511        config->user = psMemIncrRefCounter(sysConfig->user);
     512        config->system = psMemIncrRefCounter(sysConfig->system);
     513
     514        psFree (config->files);
     515        config->files = psMemIncrRefCounter(sysConfig->files);
     516        psFree (config->arguments);
     517        config->arguments = psMemIncrRefCounter(sysConfig->arguments);
     518
     519        *outConfig = config;
     520    } else {
     521        config = sysConfig;
     522    }
     523
     524    // values which are returned to calling function
     525    psString formatName = NULL; // Name of camera format
     526    psString cameraName = NULL; // Name of camera
     527    psMetadata *camera = NULL;  // Camera configuration
     528
     529    // Determine the current format from the header; determine camera if not specified already.
     530    psMetadata *format = pmConfigCameraFormatFromHeader(&camera, &cameraName, &formatName, config, phu, true);
     531    if (!format) {
     532        psError(psErrorCodeLast(), false, "Failed to determine camera format for %s", filename);
     533        psFree(camera);
     534        psFree(formatName);
     535        psFree(phu);
     536        return false;
     537    }
     538
     539    pmFPALevel fileLevel = pmFPAPHULevel(format);
     540    if (fileLevel == PM_FPA_LEVEL_NONE) {
     541        psError(PM_ERR_CONFIG, true, "Unable to determine file level for %s", filename);
     542        psFree(camera);
     543        psFree(formatName);
     544        psFree(phu);
     545        return false;
     546    }
     547
     548    // build the template fpa, set up the basic view
     549    // we supply the metaCamera name (if NULL, baseCamera name is used)
     550    pmFPA *fpa = pmFPAConstruct(camera, cameraName);
     551    psFree(camera);
     552
     553    if (!fpa) {
     554        psError(psErrorCodeLast(), false, "Failed to construct FPA from %s", filename);
     555        psFree(formatName);
     556        psFree(format);
     557        psFree(phu);
     558        return NULL;
     559    }
     560
     561    output->fpa = fpa;
     562    output->phu = phu;
     563    output->format = format;
     564    output->fileLevel = fileLevel;
     565    output->cameraName = cameraName;
     566    output->formatName = formatName;
     567
     568    return true;
     569
     570}
    448571
    449572/// Define a file from an array of filenames
    450 static pmFPAfile *fpaFileDefineFromArray(pmConfig *config, // Configuration
     573static pmFPAfile *fpaFileDefineFromArray(pmConfig **outConfig, // output configuration
     574                                         pmConfig *sysConfig, // global configuration
    451575                                         pmFPAfile *bind, // File to bind to, or NULL
    452576                                         const char *name, // Name of file
    … …  
    454578    )
    455579{
    456     PS_ASSERT_PTR_NON_NULL(config, NULL);
     580    PS_ASSERT_PTR_NON_NULL(sysConfig, NULL);
    457581    PS_ASSERT_STRING_NON_EMPTY(name, NULL);
    458582
    … …  
    463587    pmFPALevel fileLevel = PM_FPA_LEVEL_NONE; // Level for files
    464588    psMetadata *phu = NULL;             // Primary header
     589
    465590    if (bind) {
    466591        // Use the FPA we're binding to
    … …  
    468593        fileLevel = bind->fileLevel;
    469594    } else {
    470         // Need to generate an FPA
    471         psString realName = pmConfigConvertFilename(filenames->data[0], config, false, false);
    472         if (!realName) {
    473             psError(psErrorCodeLast(), false, "Failed to convert file name %s", (char *)filenames->data[0]);
    474             return NULL;
    475         }
    476 
    477         // load the header of the first image
    478         // EXTWORD (fits->extword) is not relevant to the PHU
    479         psFits *fits = psFitsOpen(realName, "r"); // FITS file
    480         if (!fits) {
    481             psError(psErrorCodeLast(), false, "Failed to open file %s", realName);
    482             psFree(realName);
    483             return NULL;
    484         }
    485         phu = psFitsReadHeader (NULL, fits); // Primary header
    486         if (!phu) {
    487             psError(psErrorCodeLast(), false, "Failed to read file header %s", realName);
    488             psFree(realName);
    489             return NULL;
    490         }
    491         if (!psFitsClose(fits)) {
    492             psError(psErrorCodeLast(), false, "Failed to close file %s", realName);
    493             psFree(realName);
    494             return NULL;
    495         }
    496 
    497         // Determine the current format from the header; determine camera if not specified already.
    498         psMetadata *camera = NULL;      // Camera configuration
    499         format = pmConfigCameraFormatFromHeader(&camera, &cameraName, &formatName, config, phu, true);
    500         if (!format) {
    501             psError(psErrorCodeLast(), false, "Failed to determine camera format for %s", realName);
    502             psFree(camera);
    503             psFree(formatName);
    504             psFree(realName);
    505             psFree(phu);
    506             return NULL;
    507         }
    508 
    509         fileLevel = pmFPAPHULevel(format);
    510         if (fileLevel == PM_FPA_LEVEL_NONE) {
    511             psError(PM_ERR_CONFIG, true, "Unable to determine file level for %s", realName);
    512             psFree(camera);
    513             psFree(formatName);
    514             psFree(realName);
    515             psFree(phu);
    516             return NULL;
    517         }
    518 
    519         // build the template fpa, set up the basic view
    520         // XXX do we want this to be the baseCamera name or the metaCamera name?
    521         fpa = pmFPAConstruct(camera, cameraName);
    522         psFree(camera);
    523         if (!fpa) {
    524             psError(psErrorCodeLast(), false, "Failed to construct FPA from %s", realName);
    525             psFree(formatName);
    526             psFree(realName);
    527             psFree(format);
    528             psFree(phu);
    529             return NULL;
    530         }
    531         psFree(realName);
    532     }
     595        pmFPAfromFilenameOutput output;
     596        if (!pmFPAfromFilename (&output, outConfig, sysConfig, filenames->data[0])) {
     597            return NULL;
     598        }
     599        fpa = output.fpa;
     600        phu = output.phu;
     601        format = output.format;
     602        fileLevel = output.fileLevel;
     603        cameraName = output.cameraName;
     604        formatName = output.formatName;
     605    }
     606
     607    pmConfig *config = outConfig ? *outConfig : sysConfig;
    533608
    534609    // load the given filerule (from config->camera) and bind it to the fpa
    … …  
    561636        // Check that the file corresponds to the same camera and format
    562637        if (!phu) {
    563             psString realName = pmConfigConvertFilename(filenames->data[i], config, false, false);
    564             if (!realName) {
    565                 psError(psErrorCodeLast(), false, "Failed to convert file name %s", (char*)filenames->data[i]);
     638            phu = readPHUfromFilename (filenames->data[i], config);
     639            if (!phu) {
     640                psError(psErrorCodeLast(), false, "Failed to read PHU for %s", (char *)filenames->data[i]);
    566641                return NULL;
    567642            }
    568             psFits *fits = psFitsOpen(realName, "r"); // FITS file
    569             if (!fits) {
    570                 psError(psErrorCodeLast(), false, "Failed to open file %s", realName);
    571                 psFree(realName);
    572                 return NULL;
    573             }
    574             phu = psFitsReadHeader(NULL, fits);
    575             if (!psFitsClose(fits)) {
    576                 psError(psErrorCodeLast(), false, "Failed to close file %s", realName);
    577                 psFree(realName);
    578                 return NULL;
    579             }
    580             if (!phu) {
    581                 psError(psErrorCodeLast(), false, "Failed to read file header %s", realName);
    582                 psFree(realName);
    583                 return NULL;
    584             }
    585             psFree(realName);
    586643        }
    587644
    … …  
    651708}
    652709
    653 
    654 pmFPAfile *pmFPAfileDefineFromArgs(bool *success, pmConfig *config,
    655                                    const char *filename, const char *argname)
     710// find the file associated with the argname & generate a pmFPAfile for it based on the filerule
     711pmFPAfile *pmFPAfileDefineFromArgs(bool *success, pmConfig *config, const char *filename, const char *argname)
    656712{
    657713    PS_ASSERT_PTR_NON_NULL(config, NULL);
    … …  
    676732    }
    677733
    678     pmFPAfile *file = fpaFileDefineFromArray(config, NULL, filename, filenames); // File of interest
     734    pmFPAfile *file = fpaFileDefineFromArray(NULL, config, NULL, filename, filenames); // File of interest
    679735
    680736    if (success) {
    … …  
    685741}
    686742
    687 pmFPAfile *pmFPAfileBindFromArgs(bool *success, pmFPAfile *input, pmConfig *config,
    688                                 const char *filename, const char *argname)
     743// find the file associated with the argname & bind it to the given pmFPAfile for it based on the filerule
     744pmFPAfile *pmFPAfileBindFromArgs(bool *success, pmFPAfile *input, pmConfig *config, const char *filename, const char *argname)
    689745{
    690746    PS_ASSERT_PTR_NON_NULL(input, NULL);
    … …  
    710766    }
    711767
    712     pmFPAfile *file = fpaFileDefineFromArray(config, input, filename, filenames); // File of interest
     768    pmFPAfile *file = fpaFileDefineFromArray(NULL, config, input, filename, filenames); // File of interest
    713769
    714770    if (success) {
    … …  
    719775}
    720776
    721 pmFPAfile *pmFPAfileDefineSingleFromArgs(bool *success, pmConfig *config, const char *filename,
    722                                         const char *argname, int entry)
     777// find the specific file associated with the argname & generate a pmFPAfile for it based on the filerule
     778pmFPAfile *pmFPAfileDefineSingleFromArgs(bool *success, pmConfig *config, const char *filename, const char *argname, int entry)
    723779{
    724780    PS_ASSERT_PTR_NON_NULL(config, NULL);
    … …  
    746802    psArray *single = psArrayAlloc(1);  // Array of single filename of interest
    747803    single->data[0] = psMemIncrRefCounter(filenames->data[entry]);
    748     pmFPAfile *file = fpaFileDefineFromArray(config, NULL, filename, single); // File of interest
     804    pmFPAfile *file = fpaFileDefineFromArray(NULL, config, NULL, filename, single); // File of interest
    749805    psFree(single);
    750806
    … …  
    756812}
    757813
     814// find the file in the config list & generate a pmFPAfile for it based on the filerule
    758815pmFPAfile *pmFPAfileDefineFromRun(bool *success, pmFPAfile *bind, pmConfig *config, const char *filename)
    759816{
    … …  
    769826    }
    770827
    771     pmFPAfile *file = fpaFileDefineFromArray(config, bind, filename, filenames); // File of interest
     828    pmFPAfile *file = fpaFileDefineFromArray(NULL, config, bind, filename, filenames); // File of interest
    772829    psFree(filenames);
    773830
    … …  
    779836}
    780837
     838// find the files in the config list & generate an array of pmFPAfiles for them based on the filerule
    781839psArray *pmFPAfileDefineMultipleFromRun(bool *success, psArray *bind, pmConfig *config, const char *filename)
    782840{
    … …  
    808866        dummy->data[0] = files->data[i];
    809867        pmFPAfile *bindFile = bind ? bind->data[i] : NULL; // File to which to bind
    810         files->data[i] = psMemIncrRefCounter(fpaFileDefineFromArray(config, bindFile, filename, dummy));
     868        files->data[i] = psMemIncrRefCounter(fpaFileDefineFromArray(NULL, config, bindFile, filename, dummy));
    811869        if (!files->data[i]) {
    812870            psError(psErrorCodeLast(), false, "Unable to define file %s %d", filename, i);
    … …  
    823881
    824882    return files;
     883}
     884
     885// find the file associated with the argname & generate a pmFPAfile for it based on the filerule
     886pmFPAfile *pmFPAfileDefineNewConfig(bool *success, pmConfig **outConfig, pmConfig *sysConfig, const char *filename, const char *argname)
     887{
     888    PS_ASSERT_PTR_NON_NULL(outConfig, NULL);
     889    PS_ASSERT_PTR_NON_NULL(sysConfig, NULL);
     890    PS_ASSERT_STRING_NON_EMPTY(filename, NULL);
     891    PS_ASSERT_STRING_NON_EMPTY(argname, NULL);
     892
     893    // Search the argument data for the named fileset (argname)
     894    bool status;                        // Status of MD lookup
     895    psArray *filenames = psMetadataLookupPtr(&status, sysConfig->arguments, argname); // Filenames for file
     896    if (!status) {
     897        if (success) {
     898            *success = true;
     899        }
     900        return NULL;
     901    }
     902    if (filenames->n == 0) {
     903        psError(PM_ERR_CONFIG, true, "No files in array in %s in arguments", argname);
     904        if (success) {
     905            *success = false;
     906        }
     907        return NULL;
     908    }
     909
     910    pmFPAfile *file = fpaFileDefineFromArray(outConfig, sysConfig, NULL, filename, filenames); // File of interest
     911
     912    if (success) {
     913        *success = file ? true : false;
     914    }
     915
     916    return file;
    825917}
    826918

  • branches/eam_branches/ipp-20130419/psModules/src/camera/pmFPAfileDefine.h

    r27657 r35421  
    8181    const char *filename                ///< Name of file
    8282    );
     83
     84// find the file associated with the argname & generate a pmFPAfile for it based on the filerule
     85pmFPAfile *pmFPAfileDefineNewConfig(
     86    bool *success,                      ///< Found files?
     87    pmConfig **outConfig,               ///< output configuration for this file
     88    pmConfig *sysConfig,                ///< existing system config info
     89    const char *filename,               ///< name of filerule
     90    const char *argname                 ///< argument entry
     91    );
     92
    8393
    8494// look for the given argname on the argument list.  find the give filename from the file rules

  • branches/eam_branches/ipp-20130419/psModules/src/objects/models/pmModel_DEV.c

    r34403 r35421  
    295295// make an initial guess for parameters
    296296// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
    297 bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
     297bool PM_MODEL_GUESS (pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal)
    298298{
    299299    psF32 *PAR  = model->params->data.F32;

  • branches/eam_branches/ipp-20130419/psModules/src/objects/models/pmModel_DEV.h

    r29004 r35421  
    55psF32 pmModelFunc_DEV(psVector *deriv, const psVector *params, const psVector *pixcoord);
    66bool pmModelLimits_DEV(psMinConstraintMode mode, int nParam, float *params, float *beta);
    7 bool pmModelGuess_DEV(pmModel *model, pmSource *source);
     7bool pmModelGuess_DEV(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
    88psF64 pmModelFlux_DEV(const psVector *params);
    99psF64 pmModelRadius_DEV(const psVector *params, psF64 flux);

  • branches/eam_branches/ipp-20130419/psModules/src/objects/models/pmModel_EXP.c

    r34403 r35421  
    287287// make an initial guess for parameters
    288288// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
    289 bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
     289bool PM_MODEL_GUESS (pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal)
    290290{
    291291    psF32 *PAR  = model->params->data.F32;

  • branches/eam_branches/ipp-20130419/psModules/src/objects/models/pmModel_EXP.h

    r29004 r35421  
    55psF32 pmModelFunc_EXP(psVector *deriv, const psVector *params, const psVector *pixcoord);
    66bool pmModelLimits_EXP(psMinConstraintMode mode, int nParam, float *params, float *beta);
    7 bool pmModelGuess_EXP(pmModel *model, pmSource *source);
     7bool pmModelGuess_EXP(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
    88psF64 pmModelFlux_EXP(const psVector *params);
    99psF64 pmModelRadius_EXP(const psVector *params, psF64 flux);

  • branches/eam_branches/ipp-20130419/psModules/src/objects/models/pmModel_GAUSS.c

    r34403 r35421  
    192192// make an initial guess for parameters
    193193// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
    194 bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
     194bool PM_MODEL_GUESS (pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal)
    195195{
    196196    psF32 *PAR  = model->params->data.F32;

  • branches/eam_branches/ipp-20130419/psModules/src/objects/models/pmModel_GAUSS.h

    r25738 r35421  
    55psF32 pmModelFunc_GAUSS(psVector *deriv, const psVector *params, const psVector *pixcoord);
    66bool pmModelLimits_GAUSS(psMinConstraintMode mode, int nParam, float *params, float *beta);
    7 bool pmModelGuess_GAUSS(pmModel *model, pmSource *source);
     7bool pmModelGuess_GAUSS(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
    88psF64 pmModelFlux_GAUSS(const psVector *params);
    99psF64 pmModelRadius_GAUSS(const psVector *params, psF64 flux);

  • branches/eam_branches/ipp-20130419/psModules/src/objects/models/pmModel_PGAUSS.c

    r34403 r35421  
    193193// make an initial guess for parameters
    194194// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
    195 bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
     195bool PM_MODEL_GUESS (pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal)
    196196{
    197197    psF32 *PAR  = model->params->data.F32;

  • branches/eam_branches/ipp-20130419/psModules/src/objects/models/pmModel_PGAUSS.h

    r25738 r35421  
    55psF32 pmModelFunc_PGAUSS(psVector *deriv, const psVector *params, const psVector *pixcoord);
    66bool pmModelLimits_PGAUSS(psMinConstraintMode mode, int nParam, float *params, float *beta);
    7 bool pmModelGuess_PGAUSS(pmModel *model, pmSource *source);
     7bool pmModelGuess_PGAUSS(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
    88psF64 pmModelFlux_PGAUSS(const psVector *params);
    99psF64 pmModelRadius_PGAUSS(const psVector *params, psF64 flux);

  • branches/eam_branches/ipp-20130419/psModules/src/objects/models/pmModel_PS1_V1.c

    r34403 r35421  
    212212// make an initial guess for parameters
    213213// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
    214 bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
     214bool PM_MODEL_GUESS (pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal)
    215215{
    216216    psF32 *PAR  = model->params->data.F32;

  • branches/eam_branches/ipp-20130419/psModules/src/objects/models/pmModel_PS1_V1.h

    r25738 r35421  
    55psF32 pmModelFunc_PS1_V1(psVector *deriv, const psVector *params, const psVector *pixcoord);
    66bool pmModelLimits_PS1_V1(psMinConstraintMode mode, int nParam, float *params, float *beta);
    7 bool pmModelGuess_PS1_V1(pmModel *model, pmSource *source);
     7bool pmModelGuess_PS1_V1(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
    88psF64 pmModelFlux_PS1_V1(const psVector *params);
    99psF64 pmModelRadius_PS1_V1(const psVector *params, psF64 flux);

  • branches/eam_branches/ipp-20130419/psModules/src/objects/models/pmModel_QGAUSS.c

    r34403 r35421  
    213213// make an initial guess for parameters
    214214// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
    215 bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
     215bool PM_MODEL_GUESS (pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal)
    216216{
    217217    psF32 *PAR  = model->params->data.F32;

  • branches/eam_branches/ipp-20130419/psModules/src/objects/models/pmModel_QGAUSS.h

    r25738 r35421  
    55psF32 pmModelFunc_QGAUSS(psVector *deriv, const psVector *params, const psVector *pixcoord);
    66bool pmModelLimits_QGAUSS(psMinConstraintMode mode, int nParam, float *params, float *beta);
    7 bool pmModelGuess_QGAUSS(pmModel *model, pmSource *source);
     7bool pmModelGuess_QGAUSS(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
    88psF64 pmModelFlux_QGAUSS(const psVector *params);
    99psF64 pmModelRadius_QGAUSS(const psVector *params, psF64 flux);

  • branches/eam_branches/ipp-20130419/psModules/src/objects/models/pmModel_RGAUSS.c

    r34403 r35421  
    202202// make an initial guess for parameters
    203203// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
    204 bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
     204bool PM_MODEL_GUESS (pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal)
    205205{
    206206    psF32 *PAR  = model->params->data.F32;

  • branches/eam_branches/ipp-20130419/psModules/src/objects/models/pmModel_RGAUSS.h

    r25738 r35421  
    55psF32 pmModelFunc_RGAUSS(psVector *deriv, const psVector *params, const psVector *pixcoord);
    66bool pmModelLimits_RGAUSS(psMinConstraintMode mode, int nParam, float *params, float *beta);
    7 bool pmModelGuess_RGAUSS(pmModel *model, pmSource *source);
     7bool pmModelGuess_RGAUSS(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
    88psF64 pmModelFlux_RGAUSS(const psVector *params);
    99psF64 pmModelRadius_RGAUSS(const psVector *params, psF64 flux);

  • branches/eam_branches/ipp-20130419/psModules/src/objects/models/pmModel_SERSIC.c

    r34403 r35421  
    307307// make an initial guess for parameters
    308308// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
    309 bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
     309bool PM_MODEL_GUESS (pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal)
    310310{
    311311    pmMoments *moments = source->moments;

  • branches/eam_branches/ipp-20130419/psModules/src/objects/models/pmModel_SERSIC.h

    r25738 r35421  
    55psF32 pmModelFunc_SERSIC(psVector *deriv, const psVector *params, const psVector *pixcoord);
    66bool pmModelLimits_SERSIC(psMinConstraintMode mode, int nParam, float *params, float *beta);
    7 bool pmModelGuess_SERSIC(pmModel *model, pmSource *source);
     7bool pmModelGuess_SERSIC(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
    88psF64 pmModelFlux_SERSIC(const psVector *params);
    99psF64 pmModelRadius_SERSIC(const psVector *params, psF64 flux);

  • branches/eam_branches/ipp-20130419/psModules/src/objects/models/pmModel_TRAIL.c

    r35215 r35421  
    168168        dPAR[PM_PAR_THETA]  = PAR[PM_PAR_I0] * dPdT;
    169169        dPAR[PM_PAR_SIGMA]  = 0;        // we don't actually allow this to vary, so we do not need to calculate it
     170
     171        for (int i = 0; i < 7; i++) {
     172          if (isnan(dPAR[i])) {
     173            fprintf (stderr, "*");
     174          }
     175        }
     176    }
     177    if (isnan(f)) {
     178      fprintf (stderr, "!");
    170179    }
    171180    return(f);
    … …  
    228237}
    229238
     239# define NA 21
     240# define NR 21
     241static float flux[NA][NR];
     242static float npix[NA][NR];
     243
     244bool pmTrailGetAngle (float *To, float *Io, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal, float sigma) {
     245
     246  float Xo, Yo;
     247  if (!pmModelSetPosition(&Xo, &Yo, source)) return false;
     248
     249  psImage *image = source->pixels;
     250  psImage *mask = source->maskObj;
     251  psF32 **imData = image->data.F32;
     252  psImageMaskType **mkData = mask->data.PS_TYPE_IMAGE_MASK_DATA;
     253
     254  // do a loop over the pixels, generating (dX,dY) dot (cos(theta),sin(theta))
     255  int dP;
     256  int dX = Xo - mask->col0;
     257  dP = mask->numCols - dX;
     258  int DX = PS_MAX(dP, dX);
     259  int NX = mask->numCols;
     260
     261  int dY = Yo - mask->row0;
     262  dP = mask->numRows - dY;
     263  int DY = PS_MAX(dP, dY);
     264  int NY = mask->numRows;
     265
     266  // just hard wire this for now...
     267  float radius = 10.0;
     268  float radius2 = PS_SQR(radius);
     269
     270  // we have an array of Angles x Radii
     271  float dT = M_PI / NA;
     272  for (int na = 0; na < NA; na++) {
     273    memset (flux[na], 0, NR*sizeof(float));
     274    memset (npix[na], 0, NR*sizeof(float));
     275  }
     276
     277  // we skip any pixels [real or virtual] outside of the specified radius (nominally the aperture radius)
     278  // ix and iy track pixels relative to the centroid
     279  for (int ix = -DX; ix < DX + 1; ix++) {
     280    if (ix > radius) continue;
     281    int mx = ix + dX;
     282    for (int iy = -DY; iy < DY + 1; iy++) {
     283      if (iy > radius) continue;
     284      if (ix*ix + iy*iy > radius2) continue;
     285      int my = iy + dY;
     286     
     287      // include count only the unmasked pixels within the image area
     288      if (mx < 0) continue;
     289      if (my < 0) continue;
     290      if (mx >= NX) continue;
     291      if (my >= NY) continue;
     292     
     293      // count pixels which are masked only with bad pixels
     294      if (mkData[my][mx] & maskVal)continue;
     295
     296      // we have defined NA to be 21
     297      int na = 0;
     298      for (float angle = 0.0; na < NA; angle += dT, na ++) {
     299
     300        // XXX optimization : pre-compute the angle sines and cosines
     301        float Rad = (ix * cos(angle)) + (iy * sin(angle));
     302        int nr = PS_MAX (PS_MIN (NR, Rad + 0.5*NR), 0);
     303
     304        flux[na][nr] += imData[my][mx];
     305        npix[na][nr] ++;
     306      }
     307    }
     308  }
     309
     310  // generate a psf model (with integral of 1.0)
     311  float Ao = 1.0 / sqrt(2*M_PI) / sigma;
     312  float psf[NR];
     313  for (int nr = 0; nr < NR; nr++) {
     314    psf[nr] = Ao*exp(-0.5*PS_SQR((nr - 0.5*NR)/sigma));
     315  }
     316
     317  float fangle[NA];
     318  for (int na = 0; na < NA; na++) {
     319    float fpsf = 0.0;
     320    for (int nr = 0; nr < NR; nr++) {
     321      fpsf += psf[nr]*flux[na][nr];
     322    }
     323    fangle[na] = fpsf;
     324    // fprintf (stderr, "fpsf: %f, theta = %f\n", fpsf, PS_DEG_RAD*dT*na);
     325  }
     326
     327  float peak = fangle[0];
     328  int pbin = 0;
     329  for (int na = 0; na < NA; na++) {
     330    if (fangle[na] > peak) {
     331      peak = fangle[na];
     332      pbin = na;
     333    }
     334  }
     335
     336  // result is peak, pbin -> angle
     337  *To = dT * pbin;
     338  *Io = peak * Ao;
     339 
     340  return true;
     341}
     342
    230343// make an initial guess for parameters
    231344// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
    232 bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
     345bool PM_MODEL_GUESS (pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal)
    233346{
    234347    psF32 *PAR  = model->params->data.F32;
    … …  
    260373    //else { size = psfAxes.major; }
    261374
     375    float theta, peak;
     376    pmTrailGetAngle (&theta, &peak, source, maskVal, markVal, psfAxes.major);
     377
    262378    // axes.major is a sigma in the major direction; scale to
    263379    PAR[PM_PAR_LENGTH] = 1.5*2.35*size; // a tophat of length L has L = 1.5 * 2.35 * sigma
    … …  
    266382
    267383    // set the model normalization
    268     if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
    269       return false;
    270     }
     384    // if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
     385    //   return false;
     386    // }
     387    PAR[PM_PAR_I0] = peak;
    271388
    272389    // set the model position

  • branches/eam_branches/ipp-20130419/psModules/src/objects/models/pmModel_TRAIL.h

    r34259 r35421  
    55psF32 pmModelFunc_TRAIL(psVector *deriv, const psVector *params, const psVector *pixcoord);
    66bool pmModelLimits_TRAIL(psMinConstraintMode mode, int nParam, float *params, float *beta);
    7 bool pmModelGuess_TRAIL(pmModel *model, pmSource *source);
     7bool pmModelGuess_TRAIL(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
    88psF64 pmModelFlux_TRAIL(const psVector *params);
    99psF64 pmModelRadius_TRAIL(const psVector *params, psF64 flux);

  • branches/eam_branches/ipp-20130419/psModules/src/objects/pmModelFuncs.h

    r34259 r35421  
    107107//  This function provides the model guess parameters based on the details of
    108108//  the given source.
    109 typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source);
     109typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
    110110
    111111//  This function constructs the PSF model for the given source based on the

  • branches/eam_branches/ipp-20130419/psModules/src/objects/pmPSF.c

    r34403 r35421  
    350350        shape.sx  = modelPar[PM_PAR_SXX];
    351351        shape.sy  = modelPar[PM_PAR_SYY];
    352         shape.sxy = modelPar[PM_PAR_SXY];
     352        shape.sxy = modelPar[PM_PAR_SXY] / 2.0;
     353        // XXX I *think* dividing by 2.0 is the right direction, but this
     354        // needs to be checked with a real test
    353355    } else {
    354356        shape.sx  = modelPar[PM_PAR_SXX] / M_SQRT2;
    … …  
    392394        modelPar[PM_PAR_SXX] = shape.sx;
    393395        modelPar[PM_PAR_SYY] = shape.sy;
    394         modelPar[PM_PAR_SXY] = shape.sxy;
     396        modelPar[PM_PAR_SXY] = shape.sxy; // XXX NEED factor of 2 here for correct angle conversion
    395397    } else {
    396398        modelPar[PM_PAR_SXX] = shape.sx * M_SQRT2;

  • branches/eam_branches/ipp-20130419/psModules/src/objects/pmPSFtryFitEXT.c

    r34403 r35421  
    8383        }
    8484
    85         source->modelEXT = pmSourceModelGuess (source, options->type);
     85        source->modelEXT = pmSourceModelGuess (source, options->type, maskVal, markVal);
    8686        if (source->modelEXT == NULL) {
    8787            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y);

  • branches/eam_branches/ipp-20130419/psModules/src/objects/pmSourceFitModel.c

    r34403 r35421  
    7575    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    7676    PS_ASSERT_PTR_NON_NULL(source->maskObj, false);
    77     PS_ASSERT_PTR_NON_NULL(source->variance, false);
     77
     78    // XXX if variance is NULL, use pixels instead (wrong, but not badly wrong)
     79    // PS_ASSERT_PTR_NON_NULL(source->variance, false);
    7880
    7981    psBool fitStatus = true;
    … …  
    9395    float Xo = model->params->data.F32[PM_PAR_XPOS];
    9496    float Yo = model->params->data.F32[PM_PAR_YPOS];
     97
     98    // if variance is NULL, we pretend pixels == variance
     99    float **vWgt = source->variance ? source->variance->data.F32 : source->pixels->data.F32;
    95100
    96101    // fill in the coordinate and value entries
    … …  
    103108            }
    104109            // skip zero-variance points
    105             if (source->variance->data.F32[i][j] == 0) {
     110            if (vWgt[i][j] == 0) {
    106111                continue;
    107112            }
    108113            // skip nan values in image
    109114            if (!isfinite(source->pixels->data.F32[i][j])) {
     115                fprintf (stderr, "WARNING: unmasked nan in image : %x vs %x\n", source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[i][j], maskVal);
    110116                continue;
    111117            }
    112118
    113119            // skip nan values in image
    114             if (!isfinite(source->variance->data.F32[i][j])) {
    115                 fprintf (stderr, "impossible! %x vs %x\n", source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[i][j], maskVal);
     120            if (!isfinite(vWgt[i][j])) {
     121                fprintf (stderr, "WARNING: unmasked nan in variance : %x vs %x\n", source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[i][j], maskVal);
    116122                continue;
    117123            }
    … …  
    140146            // source sky variance
    141147            if (options->poissonErrors) {
    142                 yErr->data.F32[nPix] = 1.0 / source->variance->data.F32[i][j];
     148                yErr->data.F32[nPix] = 1.0 / vWgt[i][j];
    143149            } else {
    144150                yErr->data.F32[nPix] = 1.0 / options->weight;

  • branches/eam_branches/ipp-20130419/psModules/src/objects/pmSourceFitPCM.c

    r34403 r35421  
    145145bool pmSourceModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal) {
    146146
    147     if (!pcm->modelConv->modelGuess(pcm->modelConv, source)) {
     147  if (!pcm->modelConv->modelGuess(pcm->modelConv, source, maskVal, markVal)) {
    148148        return false;
    149149    }

  • branches/eam_branches/ipp-20130419/psModules/src/objects/pmSourceIO.c

    r34720 r35421  
    980980            return false;
    981981        }
     982
    982983        // if this is not TRUE, the output files only contain the psf measurements.
    983984        bool XSRC_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_ANALYSIS");
    … …  
    10051006            hdu->header = psFitsReadHeader (NULL, file->fits);
    10061007        }
     1008
     1009        // "WCS" is CMF without detected objects
     1010        if (file->type == PM_FPA_FILE_WCS) {
     1011          psTrace("psModules.objects", 6, "read CMF table from %s : %s : %s", file->filename, headname, dataname);
     1012          psFree (headname);
     1013          psFree (dataname);
     1014          psFree (deteffname);
     1015          break;
     1016        }
     1017
     1018        // EXTDATA is the PSF data associated with this image header
    10071019
    10081020        // we need to find the corresponding table EXTNAME.

  • branches/eam_branches/ipp-20130419/psModules/src/objects/pmSourceMoments.c

    r34403 r35421  
    151151
    152152        float *vPix = source->pixels->data.F32[row];
    153         float *vWgt = source->variance->data.F32[row];
     153        float *vWgt = source->variance ? source->variance->data.F32[row] : source->pixels->data.F32[row];
    154154
    155155        psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    … …  
    266266
    267267    float **vPix = source->pixels->data.F32;
    268     float **vWgt = source->variance->data.F32;
    269     psImageMaskType  **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
     268    float **vWgt = source->variance ? source->variance->data.F32 : source->pixels->data.F32;
     269    psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
    270270
    271271    float RF = 0.0;
    … …  
    467467
    468468        float *vPix = source->pixels->data.F32[row];
    469         float *vWgt = source->variance->data.F32[row];
     469        float *vWgt = source->variance ? source->variance->data.F32[row] : source->pixels->data.F32[row];
    470470
    471471        psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];

  • branches/eam_branches/ipp-20130419/psModules/src/objects/pmSourceSky.c

    r34403 r35421  
    122122    psTrace("psModules.objects", 10, "---- %s() begin ----\n", __func__);
    123123    PS_ASSERT_PTR_NON_NULL(source, false);
    124     PS_ASSERT_IMAGE_NON_NULL(source->variance, false);
    125124    PS_ASSERT_IMAGE_NON_NULL(source->maskObj, false);
    126125    PS_ASSERT_PTR_NON_NULL(source->peak, false);
    127126    PS_ASSERT_INT_POSITIVE(Radius, false);
     127
     128    if (source->variance == NULL) {
     129      // XXX this is needed so psphotCheckRadiusPSF has a real value
     130      source->moments->dSky = 1.0;
     131      return true;
     132    }
    128133
    129134    // maskVal is used to test for rejected pixels, and must include markVal

  • branches/eam_branches/ipp-20130419/psModules/src/objects/pmSourceUtils.c

    r34403 r35421  
    4646
    4747/******************************************************************************
    48     pmSourceModelGuess(source, model): This function allocates a new
    49     pmModel structure based on the given modelType specified in the argument list.
    50     The corresponding pmModelGuess function is returned, and used to
    51     supply the values of the params array in the pmModel structure.
     48    pmSourceModelGuess(source, model, maskVal, markVal): This function allocates a new
     49    pmModel structure based on the given modelType specified in the argument list.  The
     50    corresponding pmModelGuess function is returned, and used to supply the values of the
     51    params array in the pmModel structure.
    5252 
    5353    XXX: Many parameters are based on the src->moments structure, which is in
    … …  
    5555    functions will be in image, not subImage coords.  Remember this.
    5656*****************************************************************************/
    57 pmModel *pmSourceModelGuess(pmSource *source,
    58                             pmModelType modelType)
     57pmModel *pmSourceModelGuess(pmSource *source, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal)
    5958{
    6059    psTrace("psModules.objects", 10, "---- %s() begin ----\n", __func__);
    … …  
    6564    pmModel *model = pmModelAlloc(modelType);
    6665
    67     if (!model->modelGuess(model, source)) {
     66    if (!model->modelGuess(model, source, maskVal, markVal)) {
    6867        psFree (model);
    6968        return NULL;

  • branches/eam_branches/ipp-20130419/psModules/src/objects/pmSourceUtils.h

    r31451 r35421  
    2727pmModel *pmSourceModelGuess(
    2828    pmSource *source,   ///< The input pmSource
    29     pmModelType model   ///< The type of model to be created.
     29    pmModelType model,   ///< The type of model to be created.
     30    psImageMaskType maskVal,
     31    psImageMaskType markVal
    3032);
    3133
Note: See TracChangeset for help on using the changeset viewer.