IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 30, 2012, 2:49:37 PM (14 years ago)
Author:
eugene
Message:

merge changes from trunk

Location:
branches/eam_branches/ipp-20111122/psModules/src
Files:
25 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20111122/psModules/src/camera/pmReadoutFake.c

    r32347 r33638  
    210210    const pmPSF *psf = args->data[7];         // PSF
    211211    float minFlux = PS_SCALAR_VALUE(args->data[8], F32); // Minimum flux
    212     float radius = PS_SCALAR_VALUE(args->data[9], F32);  // Minimum radius
     212    float radius = PS_SCALAR_VALUE(args->data[9], S32);  // Minimum radius - typecast to float from S32 outside of PS_SCALAR_VALUE otherwise sets 0.0
    213213    bool circularise = PS_SCALAR_VALUE(args->data[10], U8); // Circularise PSF?
    214214    bool normalisePeak = PS_SCALAR_VALUE(args->data[11], U8); // Normalise for peak?
    … …  
    314314                }
    315315            }
    316             if (!psThreadPoolWait(true)) {
     316            if (!psThreadPoolWait(true, true)) {
    317317                psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
    318318                psFree(groups);

  • branches/eam_branches/ipp-20111122/psModules/src/concepts/pmConceptsStandard.c

    r30049 r33638  
    751751  bool has_video_cell = false;
    752752
    753   if (concept->type != PS_DATA_STRING) {
    754     psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Type for %s (%x) is not string\n",
    755             concept->name, concept->type);
    756     return NULL;
    757   }
    758 
    759   char *Vptr = strchr(concept->data.V,'V');
    760   if (Vptr) {
    761     has_video_cell = true;
     753  if (concept->type == PS_DATA_BOOL) {
     754    has_video_cell = concept->data.B;
     755  } else {
     756    if (concept->type != PS_DATA_STRING) {
     757        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Type for %s (%x) is not string\n",
     758                concept->name, concept->type);
     759        return NULL;
     760      }
     761
     762      char *Vptr = strchr(concept->data.V,'V');
     763      if (Vptr) {
     764        has_video_cell = true;
     765      }
    762766  }
    763767

  • branches/eam_branches/ipp-20111122/psModules/src/detrend/pmBias.c

    r29833 r33638  
    154154    if (threaded) {
    155155        // wait here for the threaded jobs to finish
    156         if (!psThreadPoolWait(true)) {
     156        if (!psThreadPoolWait(true, true)) {
    157157            psError(PS_ERR_UNKNOWN, false, "Unable to apply bias correction.");
    158158            return false;

  • branches/eam_branches/ipp-20111122/psModules/src/detrend/pmDark.c

    r28405 r33638  
    601601    if (threaded) {
    602602        // wait here for the threaded jobs to finish
    603         if (!psThreadPoolWait(true)) {
     603        if (!psThreadPoolWait(true, true)) {
    604604            psError(PS_ERR_UNKNOWN, false, "Unable to apply dark.");
    605605            psFree(orders);

  • branches/eam_branches/ipp-20111122/psModules/src/detrend/pmFlatField.c

    r28405 r33638  
    161161    if (threaded) {
    162162        // wait here for the threaded jobs to finish
    163         if (!psThreadPoolWait(true)) {
     163        if (!psThreadPoolWait(true, true)) {
    164164            psError(PS_ERR_UNKNOWN, false, "Unable to flat-field image.");
    165165            return false;

  • branches/eam_branches/ipp-20111122/psModules/src/detrend/pmPattern.c

    r27676 r33638  
    66
    77#include "pmPattern.h"
     8
     9#define PATTERN_ROW_BKG_FIX 1
    810
    911
    … …  
    8991    psImageInit(corr, NAN);
    9092
     93#ifdef PATTERN_ROW_BKG_FIX
     94    // CZW: 2011-11-30
     95    // Define the vectors to hold the "x" and "y" slope trends.
     96    // Briefly, the slope trend in the y-axis is a due to variations in the 0-th order term
     97    // of the PATTERN.ROW fit between individual rows across the cell.  Similarly, the 1-st
     98    // order term of the PATTERN.ROW fit defines the trend in the x-axis (as that's what we
     99    // are fitting with PATTERN.ROW in the first place).  However, the thing we're trying to
     100    // fix with PATTERN.ROW is the detector level bias wiggles.  These should be overlaid on
     101    // the true sky level.  Therefore, simply applying the PATTERN.ROW correction will
     102    // introduce cell-to-cell sky variations as these two trends are removed.  To avoid this,
     103    // We store the 0th and 1st order values used for each row, and then fit a polynomial to
     104    // these results.  By re-adding these systematic trends back, we can remove the row-to-row
     105    // variations without improperly removing the real sky trend.
     106    psVector *yaxisData = psVectorAlloc(numRows, PS_TYPE_F32); // Data to fit to the constant term
     107    psVector *yaxisMask = psVectorAlloc(numRows, PS_TYPE_VECTOR_MASK); // Mask for rows with no fit
     108    psVector *xaxisData = psVectorAlloc(numRows, PS_TYPE_F32); // Data to fit to the linear term
     109    psVectorInit(yaxisMask, 0);
     110#endif
    91111    for (int y = 0; y < numRows; y++) {
    92112        psVectorInit(clipMask, 0);
    … …  
    105125            // Not enough points to fit
    106126            patternMaskRow(ro, y, maskBad);
     127#ifdef PATTERN_ROW_BKG_FIX
     128            // Ignore this row in our subsequent fits, because the fit failed.
     129            yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF;
     130#endif
    107131            continue;
    108132        }
    … …  
    111135            psErrorClear();
    112136            patternMaskRow(ro, y, maskBad);
    113             continue;
    114         }
    115 
    116         poly->coeff[0] -= background;
     137#ifdef PATTERN_ROW_BKG_FIX
     138            // Ignore this row in our subsequent fits, because the fit failed.
     139            yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF;
     140#endif
     141            continue;
     142        }
     143#ifndef PATTERN_ROW_BKG_FIX
     144        poly->coeff[0] -= background;
     145#else
     146        // Store the results we found for this row.
     147        yaxisData->data.F32[y] = poly->coeff[0];
     148        xaxisData->data.F32[y] = poly->coeff[1];
     149        psTrace("pattern",1,"%d %g %g\n",y,poly->coeff[0],poly->coeff[1]);
     150       
     151        //      yaxisData->data.F32[y] = 0.0;
     152/*      xaxisData->data.F32[y] = 0.0; */
     153       
     154#endif
    117155        memcpy(corr->data.F64[y], poly->coeff, (order + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
    118156        psVector *solution = psPolynomial1DEvalVector(poly, indices); // Solution vector
    … …  
    121159            psErrorClear();
    122160            patternMaskRow(ro, y, maskBad);
     161#ifdef PATTERN_ROW_BKG_FIX
     162            yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF;
     163#endif
    123164            continue;
    124165        }
    … …  
    126167        for (int x = 0; x < numCols; x++) {
    127168            image->data.F32[y][x] -= solution->data.F32[x];
     169            psTrace("pattern",5,"A: %d %d %g\n",x,y,solution->data.F32[x]);
    128170        }
    129171        psFree(solution);
    130172    }
    131173
     174#ifdef PATTERN_ROW_BKG_FIX
     175    // Put the global trends back that were removed by the PATTERN.ROW correction.
     176    // Set up the indices for the polynomial
     177    psVector *yaxisIndices = psVectorAlloc(numRows, PS_TYPE_F32);
     178    norm = 2.0 / (float)numRows;
     179    for (int y = 0; y < numRows; y++) {
     180      yaxisIndices->data.F32[y] = y * norm - 1.0;
     181      psTrace("psModules.detrend.pattern",10,"%d %f %f\n",y,yaxisIndices->data.F32[y],yaxisData->data.F32[y]);
     182    }
     183
     184    // Fit the trend of the constant term, producing the y-axis global trend
     185    psStatsInit(clip);
     186    psPolynomial1D *yaxisPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); // Polynomial to fit.
     187    if (!psVectorClipFitPolynomial1D(yaxisPoly,clip,yaxisMask,0xFF,yaxisData, NULL, yaxisIndices)) {
     188      psWarning("Unable to fit polynomial to y-axis trend");
     189      psErrorClear();
     190      // If we've failed, we need to do something, so add back in the background level, and
     191      // expect that the final image will have background mismatches.
     192      for (int y = 0; y < numRows; y++) {
     193        for (int x = 0; x < numCols; x++) {
     194          image->data.F32[y][x] += background;
     195          corr->data.F64[y][0]  -= background;
     196        }
     197      }
     198    }
     199    else {
     200      psVector *solution = psPolynomial1DEvalVector(yaxisPoly,yaxisIndices);
     201      if (!solution) {
     202        psWarning("Unable to evaluate polynomial");
     203        psErrorClear();
     204        // If we've failed, we need to do something, so add back in the background level, and
     205        // expect that the final image will have background mismatches.
     206        for (int y = 0; y < numRows; y++) {
     207          for (int x = 0; x < numCols; x++) {
     208            image->data.F32[y][x] += background;
     209            corr->data.F64[y][0]  -= background;
     210          }
     211        }
     212      }
     213      else {
     214        for (int y = 0; y < numRows; y++) {
     215          for (int x = 0; x < numCols; x++) {
     216            image->data.F32[y][x] += solution->data.F32[y];
     217            corr->data.F64[y][0]  -= solution->data.F32[y];
     218            psTrace("pattern",5,"B: %d %d %g\n",x,y,solution->data.F32[x]);
     219          }
     220        }
     221      }
     222      psFree(solution);
     223    }     
     224
     225    // Fit the trend of the linear term, producing the x-axis global trend
     226    // We can use the same mask vector, as the same rows failed the row-fit earlier.
     227    psStatsInit(clip);
     228    psPolynomial1D *xaxisPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); // Polynomial to fit.
     229    if (!psVectorClipFitPolynomial1D(xaxisPoly,clip,yaxisMask,0xFF,xaxisData, NULL, yaxisIndices)) {
     230      psWarning("Unable to fit polynomial to x-axis trend");
     231      psErrorClear();
     232    }
     233    else {
     234      psVector *solution = psPolynomial1DEvalVector(xaxisPoly,yaxisIndices);
     235      if (!solution) {
     236        psWarning("Unable to evaluate polynomial");
     237        psErrorClear();
     238      }
     239      else {
     240        for (int y = 0; y < numRows; y++) {
     241          for (int x = 0; x < numCols; x++) {
     242            image->data.F32[y][x] += solution->data.F32[y] * indices->data.F32[x];
     243            corr->data.F64[y][1]  -= solution->data.F32[y] ;
     244            psTrace("pattern",5,"C: %d %d %g %g\n",x,y,solution->data.F32[x],indices->data.F32[x]);
     245          }
     246        }
     247      }
     248      psFree(solution);
     249    }
     250    psFree(yaxisPoly);
     251    psFree(xaxisPoly);
     252    psFree(yaxisIndices);
     253    psFree(yaxisMask);
     254    psFree(yaxisData);
     255    psFree(xaxisData);
     256    // End PATTERN_ROW_BKG_FIX global trend replacement
     257#endif
     258   
    132259    psMetadataAddImage(ro->analysis, PS_LIST_TAIL, PM_PATTERN_ROW_CORRECTION, PS_META_REPLACE,
    133260                       "Pattern row correction", corr);
    … …  
    382509
    383510
     511
     512bool pmPatternContinuity(pmChip *chip, const psVector *tweak, psStatsOptions bgStat, psStatsOptions cellStat,
     513                         psImageMaskType maskVal, psImageMaskType maskBad, int edgeWidth)
     514{
     515    PS_ASSERT_PTR_NON_NULL(chip, false);
     516    PS_ASSERT_VECTOR_NON_NULL(tweak, false);
     517    PS_ASSERT_VECTOR_SIZE(tweak, chip->cells->n, false);
     518    PS_ASSERT_VECTOR_TYPE(tweak, PS_TYPE_U8, false);
     519
     520    int numCells = tweak->n;            // Number of cells
     521
     522    psVector *meanMask = psVectorAlloc(numCells, PS_TYPE_VECTOR_MASK); // Mask for means
     523    psVectorInit(meanMask, 0);
     524
     525    // Mask bits
     526    enum {
     527        PM_PATTERN_IGNORE = 0x01,       // Ignore this cell
     528        PM_PATTERN_TWEAK  = 0x02,       // Tweak this cell
     529        PM_PATTERN_ERROR  = 0x04,       // Error in calculating background
     530        PM_PATTERN_ALL    = 0xFF,       // All causes
     531    };
     532
     533    // Count number of cells to tweak
     534    int numTweak = 0;                   // Number of cells to tweak
     535    int numIgnore = 0;                  // Number of cells to ignore
     536    for (int i = 0; i < numCells; i++) {
     537        pmCell *cell = chip->cells->data[i]; // Cell of interest
     538        if (!cell || !cell->data_exists || !cell->process ||
     539            cell->readouts->n == 0 || cell->readouts->n > 1 || !cell->readouts->data[0]) {
     540            numIgnore++;
     541            meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PM_PATTERN_IGNORE;
     542            continue;
     543        }
     544        if (tweak->data.U8[i]) {
     545            meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PM_PATTERN_TWEAK;
     546            numTweak++;
     547        }
     548    }
     549    if (numTweak == 0) {
     550        // Nothing to do
     551        psFree(meanMask);
     552        return true;
     553    }
     554
     555    // Measure mean of each cell edge, and use that to determine the cell offsets.
     556
     557    psStats *bgStats = psStatsAlloc(bgStat); // Statistics on background
     558    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
     559
     560    psRegion region = {0,0,0,0};
     561
     562    /* These images hold the edge data for the OTA structure.  */
     563    psImage *A = psImageAlloc(8,8,PS_TYPE_F64); // Top edge
     564    psImage *B = psImageAlloc(8,8,PS_TYPE_F64); // Bottom edge
     565    psImage *C = psImageAlloc(8,8,PS_TYPE_F64); // Right edge
     566    psImage *D = psImageAlloc(8,8,PS_TYPE_F64); // Left edge
     567    psImageInit(A,0.0);
     568    psImageInit(B,0.0);
     569    psImageInit(C,0.0);
     570    psImageInit(D,0.0);
     571   
     572    for (int i = 0; i < numCells; i++) {
     573        if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_IGNORE) {
     574            continue;
     575        }
     576        pmCell *cell = chip->cells->data[i]; // Cell of interest
     577        pmReadout *ro = cell->readouts->data[0]; // Readout of interest
     578
     579        psStatsInit(bgStats);
     580
     581        // Convert cell iterator i into an xy coordinate on the grid of cells
     582        int y = (i % 8);
     583        int x = (i - y) / 8;
     584       
     585        for (int j = 0; j < 4; j++) {
     586          if (j == 0) {  // Region B
     587            region = psRegionSet(0,ro->image->numCols,
     588                                 0,edgeWidth);
     589          }
     590          else if (j == 1) { // Region A
     591            region = psRegionSet(0,ro->image->numCols,
     592                                 ro->image->numRows - edgeWidth,ro->image->numRows);
     593          }
     594          else if (j == 2) { // Region D
     595            region = psRegionSet(0,edgeWidth,
     596                                 0,ro->image->numRows);
     597          }
     598          else if (j == 3) { // Region C
     599            region = psRegionSet(ro->image->numCols - edgeWidth,ro->image->numCols,
     600                                 0,ro->image->numRows);
     601          }
     602          psImage *subset  = psImageSubset(ro->image,region);
     603          psImage *submask = psImageSubset(ro->mask,region);
     604
     605          if (!psImageBackground(bgStats, NULL, subset, submask, maskVal, rng)) {
     606            psWarning("Unable to measure background for cell %d on edge %d\n", i, j);
     607            psErrorClear();
     608            meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PM_PATTERN_ERROR;
     609            if (j == 0)      { B->data.F64[y][x] = NAN; }
     610            else if (j == 1) { A->data.F64[y][x] = NAN; }
     611            else if (j == 2) { C->data.F64[y][x] = NAN; }
     612            else if (j == 3) { D->data.F64[y][x] = NAN; }
     613            psFree(subset);
     614            psFree(submask);
     615            continue; // Move on to next edge, as only part of this cell may be a problem
     616          }
     617 
     618          // If the returned value is zero, assume something is wrong.  Do I still need this?
     619          if (psStatsGetValue(bgStats,bgStat) < 1e-6) {
     620            if (j == 0)      { B->data.F64[y][x] = NAN; }
     621            else if (j == 1) { A->data.F64[y][x] = NAN; }
     622            else if (j == 2) { C->data.F64[y][x] = NAN; }
     623            else if (j == 3) { D->data.F64[y][x] = NAN; }
     624          }
     625          // If we have an error for this cell/edge, make sure we mask the value
     626          if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_ERROR) {
     627            if (j == 0)      { B->data.F64[y][x] = NAN; }
     628            else if (j == 1) { A->data.F64[y][x] = NAN; }
     629            else if (j == 2) { C->data.F64[y][x] = NAN; }
     630            else if (j == 3) { D->data.F64[y][x] = NAN; }
     631          }
     632          else { // Set the value to match what we got from the edge box.
     633            if (j == 0)      { B->data.F64[y][x] = psStatsGetValue(bgStats,bgStat); }
     634            else if (j == 1) { A->data.F64[y][x] = psStatsGetValue(bgStats,bgStat); }
     635            else if (j == 2) { C->data.F64[y][x] = psStatsGetValue(bgStats,bgStat); }
     636            else if (j == 3) { D->data.F64[y][x] = psStatsGetValue(bgStats,bgStat); }
     637          }
     638
     639          for (int u = 0; u < subset->numCols; u++) {
     640            for (int v = 0; v < subset->numRows; v++) {
     641              psTrace("psModules.detrend.cont",10,"BOX: %d %d (%d %d) (%d %d) %f %d",
     642                      i,j,x,y,u,v,subset->data.F32[v][u],submask->data.PS_TYPE_IMAGE_MASK_DATA[v][u]);
     643            }
     644          }       
     645         
     646          psFree(subset);
     647          psFree(submask);
     648
     649        }
     650        psTrace("psModules.detrend.cont",5, "OTA: %d (%d %d) A: %f B: %f C: %f D: %f",
     651                i,x,y,
     652                A->data.F64[y][x],B->data.F64[y][x],C->data.F64[y][x],D->data.F64[y][x]);               
     653    }
     654    psFree(bgStats);
     655    psFree(rng);
     656
     657    // We've now allocated all the edge values, so we can now minimize the offsets.
     658    // This involves solving the equation A x = b, where
     659    // A is the (64x64 for GPC1) matrix containing the edges that match for each cell
     660    // x is the solution vector
     661    // b is the combination of offsets across each cell boundary for each cell.
     662    // Below "XX" is used as the matrix A, and "solution" is used as both b and x
     663    //   (due to the way psMatrixLUSolve operates).
     664    psVector *solution = psVectorAlloc(64,PS_TYPE_F64);
     665    psImage  *XX       = psImageAlloc(64,64,PS_TYPE_F64);
     666    psVectorInit(solution,0.0);
     667    psImageInit(XX,0.0);
     668   
     669    for (int i = 0; i < numCells; i++) {
     670      // Accumulate all the possible edge differences we can for this cell.
     671      // As we do so, make a note of the correlations by incrementing the element of the matrix.
     672      int y = (i % 8);
     673      int x = (i - y) / 8;
     674      int j;
     675      double critical_value = 0.0;
     676      if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_IGNORE) {
     677        continue;
     678      }
     679      if (x + 1 < 8) {  // We have a neighbor adjacent in the +x direction
     680        j = 8 * (x + 1) + y; // Determine that neighbor's index
     681        if (fabs(C->data.F64[y][x]) > fabs(D->data.F64[y][x+1])) {
     682          critical_value = 2.0 * fabs(D->data.F64[y][x+1]);
     683        }
     684        else {
     685          critical_value = 2.0 * fabs(C->data.F64[y][x]);
     686        }
     687        if (critical_value < 25) { critical_value = 25; }
     688        psTrace("psModules.detrend.cont",5,"CmD %d %d %d %d %g %g %g", // diagnostic
     689                i,x,y,j,
     690                C->data.F64[y][x],
     691                D->data.F64[y][x+1],
     692                critical_value
     693                );
     694        if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&&  // If there are no errors with the neighbor,
     695            (isfinite(C->data.F64[y][x]))&&(isfinite(D->data.F64[y][x+1]))&&     // and all edges have valid values,
     696            (fabs(C->data.F64[y][x] - D->data.F64[y][x+1]) < critical_value)     // and there are no large discontinuities,
     697            ) {   
     698          solution->data.F64[i] += C->data.F64[y][x] - D->data.F64[y][x+1];     // Take the difference
     699          XX->data.F64[i][i] += 1;                                              // increment our relation with ourself
     700          XX->data.F64[i][j] += -1;                                             // decrement our relation with the neighbor
     701        }
     702      }
     703      if (x - 1 > -1) { // etc.
     704        j = 8 * (x - 1) + y;
     705        if (fabs(C->data.F64[y][x-1]) > fabs(D->data.F64[y][x])) {
     706          critical_value = 2.0 * fabs(D->data.F64[y][x]);
     707        }
     708        else {
     709          critical_value = 2.0 * fabs(C->data.F64[y][x-1]);
     710        }
     711        if (critical_value < 25) { critical_value = 25; }
     712        psTrace("psModules.detrend.cont",5,"DmC %d %d %d %d %g %g %g",
     713                i,x,y,j,
     714                D->data.F64[y][x],
     715                C->data.F64[y][x-1],
     716                critical_value
     717                );
     718
     719        if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&&
     720            (isfinite(D->data.F64[y][x]))&&(isfinite(C->data.F64[y][x-1]))&&
     721            (fabs(D->data.F64[y][x] - C->data.F64[y][x-1]) < critical_value)
     722            ) {
     723          solution->data.F64[i] += D->data.F64[y][x] - C->data.F64[y][x-1];
     724          XX->data.F64[i][i] += 1;
     725          XX->data.F64[i][j] += -1;
     726        }
     727      }
     728      if (y + 1 < 8) {
     729        j = 8 * x + (y + 1);
     730        psTrace("psModules.detrend.cont",5,"AmB %d %d %d %d %g %g",
     731                i,x,y,j,
     732                A->data.F64[y][x],
     733                B->data.F64[y+1][x]
     734                );
     735        if (fabs(A->data.F64[y][x]) > fabs(B->data.F64[y+1][x])) {
     736          critical_value = 2.0 * fabs(B->data.F64[y+1][x]);
     737        }
     738        else {
     739          critical_value = 2.0 * fabs(A->data.F64[y][x]);
     740        }
     741        if (critical_value < 25) { critical_value = 25; }
     742        if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&&
     743            (isfinite(A->data.F64[y][x]))&&(isfinite(B->data.F64[y+1][x]))&&
     744            (fabs(A->data.F64[y][x] - B->data.F64[y+1][x]) < critical_value)
     745            ) {
     746          solution->data.F64[i] += A->data.F64[y][x] - B->data.F64[y+1][x];
     747          XX->data.F64[i][i] += 1;
     748          XX->data.F64[i][j] += -1;
     749        }
     750      }
     751      if (y - 1 > -1) {
     752        j = 8 * x +  (y - 1);
     753        psTrace("psModules.detrend.cont",5,"BmA %d %d %d %d %g %g",
     754                i,x,y,j,
     755                B->data.F64[y][x],
     756                A->data.F64[y-1][x]
     757                );
     758        if (fabs(A->data.F64[y-1][x]) > fabs(B->data.F64[y][x])) {
     759          critical_value = 2.0 * fabs(B->data.F64[y][x]);
     760        }
     761        else {
     762          critical_value = 2.0 * fabs(A->data.F64[y-1][x]);
     763        }
     764        if (critical_value < 25) { critical_value = 25; }
     765        if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&&
     766            (isfinite(B->data.F64[y][x]))&&(isfinite(A->data.F64[y-1][x]))&&
     767            (fabs(B->data.F64[y][x] - A->data.F64[y-1][x]) < critical_value)
     768            ) {
     769          solution->data.F64[i] += B->data.F64[y][x] - A->data.F64[y-1][x];
     770          XX->data.F64[i][i] += 1;
     771          XX->data.F64[i][j] += -1;
     772        }
     773      }
     774    }
     775    double max_XX = 0;
     776    double solution_V = 0;
     777    int i_peak = -1;
     778    for (int i = 0; i < numCells; i++) { // If any cells have no value of themself, set the matrix to 1.0.
     779      if (XX->data.F64[i][i] == 0.0) {
     780        XX->data.F64[i][i] = 1.0;
     781      }
     782      if (XX->data.F64[i][i] > max_XX) {
     783        max_XX = XX->data.F64[i][i];
     784        solution_V = solution->data.F64[i];
     785        i_peak = i;
     786      }
     787    }
     788    psTrace("psModules.detrend.cont",5,"fixed point: %d %g\n",
     789            i_peak,solution_V);
     790
     791    for (int i = 0; i < numCells; i++) {
     792/*        if (!((XX->data.F64[i][i] == 1.0)&& */
     793/*          (solution->data.F64[i] == 0.0))) { */
     794        solution->data.F64[i] -= solution_V;
     795        if (i != i_peak) {
     796          for (int j = 0; j < numCells; j++) {
     797            XX->data.F64[i][j] -= XX->data.F64[i_peak][j];
     798          }
     799        }
     800/*        } */
     801    }
     802    for (int i = 0; i < numCells; i++) {
     803      XX->data.F64[i_peak][i] = 0.0;
     804    }
     805    XX->data.F64[i_peak][i_peak] = 1.0;
     806   
     807   
     808#if (1)
     809    for (int i = 0; i < numCells; i++) { // print matrix A
     810      psTrace("psModules.detrend.cont",5,"A: %3d % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f",
     811              i,
     812              XX->data.F64[i][0],             XX->data.F64[i][1],             XX->data.F64[i][2],             XX->data.F64[i][3],
     813              XX->data.F64[i][4],             XX->data.F64[i][5],             XX->data.F64[i][6],             XX->data.F64[i][7],
     814              XX->data.F64[i][8],             XX->data.F64[i][9],             XX->data.F64[i][10],            XX->data.F64[i][11],
     815              XX->data.F64[i][12],            XX->data.F64[i][13],            XX->data.F64[i][14],            XX->data.F64[i][15],
     816              XX->data.F64[i][16],            XX->data.F64[i][17],            XX->data.F64[i][18],            XX->data.F64[i][19],
     817              XX->data.F64[i][20],            XX->data.F64[i][21],            XX->data.F64[i][22],            XX->data.F64[i][23],
     818              XX->data.F64[i][24],            XX->data.F64[i][25],            XX->data.F64[i][26],            XX->data.F64[i][27],
     819              XX->data.F64[i][28],            XX->data.F64[i][29],            XX->data.F64[i][30],            XX->data.F64[i][31],
     820              XX->data.F64[i][32],            XX->data.F64[i][33],            XX->data.F64[i][34],            XX->data.F64[i][35],
     821              XX->data.F64[i][36],            XX->data.F64[i][37],            XX->data.F64[i][38],            XX->data.F64[i][39],
     822              XX->data.F64[i][40],            XX->data.F64[i][41],            XX->data.F64[i][42],            XX->data.F64[i][43],
     823              XX->data.F64[i][44],            XX->data.F64[i][45],            XX->data.F64[i][46],            XX->data.F64[i][47],
     824              XX->data.F64[i][48],            XX->data.F64[i][49],            XX->data.F64[i][50],            XX->data.F64[i][51],
     825              XX->data.F64[i][52],            XX->data.F64[i][53],            XX->data.F64[i][54],            XX->data.F64[i][55],
     826              XX->data.F64[i][56],            XX->data.F64[i][57],            XX->data.F64[i][58],            XX->data.F64[i][59],
     827              XX->data.F64[i][60],            XX->data.F64[i][61],            XX->data.F64[i][62],            XX->data.F64[i][63]
     828              );
     829    }
     830
     831    for (int i = 0; i < numCells; i++) { // print vector b
     832      psTrace("psModules.detrend.cont",5,"b: %d %f",
     833              i,
     834              solution->data.F64[i]
     835              );
     836    }
     837#endif   
     838   
     839    // Solve the Ax=b equation
     840    //    psMatrixLUSolve(XX,solution);
     841    psMatrixGJSolve(XX,solution);
     842#if (1)
     843    for (int i = 0; i < numCells; i++) { // print vector b
     844      psTrace("psModules.detrend.cont",5,"x: %d %f",
     845              i,
     846              solution->data.F64[i]
     847              );
     848    }
     849#endif
     850   
     851    /* old code to remove the minimum solution value from the set, to give a "minimal set of offsets." Mathematically unnecessary. */
     852/*     double min = 99e99; */
     853/*     for (int i = 0; i < numCells; i++) { */
     854/*       if (solution->data.F64[i] < min) { */
     855/*      min = solution->data.F64[i]; */
     856/*       } */
     857/*       psTrace("psModules.detrend.cont",5,"x: %d %f %f ", */
     858/*            i, */
     859/*            solution->data.F64[i],min */
     860/*            ); */
     861/*     } */
     862/*     for (int i = 0; i < numCells; i++) { */
     863/*      if (solution->data.F64[i] != 0.0) { */
     864/*        solution->data.F64[i] -= min; */
     865/*      } */
     866/*     } */
     867
     868    // Cleanup
     869    psFree(XX);
     870    psFree(A);
     871    psFree(B);
     872    psFree(C);
     873    psFree(D);
     874
     875    // Correct cells based on the offsets calculated, and store the result in the analysis metadata.
     876    for (int i = 0; i < numCells; i++) {
     877        if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_IGNORE) {
     878            continue;
     879        }
     880        if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_TWEAK)) {
     881            continue;
     882        }
     883        pmCell *cell = chip->cells->data[i]; // Cell of interest
     884        pmReadout *ro = cell->readouts->data[0]; // Readout of interest
     885
     886        float correction = solution->data.F64[i];
     887        const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell
     888        psLogMsg("psModules.detrend", PS_LOG_DETAIL, "Correcting background of cell %s by %f",
     889                 cellName, correction);
     890        psBinaryOp(ro->image, ro->image, "-", psScalarAlloc(correction, PS_TYPE_F32));
     891        psMetadataAddF32(ro->analysis, PS_LIST_TAIL, PM_PATTERN_CELL_CORRECTION, PS_META_REPLACE,
     892                         "Pattern cell correction solution", correction);
     893    }
     894
     895    psFree(solution);
     896    psFree(meanMask);
     897
     898    return true;
     899}
     900
     901bool pmPatternContinuityApply(pmReadout *ro, psImageMaskType maskBad)
     902{
     903    PM_ASSERT_READOUT_NON_NULL(ro, false);
     904    PM_ASSERT_READOUT_IMAGE(ro, false);
     905
     906    bool mdok;                          // Status of MD lookup
     907    float corr = psMetadataLookupF32(&mdok, ro->analysis, PM_PATTERN_CELL_CORRECTION); // Correction to apply
     908    if (!mdok) {
     909        // No correction to apply
     910        return true;
     911    }
     912
     913    psImage *image = ro->image, *mask = ro->mask; // Image and mask of interest
     914    int numCols = image->numCols, numRows = image->numRows; // Size of image
     915
     916    if (!isfinite(corr)) {
     917        for (int y = 0; y < numRows; y++) {
     918            for (int x = 0; x < numCols; x++) {
     919                image->data.F32[y][x] = NAN;
     920            }
     921        }
     922        if (mask) {
     923            for (int y = 0; y < numRows; y++) {
     924                for (int x = 0; x < numCols; x++) {
     925                    mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskBad;
     926                }
     927            }
     928        }
     929    } else {
     930        for (int y = 0; y < numRows; y++) {
     931            for (int x = 0; x < numCols; x++) {
     932                image->data.F32[y][x] += corr;
     933            }
     934        }
     935    }
     936
     937    return true;
     938}
     939
     940

  • branches/eam_branches/ipp-20111122/psModules/src/detrend/pmPattern.h

    r26893 r33638  
    5454                        );
    5555
     56/// Fix the background on cells known to be troublesome
     57bool pmPatternContinuity(
     58    pmChip *chip,                       ///< Chip to correct
     59    const psVector *tweak,              ///< U8 vector indicating whether to tweak the corresponding cell
     60    psStatsOptions bgStat,              ///< Statistic to use for background measurement
     61    psStatsOptions cellStat,            ///< Statistic to use for combination of cell background measurements
     62    psImageMaskType maskVal,            ///< Mask value to use
     63    psImageMaskType maskBad,            ///< Mask value to give bad pixels
     64    int edgeWidth                       ///< Size of box to use
     65    );
     66
     67/// Apply previously measured cell pattern correction
     68bool pmPatternContinuityApply(pmReadout *ro,          ///< Readout to correct
     69                        psImageMaskType maskBad ///< Mask value to give bad pixels
     70                        );
     71
     72
    5673
    5774/// @}

  • branches/eam_branches/ipp-20111122/psModules/src/detrend/pmShutterCorrection.c

    r29004 r33638  
    805805        if (threaded) {
    806806            // wait here for the threaded jobs to finish
    807             if (!psThreadPoolWait(true)) {
     807            if (!psThreadPoolWait(true, true)) {
    808808                psError(PS_ERR_UNKNOWN, false, "Unable to apply shutter correction.");
    809809                psFree(shutterImage);

  • branches/eam_branches/ipp-20111122/psModules/src/extras/psVectorBracket.c

    r10610 r33638  
    5151        }
    5252    }
    53     // at this point, index[Nhi] >= key > index[Nlo]
    54     N = Nhi;
     53    N = (Nhi >= index->n) ? Nhi - 1 : Nhi;
     54    // at this point, index[N] >= key > index[Nlo]
    5555    while ((index->data.F32[N] >= key) && (N > Nlo)) {
    5656        N--;

  • branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmStackReject.c

    r31435 r33638  
    313313    }
    314314
    315     if (!psThreadPoolWait(false)) {
     315    if (!psThreadPoolWait(false, true)) {
    316316        psError(psErrorCodeLast(), false, "Unable to grow bad pixels.");
    317317        psFree(source);

  • branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmSubtraction.c

    r32695 r33638  
    796796    pmSubtractionStamp *stamp = job->args->data[0]; // List of stamps
    797797    pmSubtractionKernels *kernels = job->args->data[1]; // Kernels
    798     int footprint = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index
     798    int footprint = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index -- MEH - it is?
    799799
    800800    return pmSubtractionConvolveStamp(stamp, kernels, footprint);
    … …  
    832832
    833833#ifdef TESTING
     834    //MEH - index conflict or changed in past?
    834835    for (int j = 0; j < kernels->num; j++) {
    835836        if (stamp->convolutions1) {
    836837            psString convName = NULL;
    837             psStringAppend(&convName, "conv1_%03d_%03d.fits", index, j);
    838             psFits *fits = psFitsOpen(convName, "w");
     838            //psStringAppend(&convName, "conv1_%03d_%03d.fits", index, j);
     839            psStringAppend(&convName, "conv1_xxx_%03d.fits", j);
     840            psFits *fits = psFitsOpen(convName, "w");
    839841            psFree(convName);
    840842            psKernel *conv = stamp->convolutions1->data[j];
    … …  
    845847        if (stamp->convolutions2) {
    846848            psString convName = NULL;
    847             psStringAppend(&convName, "conv2_%03d_%03d.fits", index, j);
    848             psFits *fits = psFitsOpen(convName, "w");
     849            //psStringAppend(&convName, "conv2_%03d_%03d.fits", index, j);
     850            psStringAppend(&convName, "conv2_xxx_%03d.fits", j);
     851            psFits *fits = psFitsOpen(convName, "w");
    849852            psFree(convName);
    850853            psKernel *conv = stamp->convolutions2->data[j];
    … …  
    905908        }
    906909    }
    907     if (!psThreadPoolWait(true)) {
     910    if (!psThreadPoolWait(true, true)) {
    908911        psError(psErrorCodeLast(), false, "Error waiting for threads.");
    909912        return false;
    … …  
    14271430    }
    14281431
    1429     if (!psThreadPoolWait(false)) {
     1432    if (!psThreadPoolWait(false, true)) {
    14301433        psError(psErrorCodeLast(), false, "Error waiting for threads.");
    14311434        return false;

  • branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmSubtractionEquation.c

    r30622 r33638  
    958958    }
    959959
    960     if (!psThreadPoolWait(true)) {
     960    if (!psThreadPoolWait(true, true)) {
    961961        psError(psErrorCodeLast(), false, "Error waiting for threads.");
    962962        return false;

  • branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmSubtractionEquation.v0.c

    r30622 r33638  
    882882    }
    883883
    884     if (!psThreadPoolWait(true)) {
     884    if (!psThreadPoolWait(true, true)) {
    885885        psError(psErrorCodeLast(), false, "Error waiting for threads.");
    886886        return false;

  • branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmSubtractionMatch.c

    r31671 r33638  
    5353    fprintf(stderr, "    Memory in use: %zd\n", totalSize);
    5454    fprintf(stderr, "    Largest block: %ld\n", largest);
     55    //MEH -- osx may not like sbrk
    5556    fprintf(stderr, "    sbrk(): %zd\n", (size_t)sbrk(0));
    5657#endif
    … …  
    122123        PS_ASSERT_FLOAT_LESS_THAN(sysError, 1.0, false);
    123124    }
    124     if (isfinite(sysError)) {
     125    if (isfinite(skyError)) {
    125126        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(skyError, 0.0, false);
    126127    }
    … …  
    10911092    }
    10921093
    1093     if (!psThreadPoolWait(true)) {
     1094    if (!psThreadPoolWait(true, true)) {
    10941095        psError(psErrorCodeLast(), false, "Error waiting for threads.");
    10951096        psFree(models);

  • branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmSubtractionStamps.c

    r31543 r33638  
    12081208        // XXX this is somewhat arbitrary...
    12091209        if (source->psfMagErr > 0.05) continue;
    1210         if (fabs(source->psfMag - source->apMag) > 0.5) continue;
     1210        if (isfinite(source->apMag)) {
     1211            if (fabs(source->psfMag - source->apMag) > 0.5) continue;
     1212        } else if (isfinite(source->apMagRaw)) {
     1213            if (fabs(source->psfMag - source->apMagRaw) > 0.5) continue;
     1214        } else {
     1215            // XXX: Should we carry on or drop this source?
     1216            // drop it for now
     1217            continue;
     1218        }
    12111219
    12121220        if (source->modelPSF) {

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

    r32347 r33638  
    402402    assert (psf->params->n > PM_PAR_YPOS);
    403403    assert (psf->params->n > PM_PAR_XPOS);
     404
     405    if (! isfinite(Io)) {
     406        fprintf(stderr, "non-finite Io passed to PM_MODEL_PARAMS_FROM_PSF\n");
     407        return false;
     408    }
    404409
    405410    PAR[PM_PAR_SKY]  = 0.0;

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

    r32347 r33638  
    192192    psF32 z0 = PAR[PM_PAR_I0]*f1;
    193193    psF32 f0 = PAR[PM_PAR_SKY] + z0;
     194
     195    if (!isfinite(z0)) {
     196        fprintf(stderr, "z0 is not finite for %f %f %f %f %f.  Parameters: \n", X, Y, radius, z, f1);
     197        fprintf(stderr, "%f %f %f %f %f %f %f %f\n", PAR[0], PAR[1], PAR[2], PAR[3], PAR[4],
     198            PAR[5], PAR[6], PAR[7]);
     199    }
    194200
    195201    assert (isfinite(f2));

  • branches/eam_branches/ipp-20111122/psModules/src/objects/pmFootprintCullPeaks.c

    r32827 r33638  
    178178            psArray *myFP = pmFootprintsFind(subImg, threshold, 5);
    179179            if (!myFP) {
    180                 psWarning ("missing footprint?");
     180                psWarning ("missing footprint? threshold: %.f", threshold);
    181181                continue;
    182182            }
    183183            if (!myFP->n) {
    184                 psWarning ("empty footprint?");
     184                psWarning ("empty footprint? threshold: %.f", threshold);
    185185                psFree (myFP);
    186186                continue;

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

    r32347 r33638  
    144144
    145145    *Io = source->peak->rawFlux;
     146
     147#ifndef ALLOW_NONFINITE_PEAK
     148    // Gene says fail of peak !finite
     149    if (!isfinite(*Io)) return false;
     150#else
     151    // This is the way it used to be. Somtimes an infinite value Io made it's way down the pipeline
     152    // causing assertion failures
    146153    if (!isfinite(*Io) && !source->moments) return false;
    147154
    148155    *Io = source->moments->Peak;
    149156    if (!isfinite(*Io)) return false;
     157#endif
    150158
    151159    return true;

  • branches/eam_branches/ipp-20111122/psModules/src/objects/pmPCMdata.c

    r32725 r33638  
    136136            sum += value;
    137137        }
     138    }
     139
     140    if (!(sum > 0.0)) {
     141        // Crazy PSF image print out some debugging information ...
     142        fprintf(stderr, "invalid kernel sum %f found by pmPCMkernelFromPSF\n", sum);    for (int j = psf->yMin; j <= psf->yMax; j++) {
     143            fprintf(stderr, "Row %d\n", j);
     144            for (int i = psf->xMin; i <= psf->xMax; i++) {
     145                double value = source->psfImage->data.F32[y0 + j][x0 + i];
     146                fprintf(stderr, "  %d %f\n", i, value);
     147            }
     148        }
     149        fflush(stderr);
     150        // ... but avoid the asssertion two lines down by escaping
     151        goto escape;
    138152    }
    139153    assert (sum > 0.0);

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

    r30621 r33638  
    7373            continue;
    7474        }
     75        // If mask object does not exist, mark the source as bad.
     76        // We cannot proceed with it because psImageMaskPixels leaves an uncleared error code last which causes
     77        // psphot to exit with a fault.
     78        if (source->maskObj == NULL) {
     79            psTrace ("psModules.objects", 4, "source %d (%d,%d) : null maskObj\n", i, source->peak->x, source->peak->y);
     80            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
     81            continue;
     82        }
    7583
    7684        source->modelEXT = pmSourceModelGuess (source, options->type);
    … …  
    8997
    9098        // clear object mask to define valid pixels
    91         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
     99        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
    92100
    93101        // exclude the poor fits

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

    r32633 r33638  
    5757#define BLANK_HEADERS "BLANK.HEADERS"   // Name of metadata in camera configuration containing header names
    5858                                        // for putting values into a blank PHU
     59static bool pmReadoutReadXSRC(pmFPAfile *file, char * exttype, psMetadata *hduHeader, psString xsrcname, psArray *sources, long *sourceIndex);
     60static bool pmReadoutReadXFIT(pmFPAfile *file, char * exttype, psMetadata *hduHeader, psString xfitname, psArray *sources, long *sourceIndex);
     61static bool pmReadoutReadXRAD(pmFPAfile *file, pmReadout *readout, char * exttype, psMetadata *hduHeader, psString xfitname, psArray *sources, long *sourceIndex);
    5962
    6063// lookup the EXTNAME values used for table data and image header segments
    … …  
    961964        psString dataname = NULL;
    962965        psString deteffname = NULL;
    963         if (!pmSourceIOextnames(&headname, &dataname, &deteffname, NULL, NULL, NULL, file, view)) {
     966        psString xsrcname = NULL;
     967        psString xfitname = NULL;
     968        psString xradname = NULL;
     969
     970        // determine the output table format. Assume if we need to output extendend source
     971        // parameters that they may exist in the input.
     972        // XXX: Perhaps we should use different recipe values.
     973        // I.E. EXTENDED_SOURCE_ANALYSIS_READ or something like that
     974        psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, "PSPHOT");
     975        if (!status) {
     976            psError(PS_ERR_UNKNOWN, true, "missing recipe PSPHOT in config data");
     977            return false;
     978        }
     979        // if this is not TRUE, the output files only contain the psf measurements.
     980        bool XSRC_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_ANALYSIS");
     981        bool XFIT_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_FITS");
     982        bool XRAD_OUTPUT = psMetadataLookupBool(&status, recipe, "RADIAL_APERTURES");
     983
     984        if (!pmSourceIOextnames(&headname, &dataname, &deteffname,
     985                XSRC_OUTPUT ? &xsrcname : NULL,
     986                XFIT_OUTPUT ? &xfitname : NULL,
     987                XRAD_OUTPUT ? &xradname : NULL,
     988                file, view)) {
    964989            return false;
    965990        }
    … …  
    10391064            }
    10401065
     1066            long *sourceIndex = NULL;
     1067            if (XSRC_OUTPUT || XFIT_OUTPUT || XRAD_OUTPUT) {
     1068                long seq_max = -1;
     1069                for (long i = sources->n -1; i >= 0; i--) {
     1070                    pmSource *source = sources->data[i];
     1071                    if (source->seq < 0) {
     1072                        // This can happen cmf files that have been corrupted
     1073                        psError(PS_ERR_IO, true, "seq < 0 for source %ld: Suspect %s is corrupt", i, file->origname);
     1074                        return false;
     1075                    }
     1076                    if (source->seq > seq_max) {
     1077                        seq_max = source->seq;
     1078                    }
     1079                }
     1080                sourceIndex = psAlloc((seq_max + 1) * sizeof(long));
     1081                for (long i = 0; i < seq_max; i++) {
     1082                    sourceIndex[i] = -1;
     1083                }
     1084                for (long i = 0; i < sources->n; i++) {
     1085                    pmSource *source = sources->data[i];
     1086                    sourceIndex[source->seq] = i;
     1087                }
     1088            }
     1089            if (XSRC_OUTPUT && xsrcname) {
     1090                if (!pmReadoutReadXSRC(file, exttype, hdu->header, xsrcname, sources, sourceIndex)) {
     1091                    // XXX: is this an error?
     1092                    psErrorClear();
     1093                }
     1094                psFree(xsrcname);
     1095            }
     1096            if (XFIT_OUTPUT && xfitname) {
     1097                if (!pmReadoutReadXFIT(file, exttype, hdu->header, xfitname, sources, sourceIndex)) {
     1098                    // XXX: is this an error?
     1099                    psErrorClear();
     1100                }
     1101                psFree(xfitname);
     1102            }
     1103            if (XRAD_OUTPUT && xradname) {
     1104                if (!pmReadoutReadXRAD(file, readout, exttype, hdu->header, xradname, sources, sourceIndex)) {
     1105                    // XXX: is this an error?
     1106                    psErrorClear();
     1107                }
     1108                psFree(xradname);
     1109            }
     1110            psFree(sourceIndex);
     1111
    10411112            if (!pmReadoutReadDetEff(file->fits, readout, deteffname)) {
    10421113#if 0
    … …  
    11651236}
    11661237
    1167 
     1238// XXX: We might be able to macroize this and reuse for the other types
     1239
     1240static bool pmReadoutReadXSRC(pmFPAfile *file, char *exttype, psMetadata *hduHeader, psString xsrcname, psArray *sources, long *sourceIndex)
     1241{
     1242    if (!psFitsMoveExtName (file->fits, xsrcname)) {
     1243        psError(PS_ERR_UNKNOWN, false, "cannot find xsrc extension %s in %s", xsrcname, file->filename);
     1244        return false;
     1245    }
     1246
     1247    psMetadata *tableHeader = psFitsReadHeader(NULL, file->fits); // The FITS header
     1248    if (!tableHeader) psAbort("cannot read table header");
     1249
     1250    char *xtension = psMetadataLookupStr (NULL, tableHeader, "XTENSION");
     1251    if (!xtension) psAbort("cannot read table type");
     1252    if (strcmp (xtension, "BINTABLE")) {
     1253        psWarning ("no binary table in extension %s, skipping\n", xsrcname);
     1254        return false;
     1255    }
     1256
     1257    // XXX these are case-sensitive since the EXTYPE is case-sensitive
     1258    bool status = false;
     1259    if (file->type == PM_FPA_FILE_CMF) {
     1260        if (!strcmp (exttype, "PS1_SV1")) {
     1261            status  = pmSourcesRead_CMF_PS1_SV1_XSRC (file->fits, hduHeader, sources, sourceIndex);
     1262        }
     1263    }
     1264    psFree(tableHeader);
     1265    return status;
     1266}
     1267
     1268static bool pmReadoutReadXFIT(pmFPAfile *file, char *exttype, psMetadata *hduHeader, psString extname, psArray *sources, long *sourceIndex)
     1269{
     1270    if (!psFitsMoveExtName (file->fits, extname)) {
     1271        psError(PS_ERR_UNKNOWN, false, "cannot find extension %s in %s", extname, file->filename);
     1272        return false;
     1273    }
     1274
     1275    psMetadata *tableHeader = psFitsReadHeader(NULL, file->fits); // The FITS header
     1276    if (!tableHeader) psAbort("cannot read table header");
     1277
     1278    char *xtension = psMetadataLookupStr (NULL, tableHeader, "XTENSION");
     1279    if (!xtension) psAbort("cannot read table type");
     1280    if (strcmp (xtension, "BINTABLE")) {
     1281        psWarning ("no binary table in extension %s, skipping\n", extname);
     1282        return false;
     1283    }
     1284
     1285    // XXX these are case-sensitive since the EXTYPE is case-sensitive
     1286    bool status = false;
     1287    if (file->type == PM_FPA_FILE_CMF) {
     1288        if (!strcmp (exttype, "PS1_SV1")) {
     1289            status  = pmSourcesRead_CMF_PS1_SV1_XFIT (file->fits, hduHeader, sources, sourceIndex);
     1290        }
     1291    }
     1292    psFree(tableHeader);
     1293    return status;
     1294}
     1295static bool pmReadoutReadXRAD(pmFPAfile *file, pmReadout *readout, char *exttype, psMetadata *hduHeader, psString extname, psArray *sources, long *sourceIndex)
     1296{
     1297    if (!psFitsMoveExtName (file->fits, extname)) {
     1298        psError(PS_ERR_UNKNOWN, false, "cannot find extension %s in %s", extname, file->filename);
     1299        return false;
     1300    }
     1301
     1302    psMetadata *tableHeader = psFitsReadHeader(NULL, file->fits); // The FITS header
     1303    if (!tableHeader) psAbort("cannot read table header");
     1304
     1305    char *xtension = psMetadataLookupStr (NULL, tableHeader, "XTENSION");
     1306    if (!xtension) psAbort("cannot read table type");
     1307    if (strcmp (xtension, "BINTABLE")) {
     1308        psWarning ("no binary table in extension %s, skipping\n", extname);
     1309        return false;
     1310    }
     1311
     1312    // XXX these are case-sensitive since the EXTYPE is case-sensitive
     1313    bool status = false;
     1314    if (file->type == PM_FPA_FILE_CMF) {
     1315        if (!strcmp (exttype, "PS1_SV1")) {
     1316            status  = pmSourcesRead_CMF_PS1_SV1_XRAD (file->fits, readout, hduHeader, sources, sourceIndex);
     1317        }
     1318    }
     1319    psFree(tableHeader);
     1320    return status;
     1321}

  • branches/eam_branches/ipp-20111122/psModules/src/objects/pmSourceIO.h

    r32633 r33638  
    9393psArray *pmSourcesRead_CMF_PS1_V4 (psFits *fits, psMetadata *header);
    9494psArray *pmSourcesRead_CMF_PS1_SV1 (psFits *fits, psMetadata *header);
     95bool pmSourcesRead_CMF_PS1_SV1_XSRC (psFits *fits, psMetadata *header, psArray *sources, long *);
     96bool pmSourcesRead_CMF_PS1_SV1_XFIT (psFits *fits, psMetadata *header, psArray *sources, long *);
     97bool pmSourcesRead_CMF_PS1_SV1_XRAD (psFits *fits, pmReadout *readout, psMetadata *header, psArray *sources, long *);
    9598psArray *pmSourcesRead_CMF_PS1_DV1 (psFits *fits, psMetadata *header);
    9699psArray *pmSourcesRead_CMF_PS1_DV2 (psFits *fits, psMetadata *header);

  • branches/eam_branches/ipp-20111122/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c

    r32347 r33638  
    539539}
    540540
     541bool pmSourcesRead_CMF_PS1_SV1_XSRC(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex)
     542{
     543    PS_ASSERT_PTR_NON_NULL(fits, false);
     544    PS_ASSERT_PTR_NON_NULL(sources, false);
     545
     546    bool status;
     547    long numSources = psFitsTableSize(fits); // Number of sources in table
     548    if (numSources == 0) {
     549        psError(psErrorCodeLast(), false, "XSRC Table contains no entries\n");
     550        return false;
     551    }
     552
     553    // petrosian mags are not saved, we need to calculate fluxes. For this we need exptime and zero point
     554    float zeropt = psMetadataLookupF32(&status, hduHeader, "FPA.ZP");
     555    float exptime = psMetadataLookupF32(&status, hduHeader, "EXPTIME");
     556    float magOffset = zeropt + 2.5*log10(exptime);
     557
     558    for (long i = 0; i < numSources; i++) {
     559        psMetadata *row = psFitsReadTableRow(fits, i); // Table row
     560        if (!row) {
     561            psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i);
     562            psFree(row);
     563            return false;
     564        }
     565        // Find the source with this sequence number.
     566        // XXX: I am assuming that sources is sorted in order of seq
     567        long seq = psMetadataLookupU32 (&status, row, "IPP_IDET");
     568        pmSource *source = NULL;
     569#ifndef ASSUME_SORTED
     570        long j = seq < sources->n ? seq : sources->n - 1;
     571        for (; j >= 0; j--) {
     572            source = sources->data[j];
     573            if (source->seq == seq) {
     574                break;
     575            }
     576        }
     577#else
     578        long j = sourceIndex[seq];
     579        psAssert(j >= 0 && j < sources->n, "invalid sourceIndex");
     580        source = sources->data[j];
     581#endif
     582        if (!source) {
     583            psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq);
     584            psFree(row);
     585            return false;
     586        }
     587
     588        if (!source->extpars) {
     589            source->extpars = pmSourceExtendedParsAlloc ();
     590        }
     591        pmSourceExtendedPars *extpars = source->extpars;
     592
     593        // Assume that X_EXT Y_EXT and sigmas match the psf src so skip
     594
     595        // We don't have enough information to calculate the major and minor axis. Set major to 1. Should we scale this by
     596        // psf size or something?
     597        extpars->axes.major = 1.0;
     598        extpars->axes.minor = extpars->axes.major * psMetadataLookupF32(&status, row, "F25_ARATIO");
     599        extpars->axes.theta = psMetadataLookupF32(&status, row, "F25_THETA");
     600
     601        float mag = psMetadataLookupF32(&status, row, "PETRO_MAG");
     602        float magErr = psMetadataLookupF32(&status, row, "PETRO_MAG_ERR");
     603        if (isfinite(mag)) {
     604            extpars->petrosianFlux    = pow(10., (magOffset - mag) / 2.5);
     605            if (isfinite(magErr)) {
     606                extpars->petrosianFluxErr = extpars->petrosianFlux / magErr;
     607            }
     608        }
     609
     610        extpars->petrosianRadius   = psMetadataLookupF32(&status, row, "PETRO_RADIUS");
     611        extpars->petrosianRadiusErr= psMetadataLookupF32(&status, row, "PETRO_RADIUS_ERR");
     612        extpars->petrosianR50      = psMetadataLookupF32(&status, row, "PETRO_RADIUS_50");
     613        extpars->petrosianR50Err   = psMetadataLookupF32(&status, row, "PETRO_RADIUS_50_ERR");
     614        extpars->petrosianR90      = psMetadataLookupF32(&status, row, "PETRO_RADIUS_90");
     615        extpars->petrosianR90Err   = psMetadataLookupF32(&status, row, "PETRO_RADIUS_90_ERR");
     616        extpars->petrosianFill     = psMetadataLookupF32(&status, row, "PETRO_FILL");
     617
     618        psVector *radSB   = psMetadataLookupVector(&status, row, "PROF_SB");
     619        psVector *radFlux = psMetadataLookupVector(&status, row, "PROF_FLUX");
     620        psVector *radFill = psMetadataLookupVector(&status, row, "PROF_FILL");
     621
     622        if (radSB && radSB->n > 0) {
     623            extpars->radProfile = pmSourceRadialProfileAlloc();
     624            extpars->radProfile->binSB   = psMemIncrRefCounter(radSB);
     625            extpars->radProfile->binSum   = psMemIncrRefCounter(radFlux);
     626            extpars->radProfile->binFill = psMemIncrRefCounter(radFill);
     627        }
     628
     629        psFree(row);
     630    }
     631
     632    return true;
     633}
     634
    541635// XXX this layout is still the same as PS1_DEV_1
    542636bool pmSourcesWrite_CMF_PS1_SV1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname)
    … …  
    684778    psFree (outhead);
    685779    psFree (table);
     780    return true;
     781}
     782
     783bool pmSourcesRead_CMF_PS1_SV1_XFIT(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex)
     784{
     785    PS_ASSERT_PTR_NON_NULL(fits, false);
     786    PS_ASSERT_PTR_NON_NULL(sources, false);
     787
     788    bool status;
     789    long numSources = psFitsTableSize(fits); // Number of sources in table
     790    if (numSources == 0) {
     791        psError(psErrorCodeLast(), false, "XFIT Table contains no entries\n");
     792        return false;
     793    }
     794
     795    for (long i = 0; i < numSources; i++) {
     796        psMetadata *row = psFitsReadTableRow(fits, i); // Table row
     797        if (!row) {
     798            psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i);
     799            psFree(row);
     800            return false;
     801        }
     802        // Find the source with this sequence number.
     803        // XXX: I am assuming that sources is sorted in order of seq.
     804        long seq = psMetadataLookupU32 (&status, row, "IPP_IDET");
     805        long j = seq < sources->n ? seq : sources->n - 1;
     806        pmSource *source = NULL;
     807        for (; j >= 0; j--) {
     808            source = sources->data[j];
     809            if (source->seq == seq) {
     810                break;
     811            }
     812        }
     813        if (!source) {
     814            psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq);
     815            psFree(row);
     816            return false;
     817        }
     818        if (!source->modelFits) {
     819            // XXX: where to find the number of models to expect?
     820            source->modelFits = psArrayAllocEmpty(5);
     821        }
     822        psString modelName = psMetadataLookupStr(&status, row, "MODEL_TYPE");
     823        if (!modelName) {
     824            psError(PS_ERR_UNKNOWN, true, "Failed to find model name for row %ld\n", i);
     825            psFree(row);
     826            return false;
     827        }
     828        pmModelType modelType = pmModelClassGetType(modelName);
     829        if (modelType < 0) {
     830            psError(PS_ERR_UNKNOWN, true, "Failed to find model type for %s\n", modelName);
     831            psFree(row);
     832            return false;
     833        }
     834        pmModel *model = pmModelAlloc(modelType);
     835
     836        psF32 *PAR = model->params->data.F32;
     837        psF32 *dPAR = model->dparams->data.F32;
     838
     839        PAR[PM_PAR_XPOS] = psMetadataLookupF32(&status, row, "X_EXT");
     840        PAR[PM_PAR_YPOS] = psMetadataLookupF32(&status, row, "Y_EXT");
     841        dPAR[PM_PAR_XPOS] = psMetadataLookupF32(&status, row, "X_EXT_SIG");
     842        dPAR[PM_PAR_YPOS] = psMetadataLookupF32(&status, row, "Y_EXT_SIG");
     843
     844        model->mag = psMetadataLookupF32(&status, row, "EXT_INST_MAG");
     845        model->magErr = psMetadataLookupF32(&status, row, "EXT_INST_MAG_SIG");
     846
     847        psEllipseAxes axes;
     848        axes.major = psMetadataLookupF32(&status, row, "EXT_WIDTH_MAJ");
     849        axes.minor = psMetadataLookupF32(&status, row, "EXT_WIDTH_MIN");
     850        axes.theta = psMetadataLookupF32(&status, row, "EXT_THETA");
     851        if (!pmPSF_AxesToModel(PAR, axes, modelType)) {
     852            // Do we need to fail here or can this happen?
     853            psError(PS_ERR_UNKNOWN, false, "Failed to convert psf axes to model");
     854            psFree(model);
     855            psFree(row);
     856            return false;
     857        }
     858        // XXX: clean this up
     859        if (model->params->n > 7) {
     860            PAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07");
     861        }
     862        // read the covariance matrix
     863        int nparams = model->params->n;
     864        psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32);
     865        for (int y = 0; y < nparams; y++) {
     866            for (int x = 0; x < nparams; x++) {
     867                char name[64];
     868                snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x);
     869                covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name);
     870            }
     871        }
     872        model->covar = covar;
     873
     874        psArrayAdd(source->modelFits, 1, model);
     875        psFree(model);
     876
     877        psFree(row);
     878    }
     879
    686880    return true;
    687881}
    … …  
    8311025    return true;
    8321026}
     1027
     1028bool pmSourcesRead_CMF_PS1_SV1_XRAD(psFits *fits, pmReadout *readout, psMetadata *hduHeader, psArray *sources, long *sourceIndex)
     1029{
     1030    PS_ASSERT_PTR_NON_NULL(fits, false);
     1031    PS_ASSERT_PTR_NON_NULL(sources, false);
     1032
     1033    bool status;
     1034    long numSources = psFitsTableSize(fits); // Number of sources in table
     1035    if (numSources == 0) {
     1036        psError(psErrorCodeLast(), false, "XRAD Table contains no entries\n");
     1037        return false;
     1038    }
     1039
     1040    long       seq_first = -1;
     1041    long       seq_last = -1;
     1042    psVector   *fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32);
     1043    long       max_entries = -1;
     1044    long       num_entries = -1;
     1045
     1046    for (long i = 0; i < numSources; i++) {
     1047        psMetadata *row = psFitsReadTableRow(fits, i); // Table row
     1048        if (!row) {
     1049            psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i);
     1050            psFree(row);
     1051            return false;
     1052        }
     1053        // Find the source with this sequence number.
     1054        // XXX: I am assuming that sources is sorted in order of seq.
     1055        long seq = psMetadataLookupU32 (&status, row, "IPP_IDET");
     1056        long j = seq < sources->n ? seq : sources->n - 1;
     1057        pmSource *source = NULL;
     1058        for (; j >= 0; j--) {
     1059            source = sources->data[j];
     1060            if (source->seq == seq) {
     1061                break;
     1062            }
     1063        }
     1064        if (!source) {
     1065            psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq);
     1066            psFree(row);
     1067            return false;
     1068        }
     1069        if (seq_first == -1) {
     1070            seq_first = seq;
     1071        }
     1072        if (seq == seq_first) {
     1073            psF32 value = psMetadataLookupF32(&status, row, "PSF_FWHM");
     1074            psVectorAppend(fwhmValues, value);
     1075        }
     1076        if (seq == seq_last) {
     1077            num_entries++;
     1078        } else {
     1079            num_entries = 1;
     1080            seq_last = seq;
     1081        }
     1082        if (num_entries > max_entries) {
     1083            max_entries = num_entries;
     1084        }
     1085
     1086        if (!source->radialAper) {
     1087            // XXX: where to find the number of models to expect?
     1088            source->radialAper = psArrayAllocEmpty(5);
     1089        }
     1090        pmSourceRadialApertures *radialAper = pmSourceRadialAperturesAlloc();
     1091
     1092        radialAper->flux = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX"));
     1093        radialAper->fluxStdev = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX_STDEV"));
     1094        radialAper->fluxErr = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX_ERR"));
     1095        radialAper->fill = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FILL"));
     1096
     1097        psArrayAdd(source->radialAper, 1, radialAper);
     1098
     1099        psFree(radialAper);
     1100        psFree(row);
     1101    }
     1102
     1103    // check for consistency between the length of fwhmValues and the maximum number of entries for each row
     1104    if (fwhmValues->n != max_entries) {
     1105        psError(PS_ERR_PROGRAMMING, true, "number of PSF_FWHM values found %ld does not match expected number: %ld\n",
     1106            fwhmValues->n, max_entries);
     1107        psAssert(0, "fixme");
     1108    }
     1109
     1110    if (!readout->analysis) {
     1111        readout->analysis = psMetadataAlloc();
     1112    }
     1113
     1114    psMetadataAddVector(readout->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues);
     1115    psFree(fwhmValues);
     1116
     1117    return true;
     1118}

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

    r33085 r33638  
    569569        }
    570570    }
    571 
    572571    return (true);
    573572}
Note: See TracChangeset for help on using the changeset viewer.