Changeset 4188 for trunk/psphot/src/find_defects.c

Timestamp:

Jun 9, 2005, 6:28:24 AM (21 years ago)

Author:

eugene

Message:

exploring the defect strategy

File:

• trunk/psphot/src/find_defects.c (modified) (6 diffs)

trunk/psphot/src/find_defects.c

@@ -10 +10 @@
   float gradx, grady, gradlim;
   float N0, N1, peak, sky;
-  int x, y;
+  int x, y, Nzap;
+
+  psTimerStart ("psphot");
 
   // measure the x and y max gradients
@@ -56 +58 @@
       
   // search for pixels with (dfdy > grady) or (dfdx > gradx) allowing for Nsigma outliers
+  Nzap = 0;
   for (int i = 0; i < sources->n; i++) {
     psSource *source = sources->data[i];
@@ -66 +69 @@
     sky  = source->moments->Sky;
     peak = source->moments->Peak - sky;
+    if (peak < 0) continue;
     if (peak < peak*PS_MIN(gradx, grady) + NSIGMA*sqrt(2*sky + peak)) continue;
     psTrace (".psphot.defects", 4, "peak: %f, SN: %f\n", peak, source->moments->SN);
@@ -74 +78 @@
         N0 = image->data.F32[iy][ix];
         N1 = image->data.F32[iy][ix + 1];
-        gradlim = fabs(N1 - N0 - NSIGMA*(sqrt(N1) + sqrt(N0))) / peak;
-        // XXX this line is wrong: we should use the appropriate noise from the noise map
+        gradlim = fabs(N1 - N0) - NSIGMA*(sqrt(N1 + N0));
+        // gradient - N*sigma
+        gradlim = gradlim / peak;
         
         // if gradient is too large, zap the high pixel
         if (gradlim < gradx) continue;
-        if (N1 > N0) {
+        if (N0 > N1) {
+          zap->data.U8[iy][ix] = 1;
+        } else {
           zap->data.U8[iy][ix + 1] = 1;
-        } else {
-          zap->data.U8[iy][ix] = 1;
         }
       }
@@ -92 +97 @@
         N0 = image->data.F32[iy][ix];
         N1 = image->data.F32[iy + 1][ix];
-        gradlim = fabs(N1 - N0 - NSIGMA*(sqrt(N1) + sqrt(N0))) / peak;
-        // XXX this line is wrong: we should use the appropriate noise from the noise map
+        gradlim = fabs(N1 - N0) - NSIGMA*(sqrt(N1 + N0));
+        gradlim = gradlim / peak;
         
         // if gradient is too large, zap the high pixel
         // note the inefficiency: we are re-testing known zapped entries
-        if (gradlim > grady) continue;
-        if (N1 > N0) {
+        if (gradlim < grady) continue;
+        if (N0 > N1) {
+          zap->data.U8[iy][ix] = 1;
+        } else {
           zap->data.U8[iy + 1][ix] = 1;
-        } else {
-          zap->data.U8[iy][ix] = 1;
         }           
       }
     }
-    psFree (zap);
 
     // mask & zero the zapped pixels
@@ -114 +118 @@
         y = iy + image->row0;
         zapmask->data.F32[y][x] = 1;
+        Nzap ++;
       }
     }
+    psFree (zap);
 
   }
+  psLogMsg ("psphot", 3, "zap: %f sec, %d pixels\n", psTimerMark ("psphot"), Nzap);
+
   return (true);
 }