IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 1, 2013, 5:01:08 PM (13 years ago)
Author:
watersc1
Message:

Merging changes that implement FPA level background continuity.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r33340 r35081  
    938938}
    939939
    940 
     940// This should reuse code more efficiently
     941bool pmPatternContinuityBackground(pmFPAfile *in, pmFPAfile *out,  psStatsOptions bgStat, psStatsOptions cellStat,
     942                                   psImageMaskType maskVal, psImageMaskType maskBad, int edgeWidth)
     943{
     944    PS_ASSERT_PTR_NON_NULL(in, false);
     945    PS_ASSERT_PTR_NON_NULL(out, false);
     946
     947    int numChips = out->fpa->chips->n;            // Number of cells
     948
     949    psVector *meanMask = psVectorAlloc(numChips, PS_TYPE_VECTOR_MASK); // Mask for means
     950    psVectorInit(meanMask, 0);
     951
     952    // Mask bits
     953    enum {
     954        PM_PATTERN_IGNORE = 0x01,       // Ignore this cell
     955        PM_PATTERN_TWEAK  = 0x02,       // Tweak this cell
     956        PM_PATTERN_ERROR  = 0x04,       // Error in calculating background
     957        PM_PATTERN_ALL    = 0xFF,       // All causes
     958    };
     959
     960    // Measure mean of each cell edge, and use that to determine the cell offsets.
     961    psStatsOptions stat = cellStat;     // Define which statistic to use.
     962   
     963    psStats *bgStats = psStatsAlloc(stat); // Statistics on background
     964    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
     965
     966    psRegion region = {0,0,0,0};
     967
     968    /* These images hold the edge data for the OTA structure.  */
     969    psImage *A = psImageAlloc(8,8,PS_TYPE_F64); // Top edge
     970    psImage *B = psImageAlloc(8,8,PS_TYPE_F64); // Bottom edge
     971    psImage *C = psImageAlloc(8,8,PS_TYPE_F64); // Right edge
     972    psImage *D = psImageAlloc(8,8,PS_TYPE_F64); // Left edge
     973    psImageInit(A,0.0);
     974    psImageInit(B,0.0);
     975    psImageInit(C,0.0);
     976    psImageInit(D,0.0);
     977
     978    // Corners don't exist.
     979    A->data.F64[0][0] = NAN;
     980    B->data.F64[0][0] = NAN;
     981    C->data.F64[0][0] = NAN;
     982    D->data.F64[0][0] = NAN;
     983    A->data.F64[7][0] = NAN;
     984    B->data.F64[7][0] = NAN;
     985    C->data.F64[7][0] = NAN;
     986    D->data.F64[7][0] = NAN;
     987    A->data.F64[0][7] = NAN;
     988    B->data.F64[0][7] = NAN;
     989    C->data.F64[0][7] = NAN;
     990    D->data.F64[0][7] = NAN;
     991    A->data.F64[7][7] = NAN;
     992    B->data.F64[7][7] = NAN;
     993    C->data.F64[7][7] = NAN;
     994    D->data.F64[7][7] = NAN;
     995   
     996    for (int i = 0; i < numChips; i++) {
     997      pmChip *chip = out->fpa->chips->data[i];
     998      pmCell *cell = chip->cells->data[0]; // Cell of interest
     999      pmReadout *ro = cell->readouts->data[0]; // Readout of interest
     1000
     1001      psStatsInit(bgStats);
     1002
     1003      // Convert cell iterator i into an xy coordinate on the grid of cells
     1004      // This is wrong for chips
     1005      const char *chipName = psMetadataLookupStr(NULL,chip->concepts, "CHIP.NAME");
     1006      int xParity          = psMetadataLookupS16(NULL,chip->concepts, "CHIP.XPARITY");
     1007      int yParity          = psMetadataLookupS16(NULL,chip->concepts, "CHIP.YPARITY");
     1008      int x = chipName[2] - '0';
     1009      int y = chipName[3] - '0';
     1010
     1011     
     1012      for (int j = 0; j < 4; j++) {
     1013        if (j == 0) {  // Region B
     1014          region = psRegionSet(0,ro->image->numCols,
     1015                               0,edgeWidth);
     1016        }
     1017        else if (j == 1) { // Region A
     1018          region = psRegionSet(0,ro->image->numCols,
     1019                               ro->image->numRows - edgeWidth,ro->image->numRows);
     1020        }
     1021        else if (j == 2) { // Region D
     1022          region = psRegionSet(0,edgeWidth,
     1023                               0,ro->image->numRows);
     1024        }
     1025        else if (j == 3) { // Region C
     1026          region = psRegionSet(ro->image->numCols - edgeWidth,ro->image->numCols,
     1027                               0,ro->image->numRows);
     1028        }
     1029        psImage *subset  = psImageSubset(ro->image,region);
     1030
     1031        if (!psImageBackground(bgStats, NULL, subset, NULL, maskVal, rng)) {
     1032          psWarning("Unable to measure background for cell %d on edge %d\n", i, j);
     1033          psErrorClear();
     1034          meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PM_PATTERN_ERROR;
     1035          if (j == 0)      { B->data.F64[y][x] = NAN; }
     1036          else if (j == 1) { A->data.F64[y][x] = NAN; }
     1037          else if (j == 2) { C->data.F64[y][x] = NAN; }
     1038          else if (j == 3) { D->data.F64[y][x] = NAN; }
     1039          psFree(subset);
     1040          continue; // Move on to next edge, as only part of this cell may be a problem
     1041        }
     1042       
     1043        // If the returned value is zero, assume something is wrong.  Do I still need this?
     1044        if (psStatsGetValue(bgStats,stat) < 1e-6) {
     1045          if (j == 0)      { B->data.F64[y][x] = NAN; }
     1046          else if (j == 1) { A->data.F64[y][x] = NAN; }
     1047          else if (j == 2) { C->data.F64[y][x] = NAN; }
     1048          else if (j == 3) { D->data.F64[y][x] = NAN; }
     1049        }
     1050        // If we have an error for this cell/edge, make sure we mask the value
     1051        if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_ERROR) {
     1052          if (j == 0)      { B->data.F64[y][x] = NAN; }
     1053          else if (j == 1) { A->data.F64[y][x] = NAN; }
     1054          else if (j == 2) { C->data.F64[y][x] = NAN; }
     1055          else if (j == 3) { D->data.F64[y][x] = NAN; }
     1056        }
     1057        else { // Set the value to match what we got from the edge box.
     1058          if (xParity == -1) {
     1059            if (j == 2)      { C->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
     1060            else if (j == 3) { D->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
     1061          }
     1062          if (xParity == 1) {
     1063            if (j == 3)      { C->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
     1064            else if (j == 2) { D->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
     1065          }
     1066          if (yParity == 1) {
     1067            if (j == 0)      { B->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
     1068            else if (j == 1) { A->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
     1069          }
     1070          if (yParity == -1) {
     1071            if (j == 1)      { B->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
     1072            else if (j == 0) { A->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
     1073          }
     1074         
     1075        }
     1076#if (0)
     1077        for (int u = 0; u < subset->numCols; u++) {
     1078          for (int v = 0; v < subset->numRows; v++) {
     1079            psTrace("psModules.detrend.cont",10,"BOX: %d %d (%d %d) (%d %d) %f %d",
     1080                    i,j,x,y,u,v,subset->data.F32[v][u],0);
     1081          }
     1082        }
     1083#endif 
     1084        psFree(subset);
     1085       
     1086      }
     1087      psTrace("psModules.detrend.cont",5, "OTA: %d (%d %d) (%d %d) A: %f B: %f C: %f D: %f",
     1088              i,x,y,xParity,yParity,
     1089              A->data.F64[y][x],B->data.F64[y][x],C->data.F64[y][x],D->data.F64[y][x]);         
     1090    }
     1091    psFree(bgStats);
     1092    psFree(rng);
     1093   
     1094    // We've now allocated all the edge values, so we can now minimize the offsets.
     1095    // This involves solving the equation A x = b, where
     1096    // A is the (64x64 for GPC1) matrix containing the edges that match for each cell
     1097    // x is the solution vector
     1098    // b is the combination of offsets across each cell boundary for each cell.
     1099    // Below "XX" is used as the matrix A, and "solution" is used as both b and x
     1100    //   (due to the way psMatrixLUSolve operates).
     1101    psVector *solution = psVectorAlloc(64,PS_TYPE_F64);
     1102    psImage  *XX       = psImageAlloc(64,64,PS_TYPE_F64);
     1103    psVectorInit(solution,0.0);
     1104    psImageInit(XX,0.0);
     1105
     1106
     1107    for (int i = 0; i < numChips; i++) {
     1108      // Accumulate all the possible edge differences we can for this cell.
     1109      // As we do so, make a note of the correlations by incrementing the element of the matrix.
     1110      pmChip *chip = out->fpa->chips->data[i];
     1111      const char *chipName = psMetadataLookupStr(NULL,chip->concepts, "CHIP.NAME");
     1112      int x = chipName[2] - '0';
     1113      int y = chipName[3] - '0';
     1114
     1115      int j;
     1116      int k = 8 * x + y;
     1117      double critical_value = 0.0;
     1118
     1119      if (x + 1 < 8) {  // We have a neighbor adjacent in the +x direction
     1120        j = 8 * (x + 1) + y; // Determine that neighbor's index
     1121       
     1122        if (fabs(C->data.F64[y][x]) > fabs(D->data.F64[y][x+1])) {
     1123          critical_value = 2.0 * fabs(D->data.F64[y][x+1]);
     1124        }
     1125        else {
     1126          critical_value = 2.0 * fabs(C->data.F64[y][x]);
     1127        }
     1128        if (critical_value < 25) { critical_value = 25; }
     1129        psTrace("psModules.detrend.cont",5,"CmD %d %d %d %d %g %g %g", // diagnostic
     1130                i,x,y,j,
     1131                C->data.F64[y][x],
     1132                D->data.F64[y][x+1],
     1133                critical_value
     1134                );
     1135        if (// If there are no errors with the neighbor,
     1136            (isfinite(C->data.F64[y][x]))&&(isfinite(D->data.F64[y][x+1]))&&     // and all edges have valid values,
     1137            (fabs(C->data.F64[y][x] - D->data.F64[y][x+1]) < critical_value)     // and there are no large discontinuities,
     1138            ) {   
     1139          solution->data.F64[k] += C->data.F64[y][x] - D->data.F64[y][x+1];     // Take the difference
     1140          XX->data.F64[k][k] += 1;                                              // increment our relation with ourself
     1141          XX->data.F64[k][j] += -1;                                             // decrement our relation with the neighbor
     1142        }
     1143      }
     1144      if (x - 1 > -1) { // etc.
     1145        j = 8 * (x - 1) + y;
     1146        if (fabs(C->data.F64[y][x-1]) > fabs(D->data.F64[y][x])) {
     1147          critical_value = 2.0 * fabs(D->data.F64[y][x]);
     1148        }
     1149        else {
     1150          critical_value = 2.0 * fabs(C->data.F64[y][x-1]);
     1151        }
     1152        if (critical_value < 25) { critical_value = 25; }
     1153        psTrace("psModules.detrend.cont",5,"DmC %d %d %d %d %g %g %g",
     1154                i,x,y,j,
     1155                D->data.F64[y][x],
     1156                C->data.F64[y][x-1],
     1157                critical_value
     1158                );
     1159
     1160        if (
     1161            (isfinite(D->data.F64[y][x]))&&(isfinite(C->data.F64[y][x-1]))&&
     1162            (fabs(D->data.F64[y][x] - C->data.F64[y][x-1]) < critical_value)
     1163            ) {
     1164          solution->data.F64[k] += D->data.F64[y][x] - C->data.F64[y][x-1];
     1165          XX->data.F64[k][k] += 1;
     1166          XX->data.F64[k][j] += -1;
     1167        }
     1168      }
     1169      if (y + 1 < 8) {
     1170        j = 8 * x + (y + 1);
     1171        psTrace("psModules.detrend.cont",5,"AmB %d %d %d %d %g %g",
     1172                i,x,y,j,
     1173                A->data.F64[y][x],
     1174                B->data.F64[y+1][x]
     1175                );
     1176        if (fabs(A->data.F64[y][x]) > fabs(B->data.F64[y+1][x])) {
     1177          critical_value = 2.0 * fabs(B->data.F64[y+1][x]);
     1178        }
     1179        else {
     1180          critical_value = 2.0 * fabs(A->data.F64[y][x]);
     1181        }
     1182        if (critical_value < 25) { critical_value = 25; }
     1183        if (
     1184            (isfinite(A->data.F64[y][x]))&&(isfinite(B->data.F64[y+1][x]))&&
     1185            (fabs(A->data.F64[y][x] - B->data.F64[y+1][x]) < critical_value)
     1186            ) {
     1187          solution->data.F64[k] += A->data.F64[y][x] - B->data.F64[y+1][x];
     1188          XX->data.F64[k][k] += 1;
     1189          XX->data.F64[k][j] += -1;
     1190        }
     1191      }
     1192      if (y - 1 > -1) {
     1193        j = 8 * x +  (y - 1);
     1194        psTrace("psModules.detrend.cont",5,"BmA %d %d %d %d %g %g",
     1195                i,x,y,j,
     1196                B->data.F64[y][x],
     1197                A->data.F64[y-1][x]
     1198                );
     1199        if (fabs(A->data.F64[y-1][x]) > fabs(B->data.F64[y][x])) {
     1200          critical_value = 2.0 * fabs(B->data.F64[y][x]);
     1201        }
     1202        else {
     1203          critical_value = 2.0 * fabs(A->data.F64[y-1][x]);
     1204        }
     1205        if (critical_value < 25) { critical_value = 25; }
     1206        if (
     1207            (isfinite(B->data.F64[y][x]))&&(isfinite(A->data.F64[y-1][x]))&&
     1208            (fabs(B->data.F64[y][x] - A->data.F64[y-1][x]) < critical_value)
     1209            ) {
     1210          solution->data.F64[k] += B->data.F64[y][x] - A->data.F64[y-1][x];
     1211          XX->data.F64[k][k] += 1;
     1212          XX->data.F64[k][j] += -1;
     1213        }
     1214      }
     1215    }
     1216    double max_XX = 0;
     1217    double solution_V = 0;
     1218    int i_peak = -1;
     1219    for (int i = 0; i < numChips + 4; i++) { // If any cells have no value of themself, set the matrix to 1.0.
     1220      if (XX->data.F64[i][i] == 0.0) {
     1221        XX->data.F64[i][i] = 1.0;
     1222      }
     1223      if (XX->data.F64[i][i] > max_XX) {
     1224        max_XX = XX->data.F64[i][i];
     1225        solution_V = solution->data.F64[i];
     1226        i_peak = i;
     1227      }
     1228    }
     1229    psTrace("psModules.detrend.cont",5,"fixed point: %d %g\n",
     1230            i_peak,solution_V);
     1231   
     1232#if (1)
     1233    for (int i = 0; i < numChips + 4; i++) { // print matrix A
     1234      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",
     1235              i,
     1236              XX->data.F64[i][0],             XX->data.F64[i][1],             XX->data.F64[i][2],             XX->data.F64[i][3],
     1237              XX->data.F64[i][4],             XX->data.F64[i][5],             XX->data.F64[i][6],             XX->data.F64[i][7],
     1238              XX->data.F64[i][8],             XX->data.F64[i][9],             XX->data.F64[i][10],            XX->data.F64[i][11],
     1239              XX->data.F64[i][12],            XX->data.F64[i][13],            XX->data.F64[i][14],            XX->data.F64[i][15],
     1240              XX->data.F64[i][16],            XX->data.F64[i][17],            XX->data.F64[i][18],            XX->data.F64[i][19],
     1241              XX->data.F64[i][20],            XX->data.F64[i][21],            XX->data.F64[i][22],            XX->data.F64[i][23],
     1242              XX->data.F64[i][24],            XX->data.F64[i][25],            XX->data.F64[i][26],            XX->data.F64[i][27],
     1243              XX->data.F64[i][28],            XX->data.F64[i][29],            XX->data.F64[i][30],            XX->data.F64[i][31],
     1244              XX->data.F64[i][32],            XX->data.F64[i][33],            XX->data.F64[i][34],            XX->data.F64[i][35],
     1245              XX->data.F64[i][36],            XX->data.F64[i][37],            XX->data.F64[i][38],            XX->data.F64[i][39],
     1246              XX->data.F64[i][40],            XX->data.F64[i][41],            XX->data.F64[i][42],            XX->data.F64[i][43],
     1247              XX->data.F64[i][44],            XX->data.F64[i][45],            XX->data.F64[i][46],            XX->data.F64[i][47],
     1248              XX->data.F64[i][48],            XX->data.F64[i][49],            XX->data.F64[i][50],            XX->data.F64[i][51],
     1249              XX->data.F64[i][52],            XX->data.F64[i][53],            XX->data.F64[i][54],            XX->data.F64[i][55],
     1250              XX->data.F64[i][56],            XX->data.F64[i][57],            XX->data.F64[i][58],            XX->data.F64[i][59],
     1251              XX->data.F64[i][60],            XX->data.F64[i][61],            XX->data.F64[i][62],            XX->data.F64[i][63]
     1252              );
     1253    }
     1254
     1255    for (int i = 0; i < numChips + 4; i++) { // print vector b
     1256      psTrace("psModules.detrend.cont",5,"b: %d %f",
     1257              i,
     1258              solution->data.F64[i]
     1259              );
     1260    }
     1261#endif   
     1262   
     1263    // Solve the Ax=b equation
     1264    //    psMatrixLUSolve(XX,solution);
     1265    psMatrixGJSolve(XX,solution);
     1266#if (1)
     1267    for (int i = 0; i < numChips + 4; i++) { // print vector b
     1268      psTrace("psModules.detrend.cont",5,"x: %d %f",
     1269              i,
     1270              solution->data.F64[i]
     1271              );
     1272    }
     1273#endif
     1274   
     1275    // Cleanup
     1276    psFree(XX);
     1277    psFree(A);
     1278    psFree(B);
     1279    psFree(C);
     1280    psFree(D);
     1281
     1282    // Correct cells based on the offsets calculated, and store the result in the analysis metadata.
     1283    for (int i = 0; i < numChips; i++) {
     1284        pmChip *chip = out->fpa->chips->data[i];
     1285        pmCell *cell = chip->cells->data[0]; // Cell of interest
     1286        pmReadout *ro = cell->readouts->data[0]; // Readout of interest
     1287        const char *chipName = psMetadataLookupStr(NULL,chip->concepts, "CHIP.NAME");
     1288        int x = chipName[2] - '0';
     1289        int y = chipName[3] - '0';
     1290        int k = 8 * x + y;
     1291        float correction = solution->data.F64[k];
     1292        psLogMsg("psModules.detrend", PS_LOG_DETAIL, "Correcting background of chip %s by %f",
     1293                 chipName, correction);
     1294        psBinaryOp(ro->image, ro->image, "-", psScalarAlloc(correction, PS_TYPE_F32));
     1295        psMetadataAddF32(ro->analysis, PS_LIST_TAIL, PM_PATTERN_CELL_CORRECTION, PS_META_REPLACE,
     1296                         "Pattern chip correction solution", correction);
     1297    }
     1298
     1299    psFree(solution);
     1300    psFree(meanMask);
     1301
     1302    return true;
     1303}
Note: See TracChangeset for help on using the changeset viewer.