Changeset 17336

Timestamp:

11/21/13 17:25:07 (11 years ago)

Author:

tbretz

Message:

Added a new spike removal algorithm; changed the interface from 'all pixels' to #one pixel' to allow procesing a single pixel in debugging cases.

File:

• trunk/Mars/mcore/DrsCalib.h (modified) (3 diffs)

trunk/Mars/mcore/DrsCalib.h

@@ -569 +569 @@
     }
 
-    static void RemoveSpikes(float *vec, uint32_t roi)
+    static void RemoveSpikes(float *p, uint32_t roi)
     {
         if (roi<4)
             return;
 
-        for (size_t ch=0; ch<1440; ch++)
-        {
-            float *p = vec + ch*roi;
-
-            for (size_t i=1; i<roi-2; i++)
-            {
-                if (p[i]-p[i-1]>25 && p[i]-p[i+1]>25)
-                {
-                    p[i] = (p[i-1]+p[i+1])/2;
-                }
-
-                if (p[i]-p[i-1]>22 && fabs(p[i]-p[i+1])<4 && p[i+1]-p[i+2]>22)
-                {
-                    p[i] = (p[i-1]+p[i+2])/2;
-                    p[i+1] = p[i];
-                }
-            }
-        }
-    }
-
-    static void RemoveSpikes2(float *vec, uint32_t roi)//from Werner
+        for (size_t i=1; i<roi-2; i++)
+        {
+            if (p[i]-p[i-1]>25 && p[i]-p[i+1]>25)
+            {
+                p[i] = (p[i-1]+p[i+1])/2;
+            }
+
+            if (p[i]-p[i-1]>22 && fabs(p[i]-p[i+1])<4 && p[i+1]-p[i+2]>22)
+            {
+                p[i] = (p[i-1]+p[i+2])/2;
+                p[i+1] = p[i];
+            }
+        }
+    }
+
+    static void RemoveSpikes2(float *p, uint32_t roi)//from Werner
     {
         if (roi<4)
             return;
 
-        for (size_t ch=0; ch<1440; ch++)
-        {
-            float *p = vec + ch*roi;
-
-            std::vector<float> Ameas(p, p+roi);
-
-            std::vector<float> diff(roi);
-            for (size_t i=1; i<roi-1; i++)
-                diff[i] = (p[i-1] + p[i+1])/2 - p[i];
-
-            //std::vector<float> N1mean(roi);
-            //for (size_t i=1; i<roi-1; i++)
-            //    N1mean[i] = (p[i-1] + p[i+1])/2;
-
-            const float fract = 0.8;
-
-            for (size_t i=0; i<roi-3; i++)
-            {
-                if (diff[i]<5)
-                    continue;
-
-                if (Ameas[i+2] - (Ameas[i] + Ameas[i+3])/2 > 10)
-                {
-                    p[i+1]=   (Ameas[i+3] - Ameas[i])/3 + Ameas[i];
-                    p[i+2]= 2*(Ameas[i+3] - Ameas[i])/3 + Ameas[i];
-
-                    i += 3;
-
-                    continue;
-                }
-
-                if ( (diff[i+1]<-diff[i]*fract*2) && (diff[i+2]>10) )
-                {
-                    p[i+1]    = (Ameas[i]+Ameas[i+2])/2;
-                    diff[i+2] = (p[i+1] + Ameas[i+3])/2 - Ameas[i+2];
-
-                    i += 2;
-                }
-
-                // const float x = Ameas[i] - N1mean[i];
-                // if (x > -5.)
-                //     continue;
-
-                // if (Ameas[i+2] - (Ameas[i] + Ameas[i+3])/2. > 10.)
-                // {
-                //     p[i+1]=   (Ameas[i+3] - Ameas[i])/3 + Ameas[i];
-                //     p[i+2]= 2*(Ameas[i+3] - Ameas[i])/3 + Ameas[i];
-                //     i += 3;
-                //     continue;
-                // }
-
-                // const float xp  = Ameas[i+1] - N1mean[i+1];
-                // const float xpp = Ameas[i+2] - N1mean[i+2];
-
-                // if ( (xp > -2.*x*fract) && (xpp < -10.) )
-                // {
-                //     p[i+1] = N1mean[i+1];
-                //     N1mean[i+2] = Ameas[i+1] - Ameas[i+3]/2;
-                //
-                //     i += 2;
-                // }
-            }
+        std::vector<float> Ameas(p, p+roi);
+
+        std::vector<float> diff(roi);
+        for (size_t i=1; i<roi-1; i++)
+            diff[i] = (p[i-1] + p[i+1])/2 - p[i];
+
+        //std::vector<float> N1mean(roi);
+        //for (size_t i=1; i<roi-1; i++)
+        //    N1mean[i] = (p[i-1] + p[i+1])/2;
+
+        const float fract = 0.8;
+
+        for (size_t i=0; i<roi-3; i++)
+        {
+            if (diff[i]<5)
+                continue;
+
+            if (Ameas[i+2] - (Ameas[i] + Ameas[i+3])/2 > 10)
+            {
+                p[i+1]=   (Ameas[i+3] - Ameas[i])/3 + Ameas[i];
+                p[i+2]= 2*(Ameas[i+3] - Ameas[i])/3 + Ameas[i];
+
+                i += 3;
+
+                continue;
+            }
+
+            if ( (diff[i+1]<-diff[i]*fract*2) && (diff[i+2]>10) )
+            {
+                p[i+1]    = (Ameas[i]+Ameas[i+2])/2;
+                diff[i+2] = (p[i+1] + Ameas[i+3])/2 - Ameas[i+2];
+
+                i += 2;
+            }
+
+            // const float x = Ameas[i] - N1mean[i];
+            // if (x > -5.)
+            //     continue;
+
+            // if (Ameas[i+2] - (Ameas[i] + Ameas[i+3])/2. > 10.)
+            // {
+            //     p[i+1]=   (Ameas[i+3] - Ameas[i])/3 + Ameas[i];
+            //     p[i+2]= 2*(Ameas[i+3] - Ameas[i])/3 + Ameas[i];
+            //     i += 3;
+            //     continue;
+            // }
+
+            // const float xp  = Ameas[i+1] - N1mean[i+1];
+            // const float xpp = Ameas[i+2] - N1mean[i+2];
+
+            // if ( (xp > -2.*x*fract) && (xpp < -10.) )
+            // {
+            //     p[i+1] = N1mean[i+1];
+            //     N1mean[i+2] = Ameas[i+1] - Ameas[i+3]/2;
+            //
+            //     i += 2;
+            // }
         }
     }
@@ -666 +656 @@
     static void RemoveSpikes3(float *vec, uint32_t roi)//from Werner
     {
+        if (roi<4)
+            return;
+
         const float SingleCandidateTHR = -10.;
         const float DoubleCandidateTHR =  -5.;
@@ -725 +718 @@
             }
         } 
+    }
+
+    static void RemoveSpikes4(float *ptr, uint32_t roi)
+    {
+        if (roi<7)
+            return;
+
+        for (uint32_t i=0; i<roi-6; i++)
+        {
+            double d10, d21, d32, d43, d54;
+
+            // ============================================
+            d43 = ptr[i+4]-ptr[i+3];
+            d54 = ptr[i+5]-ptr[i+4];
+
+            if ((d43>35 && -d54>35) || (d43<-35 && -d54<-35))
+            {
+                ptr[i+4] = (ptr[i+3]+ptr[i+5])/2;
+            }
+
+            // ============================================
+            d32 = ptr[i+3]-ptr[i+2];
+            d54 = ptr[i+5]-ptr[i+4];
+
+            if ((d32>9   && -d54>13  && d32-d54>31)/* || (d32<-13 && -d54<-13 && d32+d54<-63)*/)
+            {
+                double avg0 = (ptr[i+2]+ptr[i+5])/2;
+                double avg1 = (ptr[i+3]+ptr[i+4])/2;
+
+                ptr[i+3] = ptr[i+3] - avg1+avg0;
+                ptr[i+4] = ptr[i+4] - avg1+avg0;
+            }
+
+            // ============================================
+            d21 = ptr[i+2]-ptr[i+1];
+            d54 = ptr[i+5]-ptr[i+4];
+
+            if (d21>15 && -d54>17)
+            {
+                double avg0 = (ptr[i+1]+ptr[i+5])/2;
+                double avg1 = (ptr[i+2]+ptr[i+3]+ptr[i+4])/3;
+
+                ptr[i+2] = ptr[i+2] - avg1+avg0;
+                ptr[i+3] = ptr[i+3] - avg1+avg0;
+                ptr[i+4] = ptr[i+4] - avg1+avg0;
+            }
+
+            // ============================================
+            d10 = ptr[i+1]-ptr[i];
+            d54 = ptr[i+5]-ptr[i+4];
+
+            if (d10>18 && -d54>20)
+            {
+                double avg0 = (ptr[i]+ptr[i+5])/2;
+                double avg1 = (ptr[i+1]+ptr[i+2]+ptr[i+3]+ptr[i+4])/4;
+
+                ptr[i+1] = ptr[i+1] - avg1+avg0;
+                ptr[i+2] = ptr[i+2] - avg1+avg0;
+                ptr[i+3] = ptr[i+3] - avg1+avg0;
+                ptr[i+4] = ptr[i+4] - avg1+avg0;
+            }
+        }
     }