IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 34053


Ignore:
Timestamp:
Jun 22, 2012, 3:33:32 PM (14 years ago)
Author:
eugene
Message:

read the time table in psTime.c; updates to the variance model tests; better trace options for output masks; convert to new PM_SOURCE_PHOTFIT version of weighting

Location:
branches/eam_branches/ipp-20120601
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20120601/Ohana/src/libkapa/src/DrawRotString.c

    r33982 r34053  
    7777    /* check for special characters */
    7878    if (!code && !protect) {
     79      /* subscript, starts with _ */
    7980      if (N == 94) {
    8081        SetRotFont (currentname, (int)(0.8*currentsize));
    … …  
    8384        continue;
    8485      }
     86      /* superscript, starts with ^ */
    8587      if (N == 95) {
    8688        SetRotFont (currentname, (int)(0.8*currentsize));
    … …  
    8991        continue;
    9092      }
     93      /* normal script, starts with | */
    9194      if (N == 124) {
    9295        SetRotFont (currentname, basesize);

  • branches/eam_branches/ipp-20120601/psLib/src/astro/psTime.c

    r34045 r34053  
    389389        if (i < numTables) {
    390390            table = psLookupTableAlloc(fullTableName, (const char*)tableFormat, tablesIndex->data.S32[i]);
    391             // XXX unused numLines = psLookupTableRead(table);
     391            psLookupTableRead(table);
    392392        } else {
    393393            psError(PS_ERR_BAD_PARAMETER_VALUE, no_problem,

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

    r34049 r34053  
    2929# define PM_MODEL_FUNCS_H
    3030
    31 # define HAVE_MODEL_VAR 1
    32 
    3331/// @addtogroup Objects Object Detection / Analysis Functions
    3432/// @{
    … …  
    5755    PM_MODEL_OP_NORM     = 0x20,
    5856    PM_MODEL_OP_NOISE    = 0x40,
    59 # if (HAVE_MODEL_VAR)
    6057    PM_MODEL_OP_MODELVAR = 0x80,
    61 # endif
    6258} pmModelOpMode;
    6359

  • branches/eam_branches/ipp-20120601/psModules/src/objects/pmModelUtils.c

    r33417 r34053  
    4343
    4444#include "pmErrorCodes.h"
     45
     46// XX static bool useModelVar = false;
     47// XX
     48// XX void pmModelSetModelVarOption (bool option) {
     49// XX   useModelVar = option;
     50// XX }
     51// XX bool pmModelGetModelVarOption (void) {
     52// XX   return useModelVar;
     53// XX }
    4554
    4655/*****************************************************************************

  • branches/eam_branches/ipp-20120601/psModules/src/objects/pmModelUtils.h

    r31153 r34053  
    4646bool pmModelSetShape (float *Sxx, float *Sxy, float *Syy, pmMoments *moments);
    4747
     48// XXX void pmModelSetModelVarOption (bool option);
     49// XXX bool pmModelGetModelVarOption (void);
     50
    4851/// @}
    4952# endif /* PM_MODEL_UTILS_H */

  • branches/eam_branches/ipp-20120601/psModules/src/objects/pmSource.c

    r34049 r34053  
    5151    psFree(tmp->pixels);
    5252    psFree(tmp->variance);
    53 # if (HAVE_MODEL_VAR)
    5453    psFree(tmp->modelVar);
    55 # endif
    5654    psFree(tmp->maskObj);
    5755    psFree(tmp->maskView);
    … …  
    8078    psFree (source->pixels);
    8179    psFree (source->variance);
    82 # if (HAVE_MODEL_VAR)
    8380    psFree (source->modelVar);
    84 # endif
    8581    psFree (source->maskObj);
    8682    psFree (source->maskView);
    … …  
    9086    source->pixels = NULL;
    9187    source->variance = NULL;
    92 # if (HAVE_MODEL_VAR)
    9388    source->modelVar = NULL;
    94 # endif
    9589    source->maskObj = NULL;
    9690    source->maskView = NULL;
    … …  
    122116    source->pixels = NULL;
    123117    source->variance = NULL;
    124 # if (HAVE_MODEL_VAR)
    125118    source->modelVar = NULL;
    126 # endif
    127119    source->maskObj = NULL;
    128120    source->maskView = NULL;
    … …  
    212204    source->pixels   = in->pixels   ? psImageCopyView(NULL, in->pixels)   : NULL;
    213205    source->variance = in->variance ? psImageCopyView(NULL, in->variance) : NULL;
    214 # if (HAVE_MODEL_VAR)
    215206    source->modelVar = NULL;
    216 # endif
    217207    source->maskView = in->maskView ? psImageCopyView(NULL, in->maskView) : NULL;
    218208
    … …  
    10441034    bool addNoise = mode & PM_MODEL_OP_NOISE;
    10451035
    1046 # if (HAVE_MODEL_VAR)
    10471036    bool addModelVar = mode & PM_MODEL_OP_MODELVAR;
    10481037    if (addModelVar) psAssert (source->modelVar, "programming error");
    1049 # endif
    10501038
    10511039    // require the use of pmModelAddWithOffset if we are adding noise (because the model size and norm are rescaled)
    … …  
    10701058        }
    10711059
    1072 # if (HAVE_MODEL_VAR)
    1073         psF32 **target = addModelVar ? source->modelVar->data.F32 : source->pixels->data.F32;
    1074 # else
    1075         psF32 **target = source->pixels->data.F32;
    1076 # endif
     1060        psF32 **target = NULL;
     1061        target = addModelVar ? source->modelVar->data.F32 : source->pixels->data.F32;
    10771062
    10781063        for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
    … …  
    10891074            }
    10901075        }
    1091 # if (HAVE_MODEL_VAR)
     1076        // do not change the flag here if we are adding/subtracting from modelVar
    10921077        if (!addModelVar) {
    1093             if (add) {
    1094                 source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
    1095             } else {
    1096                 source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
    1097             }
     1078          if (add) {
     1079            source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
     1080          } else {
     1081            source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
     1082          }
    10981083        }
    1099 # else
    1100         if (add) {
    1101             source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
    1102         } else {
    1103             source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
    1104         }
    1105 # endif
    11061084        return true;
    11071085    }
    … …  
    11111089        target = source->variance;
    11121090    }
    1113 # if (HAVE_MODEL_VAR)
    11141091    if (addModelVar) {
    1115         target = source->modelVar;
    1116     }
    1117 # endif
     1092      target = source->modelVar;
     1093    }
    11181094
    11191095    if (add) {
    11201096        status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
    1121 # if (HAVE_MODEL_VAR)
    11221097        if (!addNoise && !addModelVar) source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
    1123 # else
    1124         source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
    1125 # endif
    11261098    } else {
    11271099        status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
    1128 # if (HAVE_MODEL_VAR)
    11291100        if (!addNoise && !addModelVar) source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
    1130 # else
    1131         source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
    1132 # endif
    11331101    }
    11341102    if (!status) {

  • branches/eam_branches/ipp-20120601/psModules/src/objects/pmSource.h

    r34049 r34053  
    7272    psImage *pixels;                    ///< Rectangular region including object pixels.
    7373    psImage *variance;                  ///< Image variance.
    74 # if (HAVE_MODEL_VAR)
    7574    psImage *modelVar;                  ///< variance based on current models
    76 # endif
    7775    psImage *maskObj;                   ///< unique mask for this object which marks included pixels associated with objects.
    7876    psImage *maskView;                  ///< view into global image mask for this object region

  • branches/eam_branches/ipp-20120601/psModules/src/objects/pmSourcePhotometry.c

    r34049 r34053  
    899899}
    900900
    901 # if (HAVE_MODEL_VAR)
    902901double pmSourceModelWeight(const pmSource *Mi, int term, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
    903 # else
    904 double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
    905 # endif
    906902{
    907903    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    908 # if (HAVE_MODEL_VAR)
    909     double flux = 0, wt = 1.0, factor = 0;
    910 # else
    911     double flux = 0, wt = 0, factor = 0;
    912 # endif
     904    double flux = 0;
     905    double wt = 1.0;
     906    double factor = 0;
    913907
    914908    const psImage *Pi = Mi->modelFlux;
    915909    assert (Pi != NULL);
    916910
    917 # if (HAVE_MODEL_VAR)
    918911    const psImage *Wi = NULL;
    919912    switch (fitVarMode) {
    … …  
    921914        break;
    922915      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
     916      case PM_SOURCE_PHOTFIT_MODEL_SKY:
     917        Wi = Mi->variance;
     918        psAssert (Wi, "programming error");
     919        break;
     920      case PM_SOURCE_PHOTFIT_MODEL_VAR:
     921        Wi = Mi->modelVar;
     922        psAssert (Wi, "programming error");
     923        break;
     924      case PM_SOURCE_PHOTFIT_NONE:
     925        psAbort("programming error");
     926    }   
     927    const psImage *Ti = Mi->maskObj;
     928    assert (Ti != NULL);
     929
     930    for (int yi = 0; yi < Pi->numRows; yi++) {
     931        for (int xi = 0; xi < Pi->numCols; xi++) {
     932            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
     933                continue;
     934            if (fitVarMode != PM_SOURCE_PHOTFIT_CONST) {
     935                wt = covarFactor * Wi->data.F32[yi][xi];
     936                if (wt == 0) continue;
     937            }
     938            switch (term) {
     939              case 0:
     940                factor = 1;
     941                break;
     942              case 1:
     943                factor = xi + Pi->col0;
     944                break;
     945              case 2:
     946                factor = yi + Pi->row0;
     947                break;
     948              default:
     949                psAbort("invalid term for pmSourceWeight");
     950            }
     951
     952            // wt is 1.0 for CONST
     953            flux += (factor * Pi->data.F32[yi][xi]) / wt;
     954        }
     955    }
     956    return flux;
     957}
     958
     959double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
     960{
     961    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
     962    PS_ASSERT_PTR_NON_NULL(Mj, NAN);
     963    int Xs, Xe, Ys, Ye;
     964    int xi, xj, yi, yj;
     965    int xIs, xJs, yIs, yJs;
     966    int xIe, yIe;
     967    double flux;
     968    double wt = 1.0;
     969
     970    const psImage *Pi = Mi->modelFlux;
     971    assert (Pi != NULL);
     972    const psImage *Pj = Mj->modelFlux;
     973    assert (Pj != NULL);
     974
     975    const psImage *Wi = NULL;
     976    switch (fitVarMode) {
     977      case PM_SOURCE_PHOTFIT_CONST:
     978        break;
     979      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
     980      case PM_SOURCE_PHOTFIT_MODEL_SKY:
     981        Wi = Mi->variance;
     982        psAssert (Wi, "programming error");
     983        break;
     984      case PM_SOURCE_PHOTFIT_MODEL_VAR:
     985        Wi = Mi->modelVar;
     986        psAssert (Wi, "programming error");
     987        break;
     988      case PM_SOURCE_PHOTFIT_NONE:
     989        psAbort("programming error");
     990    }   
     991
     992    const psImage *Ti = Mi->maskObj;
     993    assert (Ti != NULL);
     994    const psImage *Tj = Mj->maskObj;
     995    assert (Tj != NULL);
     996
     997    Xs = PS_MAX (Pi->col0, Pj->col0);
     998    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);
     999
     1000    Ys = PS_MAX (Pi->row0, Pj->row0);
     1001    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);
     1002
     1003    xIs = Xs - Pi->col0;
     1004    xJs = Xs - Pj->col0;
     1005    yIs = Ys - Pi->row0;
     1006    yJs = Ys - Pj->row0;
     1007
     1008    xIe = Xe - Pi->col0;
     1009    yIe = Ye - Pi->row0;
     1010
     1011    // note that weight is addressing the same image pixels
     1012    flux = 0;
     1013    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
     1014        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
     1015            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
     1016                continue;
     1017            if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj] & maskVal)
     1018                continue;
     1019
     1020            float value = (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
     1021            switch (fitVarMode) {
     1022              case PM_SOURCE_PHOTFIT_CONST:
     1023                wt = 1.0;
     1024                break;
     1025              case PM_SOURCE_PHOTFIT_IMAGE_VAR:
     1026              case PM_SOURCE_PHOTFIT_MODEL_SKY:
     1027              case PM_SOURCE_PHOTFIT_MODEL_VAR:
     1028                wt = covarFactor * Wi->data.F32[yi][xi];
     1029                break;
     1030              case PM_SOURCE_PHOTFIT_NONE:
     1031                psAbort("programming error");
     1032            }
     1033            // skip pixels with nonsense weight values
     1034            if (wt <= 0) continue;
     1035
     1036            flux += value / wt;
     1037        }
     1038    }
     1039    return flux;
     1040}
     1041
     1042double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
     1043{
     1044    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
     1045    PS_ASSERT_PTR_NON_NULL(Mj, NAN);
     1046    int Xs, Xe, Ys, Ye;
     1047    int xi, xj, yi, yj;
     1048    int xIs, xJs, yIs, yJs;
     1049    int xIe, yIe;
     1050    double flux;
     1051    double wt = 1.0;
     1052
     1053    const psImage *Pi = Mi->pixels;
     1054    assert (Pi != NULL);
     1055    const psImage *Pj = Mj->modelFlux;
     1056    assert (Pj != NULL);
     1057
     1058    const psImage *Wi = NULL;
     1059    switch (fitVarMode) {
     1060      case PM_SOURCE_PHOTFIT_CONST:
     1061        break;
     1062      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
     1063      case PM_SOURCE_PHOTFIT_MODEL_SKY:
    9231064        Wi = Mi->variance;
    9241065        psAssert (Wi, "programming error");
    … …  
    9311072        psAbort("programming error");
    9321073    }   
    933 # else
    934     const psImage *Wi = Mi->variance;
    935     if (!unweighted_sum) {
    936         assert (Wi != NULL);
    937     }
    938 # endif
    939     const psImage *Ti = Mi->maskObj;
    940     assert (Ti != NULL);
    941 
    942     for (int yi = 0; yi < Pi->numRows; yi++) {
    943         for (int xi = 0; xi < Pi->numCols; xi++) {
    944             if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
    945                 continue;
    946 # if (HAVE_MODEL_VAR)
    947             if (fitVarMode != PM_SOURCE_PHOTFIT_CONST) {
    948                 wt = covarFactor * Wi->data.F32[yi][xi];
    949                 if (wt == 0) continue;
    950             }
    951 # else
    952             if (!unweighted_sum) {
    953                 wt = covarFactor * Wi->data.F32[yi][xi];
    954                 if (wt == 0)
    955                     continue;
    956             }
    957 # endif
    958 
    959             switch (term) {
    960               case 0:
    961                 factor = 1;
    962                 break;
    963               case 1:
    964                 factor = xi + Pi->col0;
    965                 break;
    966               case 2:
    967                 factor = yi + Pi->row0;
    968                 break;
    969               default:
    970                 psAbort("invalid term for pmSourceWeight");
    971             }
    972 
    973 # if (HAVE_MODEL_VAR)
    974             // wt is 1.0 for CONST
    975             flux += (factor * Pi->data.F32[yi][xi]) / wt;
    976 # else
    977             if (unweighted_sum) {
    978                 flux += (factor * Pi->data.F32[yi][xi]);
    979             } else {
    980                 flux += (factor * Pi->data.F32[yi][xi]) / wt;
    981             }
    982 # endif
    983         }
    984     }
    985     return flux;
    986 }
    987 
    988 # if (HAVE_MODEL_VAR)
    989 double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
    990 # else
    991 double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
    992 # endif
    993 {
    994     PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    995     PS_ASSERT_PTR_NON_NULL(Mj, NAN);
    996     int Xs, Xe, Ys, Ye;
    997     int xi, xj, yi, yj;
    998     int xIs, xJs, yIs, yJs;
    999     int xIe, yIe;
    1000 # if (HAVE_MODEL_VAR)
    1001     double flux;
    1002     double wt = 1.0;
    1003 # else
    1004     double flux, wt;
    1005 # endif
    1006 
    1007     const psImage *Pi = Mi->modelFlux;
    1008     assert (Pi != NULL);
    1009     const psImage *Pj = Mj->modelFlux;
    1010     assert (Pj != NULL);
    1011 
    1012 # if (HAVE_MODEL_VAR)
    1013     const psImage *Wi = NULL;
    1014     switch (fitVarMode) {
    1015       case PM_SOURCE_PHOTFIT_CONST:
    1016         break;
    1017       case PM_SOURCE_PHOTFIT_IMAGE_VAR:
    1018         Wi = Mi->variance;
    1019         psAssert (Wi, "programming error");
    1020         break;
    1021       case PM_SOURCE_PHOTFIT_MODEL_VAR:
    1022         Wi = Mi->modelVar;
    1023         psAssert (Wi, "programming error");
    1024         break;
    1025       case PM_SOURCE_PHOTFIT_NONE:
    1026         psAbort("programming error");
    1027     }   
    1028 # else
    1029     const psImage *Wi = Mi->variance;
    1030     if (!unweighted_sum) {
    1031         assert (Wi != NULL);
    1032     }
    1033 # endif
    10341074
    10351075    const psImage *Ti = Mi->maskObj;
    … …  
    10521092    yIe = Ye - Pi->row0;
    10531093
    1054     // note that weight is addressing the same image pixels
     1094    // note that weight is addressing the same image pixels,
    10551095    flux = 0;
    10561096    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
    … …  
    10611101                continue;
    10621102
    1063 # if (HAVE_MODEL_VAR)
    10641103            float value = (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
    10651104            switch (fitVarMode) {
    … …  
    10681107                break;
    10691108              case PM_SOURCE_PHOTFIT_IMAGE_VAR:
     1109              case PM_SOURCE_PHOTFIT_MODEL_SKY:
    10701110              case PM_SOURCE_PHOTFIT_MODEL_VAR:
    10711111                wt = covarFactor * Wi->data.F32[yi][xi];
    … …  
    10781118
    10791119            flux += value / wt;
    1080 # else
    1081             // XXX skip the nonsense weight pixels?
    1082             if (unweighted_sum) {
    1083                 flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
    1084             } else {
    1085                 wt = covarFactor * Wi->data.F32[yi][xi];
    1086                 if (wt > 0) {
    1087                     flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
    1088                 }
    1089             }
    1090 # endif
    10911120        }
    10921121    }
    10931122    return flux;
    10941123}
    1095 
    1096 # if (HAVE_MODEL_VAR)
    1097 double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
    1098 # else
    1099 double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
    1100 # endif
    1101 {
    1102     PS_ASSERT_PTR_NON_NULL(Mi, NAN);
    1103     PS_ASSERT_PTR_NON_NULL(Mj, NAN);
    1104     int Xs, Xe, Ys, Ye;
    1105     int xi, xj, yi, yj;
    1106     int xIs, xJs, yIs, yJs;
    1107     int xIe, yIe;
    1108 # if (HAVE_MODEL_VAR)
    1109     double flux;
    1110     double wt = 1.0;
    1111 # else
    1112     double flux, wt;
    1113 # endif
    1114 
    1115     const psImage *Pi = Mi->pixels;
    1116     assert (Pi != NULL);
    1117     const psImage *Pj = Mj->modelFlux;
    1118     assert (Pj != NULL);
    1119 
    1120 # if (HAVE_MODEL_VAR)
    1121     const psImage *Wi = NULL;
    1122     switch (fitVarMode) {
    1123       case PM_SOURCE_PHOTFIT_CONST:
    1124         break;
    1125       case PM_SOURCE_PHOTFIT_IMAGE_VAR:
    1126         Wi = Mi->variance;
    1127         psAssert (Wi, "programming error");
    1128         break;
    1129       case PM_SOURCE_PHOTFIT_MODEL_VAR:
    1130         Wi = Mi->modelVar;
    1131         psAssert (Wi, "programming error");
    1132         break;
    1133       case PM_SOURCE_PHOTFIT_NONE:
    1134         psAbort("programming error");
    1135     }   
    1136 # else
    1137     const psImage *Wi = Mi->variance;
    1138     if (!unweighted_sum) {
    1139         assert (Wi != NULL);
    1140     }
    1141 # endif
    1142 
    1143     const psImage *Ti = Mi->maskObj;
    1144     assert (Ti != NULL);
    1145     const psImage *Tj = Mj->maskObj;
    1146     assert (Tj != NULL);
    1147 
    1148     Xs = PS_MAX (Pi->col0, Pj->col0);
    1149     Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);
    1150 
    1151     Ys = PS_MAX (Pi->row0, Pj->row0);
    1152     Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);
    1153 
    1154     xIs = Xs - Pi->col0;
    1155     xJs = Xs - Pj->col0;
    1156     yIs = Ys - Pi->row0;
    1157     yJs = Ys - Pj->row0;
    1158 
    1159     xIe = Xe - Pi->col0;
    1160     yIe = Ye - Pi->row0;
    1161 
    1162     // note that weight is addressing the same image pixels,
    1163     flux = 0;
    1164     for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
    1165         for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
    1166             if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
    1167                 continue;
    1168             if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj] & maskVal)
    1169                 continue;
    1170 
    1171 # if (HAVE_MODEL_VAR)
    1172             float value = (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
    1173             switch (fitVarMode) {
    1174               case PM_SOURCE_PHOTFIT_CONST:
    1175                 wt = 1.0;
    1176                 break;
    1177               case PM_SOURCE_PHOTFIT_IMAGE_VAR:
    1178               case PM_SOURCE_PHOTFIT_MODEL_VAR:
    1179                 wt = covarFactor * Wi->data.F32[yi][xi];
    1180                 break;
    1181               case PM_SOURCE_PHOTFIT_NONE:
    1182                 psAbort("programming error");
    1183             }
    1184             // skip pixels with nonsense weight values
    1185             if (wt <= 0) continue;
    1186 
    1187             flux += value / wt;
    1188 
    1189 # else
    1190             // XXX skip the nonsense weight pixels?
    1191             if (unweighted_sum) {
    1192                 flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
    1193             } else {
    1194                 wt = covarFactor * Wi->data.F32[yi][xi];
    1195                 if (wt > 0) {
    1196                     flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
    1197                 }
    1198             }
    1199 # endif
    1200         }
    1201     }
    1202     return flux;
    1203 }

  • branches/eam_branches/ipp-20120601/psModules/src/objects/pmSourcePhotometry.h

    r34049 r34053  
    3838} pmSourcePhotometryMode;
    3939
    40 # if (HAVE_MODEL_VAR)
    4140typedef enum {
    4241    PM_SOURCE_PHOTFIT_NONE       = 0,
    … …  
    4443    PM_SOURCE_PHOTFIT_IMAGE_VAR  = 2,
    4544    PM_SOURCE_PHOTFIT_MODEL_VAR  = 3,
     45    PM_SOURCE_PHOTFIT_MODEL_SKY  = 4,   // XXX bad name: set variance floor based on mean variance image (variance of sky)
    4646} pmSourceFitVarMode;
    47 # endif
    4847
    4948bool pmSourcePhotometryModel(
    … …  
    8483bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
    8584
    86 # if (HAVE_MODEL_VAR)
    8785double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal);
    8886double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal);
    8987double pmSourceModelWeight(const pmSource *Mi, int term, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal);
    90 # else
    91 double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);
    92 double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);
    93 double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);
    94 # endif
    9588
    9689bool pmSourceNeighborFlags (pmSource *source);

  • branches/eam_branches/ipp-20120601/psphot/src/psphot.h

    r34049 r34053  
    6161bool            psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
    6262
     63bool            psphotUpdateVariance (pmConfig *config, const pmFPAview *view, const char *filerule);
     64bool            psphotUpdateVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
     65
    6366bool            psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filerule);
    6467bool            psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
    … …  
    99102
    100103bool            psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final);
    101 
    102 # if (HAVE_MODEL_VAR)
    103104bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode);
    104 # else
    105 bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final);
    106 # endif
    107105
    108106bool            psphotSourceSize (pmConfig *config, const pmFPAview *view, const char *filerule, bool getPSFsize);

  • branches/eam_branches/ipp-20120601/psphot/src/psphotEfficiency.c

    r34049 r34053  
    420420
    421421    // psphotFitSourcesLinearReadout subtracts the model fits
    422 # if (HAVE_MODEL_VAR)
    423422    if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true, PM_SOURCE_PHOTFIT_CONST)) {
    424 # else
    425     if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true)) {
    426 # endif
    427423        psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources.");
    428424        psFree(fakeSources);

  • branches/eam_branches/ipp-20120601/psphot/src/psphotFitSourcesLinear.c

    r34049 r34053  
    1212static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal);
    1313
    14 # if (HAVE_MODEL_VAR)
    1514bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources);
    1615pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe);
    1716bool psphotFreeModelVariance (pmReadout *readout, psArray *sources);
    18 # endif
    1917
    2018// for now, let's store the detections on the readout->analysis for each readout
    … …  
    3028    assert (recipe);
    3129
    32 # if (HAVE_MODEL_VAR)
    3330    pmSourceFitVarMode fitVarMode = psphotGetFitVarMode (recipe);
    3431    if (!fitVarMode) {
    … …  
    4037    // do a single pass.
    4138    pmSourceFitVarMode fitVarModePass1 = (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) ? PM_SOURCE_PHOTFIT_CONST : fitVarMode;
    42 # endif
    4339
    4440    int num = psphotFileruleCount(config, filerule);
    … …  
    6864        psAssert (psf, "missing psf?");
    6965
    70 # if (HAVE_MODEL_VAR)
    71         if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1))
    72 # else
    73         if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final))
    74 # endif
    75         {
     66        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1)) {
    7667            psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
    7768            return false;
    7869        }
    7970
    80 # if (HAVE_MODEL_VAR)
    8171        // the MODEL_VAR weighting scheme requires knowledge of the model fluxes to generate the variance
    8272        // after we have determined the initial set of fits, then we can generate the variance image and
    … …  
    10494            }
    10595        }
    106 # endif
    10796
    10897        psphotVisualShowResidualImage (readout, (num > 0));
    … …  
    113102}
    114103
    115 # if (HAVE_MODEL_VAR)
    116104// look up the fit variance mode from the recipe; older recipes do not have the value
    117105// 'LINEAR_FIT_VARIANCE_MODE'; in those cases, look for 'CONSTANT_PHOTOMETRIC_WEIGHTS' as a boolean and
    … …  
    136124        return PM_SOURCE_PHOTFIT_IMAGE_VAR;
    137125    }
     126    if (!strcasecmp(fitVarModeString, "SKY")   || !strcasecmp(fitVarModeString, "MODEL_SKY")) {
     127        return PM_SOURCE_PHOTFIT_MODEL_SKY;
     128    }
    138129    if (!strcasecmp(fitVarModeString, "MODEL") || !strcasecmp(fitVarModeString, "MODEL_VAR")) {
    139130        return PM_SOURCE_PHOTFIT_MODEL_VAR;
    … …  
    142133    return PM_SOURCE_PHOTFIT_NONE;
    143134}
    144 # endif
    145 
    146 # if (HAVE_MODEL_VAR)
    147 bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode)
    148 # else
    149 bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final)
    150 # endif
    151 {
     135
     136bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode) {
    152137    bool status;
    153138    float x;
    … …  
    185170    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
    186171
    187 # if (!HAVE_MODEL_VAR)
    188     bool CONSTANT_PHOTOMETRIC_WEIGHTS =
    189         psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
    190     if (!status) {
    191         psAbort("You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
    192     }
    193 # endif
    194172    int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER");
    195173    if (!status) {
    … …  
    319297    psSparseBorder *border = psSparseBorderAlloc (sparse, nBorder);
    320298
    321 # if (HAVE_MODEL_VAR)
    322299    // if fitVarMode is MODEL_VAR, then we need to generate the model image variance
    323300    // XXX we have two possibilities here:
    … …  
    327304
    328305    // 2) do a single pass, and use the model guess to define the model variance (but do I trust the Model Guess?)
    329 # endif
    330306
    331307    // fill out the sparse matrix elements and border elements (B)
    … …  
    336312
    337313        // diagonal elements of the sparse matrix (auto-cross-product)
    338 # if (HAVE_MODEL_VAR)
    339314        f = pmSourceModelDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
    340 # else
    341         f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    342 # endif
    343315        psSparseMatrixElement (sparse, i, i, f);
    344316
    345 # if (HAVE_MODEL_VAR)
    346         // if we have used CONSTANT errors, then we need to calculate the value of the parameter error
    347         if (fitVarMode != PM_SOURCE_PHOTFIT_IMAGE_VAR) {
     317        // if we have used CONSTANT errors, then we need to re-calculate the value of the parameter error
     318        if (fitVarMode == PM_SOURCE_PHOTFIT_CONST) {
    348319            float var = pmSourceModelDotModel (SRCi, SRCi, PM_SOURCE_PHOTFIT_IMAGE_VAR, covarFactor, maskVal);
    349320            errors->data.F32[i] = 1.0 / sqrt(var);
    … …  
    351322            errors->data.F32[i] = 1.0 / sqrt(f);
    352323        }
    353 # else
    354         // the formal error depends on the weighting scheme
    355         if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
    356             float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor, maskVal);
    357             errors->data.F32[i] = 1.0 / sqrt(var);
    358         } else {
    359             errors->data.F32[i] = 1.0 / sqrt(f);
    360         }
    361 # endif
    362324
    363325        // find the image x model value
    364 # if (HAVE_MODEL_VAR)
    365326        f = pmSourceDataDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
    366 # else
    367         f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    368 # endif
    369327        psSparseVectorElement (sparse, i, f);
    370328
    … …  
    372330        switch (SKY_FIT_ORDER) {
    373331          case 1:
    374 # if (HAVE_MODEL_VAR)
    375332            f = pmSourceModelWeight (SRCi, 1, fitVarMode, covarFactor, maskVal);
    376333            psSparseBorderElementB (border, i, 1, f);
    377334            f = pmSourceModelWeight (SRCi, 2, fitVarMode, covarFactor, maskVal);
    378335            psSparseBorderElementB (border, i, 2, f);
    379 # else
    380             f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    381             psSparseBorderElementB (border, i, 1, f);
    382             f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    383             psSparseBorderElementB (border, i, 2, f);
    384 # endif
    385336
    386337          case 0:
    387 # if (HAVE_MODEL_VAR)
    388338            f = pmSourceModelWeight (SRCi, 0, fitVarMode, covarFactor, maskVal);
    389 # else
    390             f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    391 # endif
    392339            psSparseBorderElementB (border, i, 0, f);
    393340            break;
    … …  
    409356
    410357            // got an overlap; calculate cross-product and add to output array
    411 # if (HAVE_MODEL_VAR)
    412358            f = pmSourceModelDotModel (SRCi, SRCj, fitVarMode, covarFactor, maskVal);
    413 # else
    414             f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
    415 # endif
    416359            psSparseMatrixElement (sparse, j, i, f);
    417360        }
    … …  
    451394
    452395    // set the sky, sky_x, sky_y components of border matrix
    453 # if (HAVE_MODEL_VAR)
    454396    SetBorderMatrixElements (border, readout, fitSources, (fitVarMode == PM_SOURCE_PHOTFIT_CONST), SKY_FIT_ORDER, markVal);
    455 # else
    456     SetBorderMatrixElements (border, readout, fitSources, CONSTANT_PHOTOMETRIC_WEIGHTS, SKY_FIT_ORDER, markVal);
    457 # endif
    458397
    459398    psSparseConstraint constraint;
    … …  
    614553}
    615554
    616 # if (HAVE_MODEL_VAR)
    617555bool psphotModelBackgroundReadout(psImage *model,  // Model image
    618556                                  psImage *modelStdev, // Model stdev image
    … …  
    650588    psImage *varModelStdev = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
    651589
     590    // generate an image of the mean variance image in DN
    652591    if (!psphotModelBackgroundReadout(varModel, varModelStdev, NULL, readout, backBinning, config, true)) {
    653592        psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
    … …  
    668607    psFree (varModelStdev);
    669608
    670     // XXX for a test:
    671     psphotSaveImage (NULL, modelVar, "model.bck.fits");
     609    float gain = 1.0;  // accept 1.0 as a default since it is not critical to the analysis
     610    pmCell *cell = readout->parent; // The parent cell
     611    if (cell) {
     612      gain = psMetadataLookupF32(&status, cell->concepts, "CELL.GAIN"); // Cell gain
     613      if (!status) {
     614        gain = 1.0;           // set note above
     615      }
     616    }
     617    if (gain > 2.0) { /* warn? */ }
    672618
    673619    // insert all of the source models
    … …  
    690636
    691637        // add the source model to the model variance image
     638        // XXX note that this should be added with gain applied
     639        // var_DN = flux_DN / gain [e/DN]
     640        // to do this requires an API upgrade...
    692641        pmSourceAdd (source, PM_MODEL_OP_MODELVAR, maskVal);
    693642    }
    694 
    695     // XXX for a test:
    696     psphotSaveImage (NULL, modelVar, "model.var.fits");
    697     psphotSaveImage (NULL, readout->variance, "image.var.fits");
    698643
    699644    // we save the model variance for future reference
    … …  
    724669    return true;
    725670}
    726 # endif

  • branches/eam_branches/ipp-20120601/psphot/src/psphotMaskReadout.c

    r29936 r34053  
    9494
    9595    // test output of files at this stage
    96     if (psTraceGetLevel("psphot") >= 5) {
     96    if (psTraceGetLevel("psphot.imsave") >= 5) {
    9797        psphotSaveImage (NULL, readout->image,  "image.fits");
    9898        psphotSaveImage (NULL, readout->mask,   "mask.fits");
    … …  
    105105    return true;
    106106}
     107
     108// XXX this function and support below was created to test the theory that the faint-end
     109// bias results from the Poisson variation of the background pixels.  This is NOT the
     110// case.  Using the code below maintains the faint-end bias.
     111bool psphotUpdateVariance (pmConfig *config, const pmFPAview *view, const char *filerule) {
     112
     113    bool status = false;
     114
     115    // select the appropriate recipe information
     116    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     117    psAssert (recipe, "missing recipe?");
     118
     119    int num = psphotFileruleCount(config, filerule);
     120
     121    // loop over the available readouts
     122    for (int i = 0; i < num; i++) {
     123
     124        // Generate the mask and weight images, including the user-defined analysis region of interest
     125        if (!psphotUpdateVarianceReadout (config, view, filerule, i, recipe)) {
     126            psError (PSPHOT_ERR_CONFIG, false, "failed to generate mask for %s entry %d", filerule, i);
     127            return false;
     128        }
     129    }
     130    return true;
     131}
     132
     133// determine the mean variance image (equivalent to the background model, but for the variance image)
     134// set the variance image to the MAX(input, mean)
     135bool psphotUpdateVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
     136
     137    bool status;
     138
     139    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
     140    psAssert (file, "missing file?");
     141
     142    // find the currently selected readout
     143    pmReadout  *readout = pmFPAviewThisReadout (view, file->fpa);
     144    psAssert (readout, "missing readout?");
     145
     146    pmSourceFitVarMode varMode = psphotGetFitVarMode (recipe);
     147    if (varMode == PM_SOURCE_PHOTFIT_NONE) {
     148      psError (PSPHOT_ERR_CONFIG, false, "need valid LINEAR_FIT_VARIANCE_MODE");
     149      return false;
     150    }
     151
     152    // make this an option via the recipe
     153    if (varMode != PM_SOURCE_PHOTFIT_MODEL_SKY) return true;
     154
     155    // create a model variance image (full-scale image to take result of psImageUnbin below)
     156    psImage *modelVar = psImageCopy (NULL, readout->variance, PS_TYPE_F32);
     157
     158    // find the binning information
     159    psImageBinning *backBinning = psphotBackgroundBinning (modelVar, config);
     160    assert (backBinning);
     161   
     162    psImage *varModel = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
     163    psImage *varModelStdev = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
     164
     165    if (!psphotModelBackgroundReadout(varModel, varModelStdev, NULL, readout, backBinning, config, true)) {
     166        psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
     167        psFree (varModel);
     168        psFree (varModelStdev);
     169        return false;
     170    }
     171
     172    // linear interpolation to full-scale
     173    if (!psImageUnbin (modelVar, varModel, backBinning)) {
     174        psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning");
     175        psFree (varModel);
     176        psFree (varModelStdev);
     177        return false;
     178    }
     179
     180    // XXX save these?
     181    psFree (varModel);
     182    psFree (varModelStdev);
     183
     184    psImage *im = readout->image;
     185    psImage *wt = readout->variance;
     186    for (int j = 0; j < im->numRows; j++) {
     187      for (int i = 0; i < im->numCols; i++) {
     188        if (!isfinite(im->data.F32[j][i])) continue;
     189        if (!isfinite(wt->data.F32[j][i])) continue;
     190        // XXX for a test, make variance constant wt->data.F32[j][i] = PS_MAX(wt->data.F32[j][i], modelVar->data.F32[j][i]);
     191        wt->data.F32[j][i] = modelVar->data.F32[j][i];
     192      }
     193    }
     194
     195    // test output of files at this stage
     196    if (psTraceGetLevel("psphot.imsave") >= 5) {
     197        psphotSaveImage (NULL, readout->image,  "image.varsky.fits");
     198        psphotSaveImage (NULL, readout->mask,   "mask.varsky.fits");
     199        psphotSaveImage (NULL, readout->variance, "variance.varsky.fits");
     200    }
     201
     202    psFree (modelVar);
     203
     204    return true;
     205}

  • branches/eam_branches/ipp-20120601/psphot/test/tap_psphot_varmodel.pro

    r33963 r34053  
    22# -*-sh-*-
    33
    4 # $KAPA = kapa -noX
    5 
    6 # PSF.CONVOLVE : if true, we insert delta functions (and optionally
    7 #                galaxies) and smooth the image with the psf model
    8 #                (uses a GAUSS regardless of the model). Note that
    9 #                PSF.CONVOLVE = T is faster than F, but (a) only
    10 #                allows Gauss models and (b) only yields quantized
    11 #                locations
    12 
    13 # config for ppImage to generate chip, mask, weight
    14 $ppImageConfig = -recipe PPIMAGE PPIMAGE_N
    15 $ppImageConfig = $ppImageConfig -Db BACKGROUND T
    16 $ppImageConfig = $ppImageConfig -Db CHIP.FITS T
    17 $ppImageConfig = $ppImageConfig -Db CHIP.MASK.FITS T
    18 $ppImageConfig = $ppImageConfig -Db CHIP.VARIANCE.FITS T
    19 $ppImageConfig = $ppImageConfig -Db BASE.FITS F
    20 $ppImageConfig = $ppImageConfig -Db VARIANCE.BUILD T
    21 $ppImageConfig = $ppImageConfig -Db PHOTOM F
    22 
    23 # basic options for the these images (filter, location, obstype)
    24 $BaseOptions = -type OBJECT -filter r -skymags 20.86 -ra 270.70 -dec -23.70 -pa 0.0
    25 $BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0
    26 
    27 # create an image with fake sources and insert the resulting cmf file into a dvodb
    28 $RefConfig = -camera SIMTEST
    29 $RefConfig = $RefConfig -recipe PPSIM STACKTEST.MAKE
    30 $RefConfig = $RefConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
    31 $RefConfig = $RefConfig -Db PSF.CONVOLVE F
    32 
    33 # options for the reference image
    34 $RefOptions = $BaseOptions -exptime 100.0
    35 $RefOptions = $RefOptions -seeing 1.0
    36 $RefOptions = $RefOptions -D PSF.MODEL PS_MODEL_GAUSS
    37 $RefOptions = $RefOptions -Df STARS.DENSITY 10.0
    38 $RefOptions = $RefOptions -Df STARS.SIGMA.LIM 0.5
    39 
    40 # basic config for ppSim with randomly distributed stars and NO galaxies
    41 $RealConfig = -camera SIMTEST
    42 $RealConfig = $RealConfig -recipe PPSIM STACKTEST.RUN
    43 $RealConfig = $RealConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
    44 $RealConfig = $RealConfig -Db STARS.FAKE F
    45 $RealConfig = $RealConfig -Db STARS.REAL T
    46 $RealConfig = $RealConfig -Db MATCH.DENSITY F
    47 $RealConfig = $RealConfig -Db PSF.CONVOLVE F
    48 $RealConfig = $RealConfig -Df STARS.DENSITY 10.0
    49 $RealConfig = $RealConfig -Df STARS.SIGMA.LIM 2.5
    50 $RealConfig = $RealConfig -Db GALAXY.FAKE F
    51 $RealConfig = $RealConfig -Db GALAXY.GRID F
    52 
    53 # options for the repeated images
    54 $RealOptions = $BaseOptions -exptime 30.0
    55  
    56 $ExtraOptions = -D PSF.MODEL PS_MODEL_GAUSS
    57 
    58 # sample alternate options:
    59 # $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_PS1_V1
    60 # $ppSimOptions = $FakeOptions -Df PSF.ARATIO 1.2
    61 # $ppSimOptions = $FakeOptions -Df PSF.THETA +30.0
    62 # $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS
    63 
    64 list fwhm
    65  1.0
    66  1.1
    67  1.2
    68  1.5
    69 end
     4# This script includes a set of tests to demonstrate the dependence of the faint-end bias on the weighting scheme
     5# We have 3 weighting schemes :
     6# CONSTANT -- per-pixel weight is fixed (disadvantage: lower S/N, especially for higher sky?)
     7# IMAGE_VAR - per-pixel weight is Poisson from image (disadvantage: faint-end bias)
     8# MODEL_VAR - per-pixel weight is Poisson from model (disadvantage: 2 linear-fit passes)
     9
     10# Functions:
     11#
     12# init  : initialize variables
     13# mkref : generate a fake reference catalog in a DVO database
     14# mkexp : generate a fake exposure from fake catalog and detrend it (saves basename.in.cmf & basename.fits)
     15# runphot : run psphot on an exposure (saves basename.cmf)
     16# ckchip.mags : generate a set of summary plots for the psphot analysis
     17
     18# I would like to produce a grid of tests:
     19# sky (bright, middle, dark)
     20# fwhm (1.0, 1.3, 1.6)
     21# PSF_MODEL input (GAUSS PS1_V1)
     22# PSF_MODEL apply (GAUSS PS1_V1)
     23# variance mode: CONSTANT, IMAGE_VAR, MODEL_VAR
    7024
    7125macro go
    72 
    7326  mkdir test
    7427
    75   $ExtraOptions = -D PSF.MODEL PS_MODEL_GAUSS
    76   mkexp test/image.00 1.0
    77   $ExtraOptions = -D PSF.MODEL PS_MODEL_PGAUSS
    78   mkexp test/image.01 1.0
    79   $ExtraOptions = -D PSF.MODEL PS_MODEL_PS1_V1
    80   mkexp test/image.02 1.0
     28  local sky fwhm psf_in psf_out
     29
     30  # mkref creates refimage.* and catdir.ref
     31  file catdir.ref found
     32  if (not($found))
     33    mkref
     34  end
     35
     36  $KAPA = kapa -noX
     37  if (not($?RefConfig)) init
     38
     39  foreach sky 19.0 20.0 21.0
     40    foreach fwhm 1.0 1.3 1.6
     41      foreach psf_in GAUSS PS1_V1
     42
     43        sprintf name "test/test.%02d.%02d.%s" $sky {10*$fwhm} $psf_in
     44        # mkexp $name $sky $fwhm $psf_in
     45       
     46        foreach psf_out GAUSS PS1_V1
     47          foreach mode CONSTANT IMAGE_VAR MODEL_VAR
     48            # runphot $name $name.$psf_out.$mode "-D LINEAR_FIT_VARIANCE_MODE $mode -D PSF_MODEL PS_MODEL_$psf_out"
     49            ckchip.mags $name.in.cmf $name.$psf_out.$mode.cmf $name.$psf_out.$mode 0.0
     50          end
     51        end
     52      end
     53    end
     54  end
    8155end
    8256
    8357# create a reference database of fake stars to be used by ppSim below
    8458macro mkref
     59  if (not($?RefConfig)) init
     60
    8561  exec rm -rf catdir.ref
    8662  exec rm -f refimage.fits
    … …  
    10076# create a realistic distribution of fake stars, GAUSS PSF
    10177macro mkexp
     78  if ($0 != 5)
     79    echo "USAGE: mkexp basename sky fwhm psf_model"
     80    break
     81  end
     82
     83  local fwhm basename psf_model
     84  $basename = $1
     85  $sky = $2
     86  $fwhm = $3
     87  $psf_model = $4
     88
     89  $ExtraOptions = -D PSF.MODEL PS_MODEL_$psf_model
     90
     91  # create the raw image
     92  echo ppSim -seeing $fwhm -skymags $sky -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename
     93  exec ppSim -seeing $fwhm -skymags $sky -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename
     94  exec /bin/mv -f $basename.cmf $basename.in.cmf
     95
     96  # create the chip output
     97  echo ppImage $ppImageConfig -file $basename.fits $basename
     98  exec ppImage $ppImageConfig -file $basename.fits $basename
     99end
     100
     101# create a realistic distribution of fake stars, GAUSS PSF
     102macro mkexp.deep
    102103  if ($0 != 3)
    103104    echo "USAGE: mkexp basename fwhm"
    … …  
    109110  $fwhm = $2
    110111
     112  $RealOptionsDeep = $BaseOptions -exptime 100
     113
    111114  # create the raw image
    112   echo ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename
    113   exec ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename
     115  echo ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptionsDeep $ExtraOptions $basename
     116  exec ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptionsDeep $ExtraOptions $basename
     117  exec /bin/mv -f $basename.cmf $basename.in.cmf
     118
     119  # create the chip output
     120  echo ppImage $ppImageConfig -file $basename.fits $basename
     121  exec ppImage $ppImageConfig -file $basename.fits $basename
     122end
     123
     124# create a realistic distribution of fake stars, GAUSS PSF
     125macro mkexp.bright
     126  if ($0 != 3)
     127    echo "USAGE: mkexp basename fwhm"
     128    break
     129  end
     130
     131  local fwhm basename
     132  $basename = $1
     133  $fwhm = $2
     134
     135  # basic options for the these images (filter, location, obstype)
     136  $BaseOptions = -type OBJECT -filter r -skymags 20.0 -ra 270.70 -dec -23.70 -pa 0.0
     137  $BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0
     138  $RealOptionsDeep = $BaseOptions -exptime 100
     139
     140  # create the raw image
     141  echo ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptionsDeep $ExtraOptions $basename
     142  exec ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptionsDeep $ExtraOptions $basename
    114143  exec /bin/mv -f $basename.cmf $basename.in.cmf
    115144
    … …  
    133162  echo psphot -threads 4 -file $basename.ch.fits -mask $basename.ch.mk.fits -variance $basename.ch.wt.fits $outname $options
    134163  exec psphot -threads 4 -file $basename.ch.fits -mask $basename.ch.mk.fits -variance $basename.ch.wt.fits $outname $options
     164end
     165
     166# compare two cmf files with extname Chip.psf
     167# things to compare:
     168# * completeness (which sources in (1) are not detected in (2)
     169# * positions (X_PSF, Y_PSF)
     170# * instrumental psf mags
     171# * position errors (no input errors; use a model?)
     172# * measured FWHM?
     173# * kron mags (fluxes)
     174# * etc, etc
     175macro ckchip.mags
     176  if ($0 != 5)
     177    echo "USAGE: ckchip.mags (raw) (out) (output) (zpt_off)"
     178    break
     179  end
     180
     181list pairs
     182  PSF_INST_MAG_out      M_raw                 PSF_INST_MAG_raw 2 -0.21 0.21 V
     183  AP_MAG_out            M_raw                 PSF_INST_MAG_raw 2 -0.21 0.21 V
     184end
     185
     186  load.cmf $1 Chip.psf raw
     187  load.cmf $2 Chip.psf out
     188
     189  # images generated with convolution will not have the right output positions
     190  set X_raw = int(X_PSF_raw) + 0.5
     191  set Y_raw = int(Y_PSF_raw) + 0.5
     192  set M_raw = PSF_INST_MAG_raw + $4
     193  set K_out = -2.5*log(KRON_FLUX_out)
     194  match2d X_PSF_out Y_PSF_out X_PSF_raw Y_PSF_raw 1.5 -index1 index1 -index2 index2
     195
     196  local i NX NY nx ny N
     197
     198  resize 1000 1000
     199
     200  clear
     201  section a0 0.0 0.0 1.0 0.5
     202  label -fn courier 14;
     203  style -pt 0 -sz 0.4
     204  show.pair 0
     205  delete -q imag_V dmag_V smag_V Dmag_V dmag_Vraw
     206  for imag -11 -6 0.2
     207    subset dmsub = delta if (rv > $imag) && (rv < $imag + 0.2)
     208    vstat -q dmsub
     209    concat $imag imag_V
     210    concat $MEAN dmag_V
     211    concat $MEDIAN Dmag_V
     212    concat $SIGMA smag_V
     213  end
     214
     215  vstat -q dmag_V
     216  $offset = $MEDIAN
     217  set dmag_V = dmag_V - $offset
     218  set Dmag_V = Dmag_V - $offset
     219
     220  section a1 0.0 0.50 1.0 0.25
     221  lim rv -0.025 0.075; label -fn courier 14; box;
     222  plot -c red   -pt 7 -sz 2.0 imag_V dmag_V
     223  plot -c darkgreen -pt 3 -sz 2.0 imag_V Dmag_V
     224  # plot -c gold -pt 4 -sz 2.0 imag_V dmag_Vraw
     225  plot -c blue  -pt 2 -sz 2.0 imag_V smag_V
     226  # label -y "mean (red), median (green), unclipped mean (gold), sigma (blue) of delta"
     227  label -y "mean (red), median (green), sigma (blue) of delta"
     228  sprintf line "OFFSET: %6.3f" $offset
     229  textline -frac 0.1 0.8 -fn courier 24 "$line"
     230
     231  set lChiNorm = log(PSF_CHISQ_out / PSF_NDOF_out)
     232  reindex chi = lChiNorm using index1
     233  section a2 0.0 0.75 1.0 0.25
     234  label -fn courier 14;
     235  lim rv chi; box;
     236  plot -c red -pt 0 -sz 0.5 rv chi
     237
     238  label -y "chisq"
     239  png -name $3.png
     240
     241  # section a1 0.0 0.5 1.0 0.5
     242  # style -pt 0 -sz 0.4
     243  # show.pair 1
     244end
     245
     246macro go.phot
     247 runphot test.00 test.00.varmode.C "-Db SAVE.RESID T -D LINEAR_FIT_VARIANCE_MODE CONSTANT"
     248 runphot test.00 test.00.varmode.I "-Db SAVE.RESID T -D LINEAR_FIT_VARIANCE_MODE IMAGE_VAR"
     249 runphot test.00 test.00.varmode.M "-Db SAVE.RESID T -D LINEAR_FIT_VARIANCE_MODE MODEL_VAR"
     250 runphot test.00 test.00.varmode.S "-Db SAVE.RESID T -D LINEAR_FIT_VARIANCE_MODE MODEL_SKY"
     251end
     252
     253macro go.vars
     254 dev -n varC; ckchip.mags test.00.in.cmf test.00.varmode.C.cmf test.00.varmode.C 0.0
     255 dev -n varI; ckchip.mags test.00.in.cmf test.00.varmode.I.cmf test.00.varmode.I 0.0
     256 dev -n varM; ckchip.mags test.00.in.cmf test.00.varmode.M.cmf test.00.varmode.M 0.0
     257 dev -n varS; ckchip.mags test.00.in.cmf test.00.varmode.S.cmf test.00.varmode.S 0.0
    135258end
    136259
    … …  
    347470    lim rv $word:4 $word:5; box; plot rv delta
    348471  end
    349   label -y '$word:0' -x '$word:2'
     472  $line = '$word:0' - '$word:1'
     473  label -y "$line" -x '$word:2'
    350474end
    351475
    … …  
    637761end
    638762
     763macro init
     764  # config for ppImage to generate chip, mask, weight
     765  $ppImageConfig = -recipe PPIMAGE PPIMAGE_N
     766  $ppImageConfig = $ppImageConfig -Db BACKGROUND T
     767  $ppImageConfig = $ppImageConfig -Db CHIP.FITS T
     768  $ppImageConfig = $ppImageConfig -Db CHIP.MASK.FITS T
     769  $ppImageConfig = $ppImageConfig -Db CHIP.VARIANCE.FITS T
     770  $ppImageConfig = $ppImageConfig -Db BASE.FITS F
     771  $ppImageConfig = $ppImageConfig -Db VARIANCE.BUILD T
     772  $ppImageConfig = $ppImageConfig -Db PHOTOM F
     773 
     774  # basic options for the these images (filter, location, obstype)
     775  $BaseOptions = -type OBJECT -filter r -ra 270.70 -dec -23.70 -pa 0.0
     776  $BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0
     777 
     778  # PSF.CONVOLVE : if true, we insert delta functions (and optionally
     779  #                galaxies) and smooth the image with the psf model
     780  #                (uses a GAUSS regardless of the model). Note that
     781  #                PSF.CONVOLVE = T is faster than F, but (a) only
     782  #                allows Gauss models and (b) only yields quantized
     783  #                locations
     784
     785  # create an image with fake sources and insert the resulting cmf file into a dvodb
     786  $RefConfig = -camera SIMTEST
     787  $RefConfig = $RefConfig -recipe PPSIM STACKTEST.MAKE
     788  $RefConfig = $RefConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
     789  $RefConfig = $RefConfig -Db PSF.CONVOLVE F
     790 
     791  # options for the reference image
     792  $RefOptions = $BaseOptions
     793  $RefOptions = $RefOptions -exptime 100.0
     794  $RefOptions = $RefOptions -seeing 1.0
     795  $RefOptions = $RefOptions -skymags 21.0 
     796  $RefOptions = $RefOptions -D PSF.MODEL PS_MODEL_GAUSS
     797  $RefOptions = $RefOptions -Df STARS.DENSITY 10.0
     798  $RefOptions = $RefOptions -Df STARS.SIGMA.LIM 0.5
     799
     800  # basic config for ppSim with randomly distributed stars and NO galaxies
     801  $RealConfig = -camera SIMTEST
     802  $RealConfig = $RealConfig -recipe PPSIM STACKTEST.RUN
     803  $RealConfig = $RealConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
     804  $RealConfig = $RealConfig -Db STARS.FAKE F
     805  $RealConfig = $RealConfig -Db STARS.REAL T
     806  $RealConfig = $RealConfig -Db MATCH.DENSITY F
     807  $RealConfig = $RealConfig -Db PSF.CONVOLVE F
     808  $RealConfig = $RealConfig -Df STARS.DENSITY 10.0
     809  $RealConfig = $RealConfig -Df STARS.SIGMA.LIM 2.5
     810  $RealConfig = $RealConfig -Db GALAXY.FAKE F
     811  $RealConfig = $RealConfig -Db GALAXY.GRID F
     812 
     813  # options for the repeated images
     814  $RealOptions = $BaseOptions -exptime 30.0
     815   
     816  # sample alternate options:
     817  # $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_PS1_V1
     818  # $ppSimOptions = $FakeOptions -Df PSF.ARATIO 1.2
     819  # $ppSimOptions = $FakeOptions -Df PSF.THETA +30.0
     820  # $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS
     821 
     822  list fwhm
     823   1.0
     824   1.1
     825   1.2
     826   1.5
     827  end
     828end
     829
    639830if ($SCRIPT)
    640831  fulltest 4
Note: See TracChangeset for help on using the changeset viewer.