IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 33243 for trunk/psModules


Ignore:
Timestamp:
Feb 10, 2012, 6:07:30 PM (14 years ago)
Author:
watersc1
Message:

Implementation of PATTERN.CONTINUITY correction, which replaces the PATTERN.CELL correction and ensures that the cell-to-cell variations are minimized, and that background fitting code has a smooth background to work with. Includes ippconfig files to enable this by default across all of GPC1. Tested, and memory leaks fixed.

Location:
trunk/psModules/src/detrend
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/detrend/pmPattern.c

    r32970 r33243  
    146146        // Store the results we found for this row.
    147147        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       
    148154#endif
    149155        memcpy(corr->data.F64[y], poly->coeff, (order + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
    … …  
    161167        for (int x = 0; x < numCols; x++) {
    162168            image->data.F32[y][x] -= solution->data.F32[x];
     169            psTrace("pattern",5,"A: %d %d %g\n",x,y,solution->data.F32[x]);
    163170        }
    164171        psFree(solution);
    … …  
    177184    // Fit the trend of the constant term, producing the y-axis global trend
    178185    psStatsInit(clip);
    179     psPolynomial1D *yaxisPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); // Polynomial to fit.
     186    psPolynomial1D *yaxisPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); // Polynomial to fit.
    180187    if (!psVectorClipFitPolynomial1D(yaxisPoly,clip,yaxisMask,0xFF,yaxisData, NULL, yaxisIndices)) {
    181188      psWarning("Unable to fit polynomial to y-axis trend");
    … …  
    209216            image->data.F32[y][x] += solution->data.F32[y];
    210217            corr->data.F64[y][0]  -= solution->data.F32[y];
     218            psTrace("pattern",5,"B: %d %d %g\n",x,y,solution->data.F32[x]);
    211219          }
    212220        }
    … …  
    218226    // We can use the same mask vector, as the same rows failed the row-fit earlier.
    219227    psStatsInit(clip);
    220     psPolynomial1D *xaxisPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); // Polynomial to fit.
     228    psPolynomial1D *xaxisPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); // Polynomial to fit.
    221229    if (!psVectorClipFitPolynomial1D(xaxisPoly,clip,yaxisMask,0xFF,xaxisData, NULL, yaxisIndices)) {
    222230      psWarning("Unable to fit polynomial to x-axis trend");
    … …  
    233241          for (int x = 0; x < numCols; x++) {
    234242            image->data.F32[y][x] += solution->data.F32[y] * indices->data.F32[x];
    235             corr->data.F64[y][0]  -= 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]);
    236245          }
    237246        }
    … …  
    500509
    501510
     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   
     568    for (int i = 0; i < numCells; i++) {
     569        if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_IGNORE) {
     570            continue;
     571        }
     572        pmCell *cell = chip->cells->data[i]; // Cell of interest
     573        pmReadout *ro = cell->readouts->data[0]; // Readout of interest
     574
     575        psStatsInit(bgStats);
     576
     577        // Convert cell iterator i into an xy coordinate on the grid of cells
     578        int y = (i % 8);
     579        int x = (i - y) / 8;
     580       
     581        for (int j = 0; j < 4; j++) {
     582          if (j == 0) {  // Region B
     583            region = psRegionSet(0,ro->image->numCols,
     584                                 0,edgeWidth);
     585          }
     586          else if (j == 1) { // Region A
     587            region = psRegionSet(0,ro->image->numCols,
     588                                 ro->image->numRows - edgeWidth,ro->image->numRows);
     589          }
     590          else if (j == 2) { // Region D
     591            region = psRegionSet(0,edgeWidth,
     592                                 0,ro->image->numRows);
     593          }
     594          else if (j == 3) { // Region C
     595            region = psRegionSet(ro->image->numCols - edgeWidth,ro->image->numCols,
     596                                 0,ro->image->numRows);
     597          }
     598          psImage *subset  = psImageSubset(ro->image,region);
     599          psImage *submask = psImageSubset(ro->mask,region);
     600
     601          if (!psImageBackground(bgStats, NULL, subset, submask, maskVal, rng)) {
     602            psWarning("Unable to measure background for cell %d on edge %d\n", i, j);
     603            psErrorClear();
     604            meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PM_PATTERN_ERROR;
     605            if (j == 0)      { B->data.F64[y][x] = NAN; }
     606            else if (j == 1) { A->data.F64[y][x] = NAN; }
     607            else if (j == 2) { C->data.F64[y][x] = NAN; }
     608            else if (j == 3) { D->data.F64[y][x] = NAN; }
     609            psFree(subset);
     610            psFree(submask);
     611            continue; // Move on to next edge, as only part of this cell may be a problem
     612          }
     613
     614          // If the returned value is zero, assume something is wrong.  Do I still need this?
     615          if (psStatsGetValue(bgStats,bgStat) < 1e-6) {
     616            if (j == 0)      { B->data.F64[y][x] = NAN; }
     617            else if (j == 1) { A->data.F64[y][x] = NAN; }
     618            else if (j == 2) { C->data.F64[y][x] = NAN; }
     619            else if (j == 3) { D->data.F64[y][x] = NAN; }
     620          }
     621          // If we have an error for this cell/edge, make sure we mask the value
     622          if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_ERROR) {
     623            if (j == 0)      { B->data.F64[y][x] = NAN; }
     624            else if (j == 1) { A->data.F64[y][x] = NAN; }
     625            else if (j == 2) { C->data.F64[y][x] = NAN; }
     626            else if (j == 3) { D->data.F64[y][x] = NAN; }
     627          }
     628          else { // Set the value to match what we got from the edge box.
     629            if (j == 0)      { B->data.F64[y][x] = psStatsGetValue(bgStats,bgStat); }
     630            else if (j == 1) { A->data.F64[y][x] = psStatsGetValue(bgStats,bgStat); }
     631            else if (j == 2) { C->data.F64[y][x] = psStatsGetValue(bgStats,bgStat); }
     632            else if (j == 3) { D->data.F64[y][x] = psStatsGetValue(bgStats,bgStat); }
     633          }
     634          psFree(subset);
     635          psFree(submask);
     636
     637        }
     638        psTrace("psModules.detrend.cont",5, "OTA: %d (%d %d) A: %f B: %f C: %f D: %f",
     639                i,x,y,
     640                A->data.F32[y][x],B->data.F32[y][x],C->data.F32[y][x],D->data.F32[y][x]);               
     641    }
     642    psFree(bgStats);
     643    psFree(rng);
     644
     645    // We've now allocated all the edge values, so we can now minimize the offsets.
     646    // This involves solving the equation A x = b, where
     647    // A is the (64x64 for GPC1) matrix containing the edges that match for each cell
     648    // x is the solution vector
     649    // b is the combination of offsets across each cell boundary for each cell.
     650    // Below "XX" is used as the matrix A, and "solution" is used as both b and x
     651    //   (due to the way psMatrixLUSolve operates).
     652    psVector *solution = psVectorAlloc(64,PS_TYPE_F64);
     653    psImage  *XX       = psImageAlloc(64,64,PS_TYPE_F64);
     654    psVectorInit(solution,0.0);
     655    psImageInit(XX,0.0);
     656   
     657    for (int i = 0; i < numCells; i++) {
     658      // Accumulate all the possible edge differences we can for this cell.
     659      // As we do so, make a note of the correlations by incrementing the element of the matrix.
     660      int y = (i % 8);
     661      int x = (i - y) / 8;
     662      int j;
     663      if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_IGNORE) {
     664        continue;
     665      }
     666      if (x + 1 < 8) {  // We have a neighbor adjacent in the +x direction
     667        j = 8 * (x + 1) + y; // Determine that neighbor's index
     668        psTrace("psModules.detrend.cont",5,"CmD %d %d %d %d %g %g", // diagnostic
     669                i,x,y,j,
     670                C->data.F64[y][x],
     671                D->data.F64[y][x+1]
     672                );
     673        if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&&  // If there are no errors with the neighbor,
     674            (isfinite(C->data.F64[y][x]))&&(isfinite(D->data.F64[y][x+1]))) {    // and all edges have valid values:
     675          solution->data.F64[i] += C->data.F64[y][x] - D->data.F64[y][x+1];     // Take the difference
     676          XX->data.F64[i][i] += 1;                                              // increment our relation with ourself
     677          XX->data.F64[i][j] += -1;                                             // decrement our relation with the neighbor
     678        }
     679      }
     680      if (x - 1 > -1) { // etc.
     681        j = 8 * (x - 1) + y;
     682        psTrace("psModules.detrend.cont",5,"DmC %d %d %d %d %g %g",
     683                i,x,y,j,
     684                D->data.F64[y][x],
     685                C->data.F64[y][x-1]
     686                );
     687
     688        if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&&
     689            (isfinite(D->data.F64[y][x]))&&(isfinite(C->data.F64[y][x-1]))) {
     690          solution->data.F64[i] += D->data.F64[y][x] - C->data.F64[y][x-1];
     691          XX->data.F64[i][i] += 1;
     692          XX->data.F64[i][j] += -1;
     693        }
     694      }
     695      if (y + 1 < 8) {
     696        j = 8 * x + (y + 1);
     697        psTrace("psModules.detrend.cont",5,"AmB %d %d %d %d %g %g",
     698                i,x,y,j,
     699                A->data.F64[y][x],
     700                B->data.F64[y+1][x]
     701                );
     702
     703        if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&&
     704            (isfinite(A->data.F64[y][x]))&&(isfinite(B->data.F64[y+1][x]))) {
     705          solution->data.F64[i] += A->data.F64[y][x] - B->data.F64[y+1][x];
     706          XX->data.F64[i][i] += 1;
     707          XX->data.F64[i][j] += -1;
     708        }
     709      }
     710      if (y - 1 > -1) {
     711        j = 8 * x +  (y - 1);
     712        psTrace("psModules.detrend.cont",5,"BmA %d %d %d %d %g %g",
     713                i,x,y,j,
     714                B->data.F64[y][x],
     715                A->data.F64[y-1][x]
     716                );
     717
     718        if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&&
     719            (isfinite(B->data.F64[y][x]))&&(isfinite(A->data.F64[y-1][x]))) {
     720          solution->data.F64[i] += B->data.F64[y][x] - A->data.F64[y-1][x];
     721          XX->data.F64[i][i] += 1;
     722          XX->data.F64[i][j] += -1;
     723        }
     724      }
     725    }
     726    for (int i = 0; i < numCells; i++) { // If any cells have no value of themself, set the matrix to 1.0.
     727      if (XX->data.F64[i][i] == 0.0) {
     728        XX->data.F64[i][i] = 1.0;
     729      }
     730    }
     731#if (0)
     732    for (int i = 0; i < numCells; i++) { // print matrix A
     733      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",
     734              i,
     735              XX->data.F64[i][0],             XX->data.F64[i][1],             XX->data.F64[i][2],             XX->data.F64[i][3],
     736              XX->data.F64[i][4],             XX->data.F64[i][5],             XX->data.F64[i][6],             XX->data.F64[i][7],
     737              XX->data.F64[i][8],             XX->data.F64[i][9],             XX->data.F64[i][10],            XX->data.F64[i][11],
     738              XX->data.F64[i][12],            XX->data.F64[i][13],            XX->data.F64[i][14],            XX->data.F64[i][15],
     739              XX->data.F64[i][16],            XX->data.F64[i][17],            XX->data.F64[i][18],            XX->data.F64[i][19],
     740              XX->data.F64[i][20],            XX->data.F64[i][21],            XX->data.F64[i][22],            XX->data.F64[i][23],
     741              XX->data.F64[i][24],            XX->data.F64[i][25],            XX->data.F64[i][26],            XX->data.F64[i][27],
     742              XX->data.F64[i][28],            XX->data.F64[i][29],            XX->data.F64[i][30],            XX->data.F64[i][31],
     743              XX->data.F64[i][32],            XX->data.F64[i][33],            XX->data.F64[i][34],            XX->data.F64[i][35],
     744              XX->data.F64[i][36],            XX->data.F64[i][37],            XX->data.F64[i][38],            XX->data.F64[i][39],
     745              XX->data.F64[i][40],            XX->data.F64[i][41],            XX->data.F64[i][42],            XX->data.F64[i][43],
     746              XX->data.F64[i][44],            XX->data.F64[i][45],            XX->data.F64[i][46],            XX->data.F64[i][47],
     747              XX->data.F64[i][48],            XX->data.F64[i][49],            XX->data.F64[i][50],            XX->data.F64[i][51],
     748              XX->data.F64[i][52],            XX->data.F64[i][53],            XX->data.F64[i][54],            XX->data.F64[i][55],
     749              XX->data.F64[i][56],            XX->data.F64[i][57],            XX->data.F64[i][58],            XX->data.F64[i][59],
     750              XX->data.F64[i][60],            XX->data.F64[i][61],            XX->data.F64[i][62],            XX->data.F64[i][63]
     751              );
     752    }
     753
     754    for (int i = 0; i < numCells; i++) { // print vector b
     755      psTrace("psModules.detrend.cont",5,"b: %d %f",
     756              i,
     757              solution->data.F64[i]
     758              );
     759    }
     760#endif   
     761   
     762    // Solve the Ax=b equation
     763    psMatrixLUSolve(XX,solution);
     764
     765    /* old code to remove the minimum solution value from the set, to give a "minimal set of offsets." Mathematically unnecessary. */
     766/*     double min = 99e99; */
     767/*     for (int i = 0; i < numCells; i++) { */
     768/*       if (solution->data.F64[i] < min) { */
     769/*      min = solution->data.F64[i]; */
     770/*       } */
     771/*       psTrace("psModules.detrend.cont",5,"x: %d %f %f ", */
     772/*            i, */
     773/*            solution->data.F64[i],min */
     774/*            ); */
     775/*     } */
     776/*     for (int i = 0; i < numCells; i++) { */
     777/*      if (solution->data.F64[i] != 0.0) { */
     778/*        solution->data.F64[i] -= min; */
     779/*      } */
     780/*     } */
     781
     782    // Cleanup
     783    psFree(XX);
     784    psFree(A);
     785    psFree(B);
     786    psFree(C);
     787    psFree(D);
     788
     789    // Correct cells based on the offsets calculated, and store the result in the analysis metadata.
     790    for (int i = 0; i < numCells; i++) {
     791        if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_IGNORE) {
     792            continue;
     793        }
     794        if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_TWEAK)) {
     795            continue;
     796        }
     797        pmCell *cell = chip->cells->data[i]; // Cell of interest
     798        pmReadout *ro = cell->readouts->data[0]; // Readout of interest
     799
     800        float correction = solution->data.F64[i];
     801        const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell
     802        psLogMsg("psModules.detrend", PS_LOG_DETAIL, "Correcting background of cell %s by %f",
     803                 cellName, correction);
     804        psBinaryOp(ro->image, ro->image, "-", psScalarAlloc(correction, PS_TYPE_F32));
     805        psMetadataAddF32(ro->analysis, PS_LIST_TAIL, PM_PATTERN_CELL_CORRECTION, PS_META_REPLACE,
     806                         "Pattern cell correction solution", correction);
     807    }
     808
     809    psFree(solution);
     810    psFree(meanMask);
     811
     812    return true;
     813}
     814
     815bool pmPatternContinuityApply(pmReadout *ro, psImageMaskType maskBad)
     816{
     817    PM_ASSERT_READOUT_NON_NULL(ro, false);
     818    PM_ASSERT_READOUT_IMAGE(ro, false);
     819
     820    bool mdok;                          // Status of MD lookup
     821    float corr = psMetadataLookupF32(&mdok, ro->analysis, PM_PATTERN_CELL_CORRECTION); // Correction to apply
     822    if (!mdok) {
     823        // No correction to apply
     824        return true;
     825    }
     826
     827    psImage *image = ro->image, *mask = ro->mask; // Image and mask of interest
     828    int numCols = image->numCols, numRows = image->numRows; // Size of image
     829
     830    if (!isfinite(corr)) {
     831        for (int y = 0; y < numRows; y++) {
     832            for (int x = 0; x < numCols; x++) {
     833                image->data.F32[y][x] = NAN;
     834            }
     835        }
     836        if (mask) {
     837            for (int y = 0; y < numRows; y++) {
     838                for (int x = 0; x < numCols; x++) {
     839                    mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskBad;
     840                }
     841            }
     842        }
     843    } else {
     844        for (int y = 0; y < numRows; y++) {
     845            for (int x = 0; x < numCols; x++) {
     846                image->data.F32[y][x] += corr;
     847            }
     848        }
     849    }
     850
     851    return true;
     852}
     853
     854

  • trunk/psModules/src/detrend/pmPattern.h

    r26893 r33243  
    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/// @}
Note: See TracChangeset for help on using the changeset viewer.