[libcvd-members] libcvd cvd/vision.h cvd_src/vision.cc

[Ethan Eade]

CVSROOT:        /cvsroot/libcvd
Module name:    libcvd
Changes by:     Ethan Eade <ethaneade>  07/09/28 13:27:05

Modified files:
        cvd            : vision.h 
        cvd_src        : vision.cc 

Log message:
        Added 3x3 median filtering to vision.h, as:
        
        template <class T> void median_filter_3x3(const SubImage<T>& in, 
SubImage<T> out);
        
        The one-pixel border is left unspecified.  A specialized version for 
byte
        images is included, using SSE2.  On my workstation, this filters a 
640x480
        frame in 2.8 ms.

CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/libcvd/cvd/vision.h?cvsroot=libcvd&r1=1.25&r2=1.26
http://cvs.savannah.gnu.org/viewcvs/libcvd/cvd_src/vision.cc?cvsroot=libcvd&r1=1.7&r2=1.8

Patches:
Index: cvd/vision.h
===================================================================
RCS file: /cvsroot/libcvd/libcvd/cvd/vision.h,v
retrieving revision 1.25
retrieving revision 1.26
diff -u -b -r1.25 -r1.26
--- cvd/vision.h        15 Aug 2007 02:05:33 -0000      1.25
+++ cvd/vision.h        28 Sep 2007 13:27:04 -0000      1.26
@@ -391,57 +391,73 @@
 template <class T> 
 int transform(const BasicImage<T>& in, BasicImage<T>& out, const 
TooN::Matrix<2>& M, const TooN::Vector<2>& inOrig, const TooN::Vector<2>& 
outOrig)
 {
-   int i,j;
    const int w = out.size().x, h = out.size().y, iw = in.size().x, ih = 
in.size().y; 
-   TooN::Vector<2> base; 
-   TooN::Vector<2> p;
-
-   TooN::Vector<2> across = M.T()[0];
-   TooN::Vector<2> down =   M.T()[1];
-   
-   //min and max x and y
-   base[0] = std::min( across[0]*(-outOrig[0]), 
across[0]*(-outOrig[0]+(double)w) ) +           
-             std::min( down[0]*(-outOrig[1]), down[0]*(-outOrig[1]+(double)h) 
) + inOrig[0];
-   base[1] = std::min( across[1]*(-outOrig[0]), 
across[1]*(-outOrig[0]+(double)w) ) +           
-             std::min( down[1]*(-outOrig[1]), down[1]*(-outOrig[1]+(double)h) 
) + inOrig[1];
-   p[0]    = std::max( across[0]*(-outOrig[0]), 
across[0]*(-outOrig[0]+(double)w) ) +           
-             std::max( down[0]*(-outOrig[1]), down[0]*(-outOrig[1]+(double)h) 
) + inOrig[0];
-   p[1]    = std::max( across[1]*(-outOrig[0]), 
across[1]*(-outOrig[0]+(double)w) ) +           
-             std::max( down[1]*(-outOrig[1]), down[1]*(-outOrig[1]+(double)h) 
) + inOrig[1];
+    const TooN::Vector<2> across = M.T()[0];
+    const TooN::Vector<2> down =   M.T()[1];
 
+    const TooN::Vector<2> p0 = inOrig - M*outOrig;
+    const TooN::Vector<2> p1 = p0 + w*across;
+    const TooN::Vector<2> p2 = p0 + h*down;
+    const TooN::Vector<2> p3 = p0 + w*across + h*down;
+        
+    // ul --> p0
+    // ur --> w*across + p0
+    // ll --> h*down + p0
+    // lr --> w*across + h*down + p0
+    double min_x = p0[0], min_y = p0[1];
+    double max_x = min_x, max_y = min_y;
+   
+    // Minimal comparisons needed to determine bounds
+    if (across[0] < 0)
+       min_x += w*across[0];
+    else
+       max_x += w*across[0];
+    if (down[0] < 0)
+       min_x += h*down[0];
+    else
+       max_x += h*down[0];
+    if (across[1] < 0)
+       min_y += w*across[1];
+    else
+       max_y += w*across[1];
+    if (down[1] < 0)
+       min_y += h*down[1];
+    else
+       max_y += h*down[1];
    
+    // This gets from the end of one row to the beginning of the next
+    const TooN::Vector<2> carriage_return = down - w*across;
 
    //If the patch being extracted is completely in the image then no 
    //check is needed with each point.
-   if ( p[0] < iw-1 && p[1] < ih-1 && base[0] >=0 && base[1] >=0 )
-   {
-      base = M * -outOrig  + inOrig;
-      for (j=0;j<h;++j,base+=down) 
+    if (min_x >= 0 && min_y >= 0 && max_x < iw-1 && max_y < ih-1) 
       {
-         p = base;    
-         for (i=0;i<w;++i,p+=across) 
-           sample(in,p[0],p[1],out[j][i]);
-      }
-   } else {
-      int tmp = 0;
-         base = M * -outOrig  + inOrig;
-      for (j=0;j<h;++j,base+=down) 
-      {
-         p = base;    
-         for (i=0;i<w;++i, p+=across) 
+       TooN::Vector<2> p = p0;
+       for (int i=0; i<h; ++i, p+=carriage_return)
+           for (int j=0; j<w; ++j, p+=across) 
+               sample(in,p[0],p[1],out[i][j]);
+       return 0;
+    } 
+    else // Check each source location
                 {
+       // Store as doubles to avoid conversion cost for comparison
+       const double x_bound = iw-1;
+       const double y_bound = ih-1;
+       int count = 0;
+       TooN::Vector<2> p = p0;
+       for (int i=0; i<h; ++i, p+=carriage_return) {
+           for (int j=0; j<w; ++j, p+=across) {
             //Make sure that we are extracting pixels in the image
-            if ( p[0] < 0 ||  p[1] < 0 || p[0] >= iw-1 ||  p[1] >= ih-1)
-            {
-               zeroPixel(out[j][i]);
-                          tmp++;
-            } else
-               sample(in,p[0],p[1],out[j][i]);
+               if (0 <= p[0] && 0 <= p[1] &&  p[0] < x_bound && p[1] < y_bound)
+                   sample(in,p[0],p[1],out[i][j]);
+               else {
+                   zeroPixel(out[i][j]);
+                   ++count;
          }
       }
-         return tmp;
    }
-   return 0;
+       return count;
+    }
 }
 
   template <class T>  void transform(const BasicImage<T>& in, BasicImage<T>& 
out, const TooN::Matrix<3>& Minv /* <-- takes points in "out" to points in "in" 
*/)
@@ -490,6 +506,63 @@
 }
 
 
+namespace median {
+    template <class T> inline T median3(T a, T b, T c) {
+       if (b<a)
+           return std::max(b,std::min(a,c));
+       else
+           return std::max(a,std::min(b,c));   
+    }
+    
+    template <class T> inline void sort3(T& a, T& b, T& c) {
+       using std::swap;
+       if (b<a) swap(a,b);
+       if (c<b) swap(b,c);
+       if (b<a) swap(a,b);
+    }
+    
+    template <class T> T median_3x3(const T* p, const int w) {
+       T a = p[-w-1], b = p[-w], c = p[-w+1], d=p[-1], e=p[0], f=p[1], 
g=p[w-1], h=p[w], i=p[w+1];
+       sort3(a,b,c);
+       sort3(d,e,f);
+       sort3(g,h,i);
+       e = median3(b,e,h);
+       g = std::max(std::max(a,d),g);
+       c = std::min(c,std::min(f,i));
+       return median3(c,e,g);
+    }
+    
+    template <class T> void median_filter_3x3(const T* p, const int w, const 
int n, T* out)
+    {
+       T a = p[-w-1], b = p[-w], d=p[-1], e=p[0], g=p[w-1], h=p[w];
+       sort3(a,d,g);
+       sort3(b,e,h);
+       for (int j=0; j<n; ++j, ++p, ++out) {
+           T c = p[-w+1], f = p[1], i = p[w+1];
+           sort3(c,f,i);
+           *out = median3(std::min(std::min(g,h),i), 
+                          median3(d,e,f), 
+                          std::max(std::max(a,b),c));
+           a=b; b=c; d=e; e=f; g=h; h=i;
+       }
+    }
+}
+
+    template <class T> void median_filter_3x3(const SubImage<T>& I, 
SubImage<T> out)
+    {
+       assert(out.size() == I.size());
+       const int s = I.row_stride();
+       const int n = I.size().x - 2;
+       for (int i=1; i+1<I.size().y; ++i)
+           median::median_filter_3x3(I[i]+1, s, n, out[i]+1);
+    }
+
+#if (CVD_HAVE_SSE2 && CVD_HAVE_EMMINTRIN)
+
+void median_filter_3x3(const SubImage<byte>& I, SubImage<byte> out);
+
+#endif
+
 //template<class T>
 
 

Index: cvd_src/vision.cc
===================================================================
RCS file: /cvsroot/libcvd/libcvd/cvd_src/vision.cc,v
retrieving revision 1.7
retrieving revision 1.8
diff -u -b -r1.7 -r1.8
--- cvd_src/vision.cc   14 Apr 2006 13:28:39 -0000      1.7
+++ cvd_src/vision.cc   28 Sep 2007 13:27:05 -0000      1.8
@@ -119,4 +119,74 @@
 #endif
 
 
+#if (CVD_HAVE_EMMINTRIN && CVD_HAVE_SSE2)
+    
+#define SORT2_SIMD_UB(a,b)                     \
+    {                                          \
+       const __m128i ta = a;                   \
+       a = _mm_min_epu8(a,b);                  \
+       b = _mm_max_epu8(ta,b);                 \
+    }
+    
+#define SORT3_SIMD_UB(a,b,c)                   \
+    {                                          \
+       SORT2_SIMD_UB(a,b);                     \
+       SORT2_SIMD_UB(b,c);                     \
+       SORT2_SIMD_UB(a,b);                     \
+    }
+    
+    namespace median {
+       inline __m128i median3(__m128i a, __m128i b, __m128i c) {
+           SORT2_SIMD_UB(a,b);
+           return _mm_max_epu8(a,_mm_min_epu8(b,c));
+       }
+       
+       __m128i median_3x3_simd(const byte* p, const int w)
+       {
+           __m128i a = _mm_loadu_si128((const __m128i*)(p-w-1));
+           __m128i b = _mm_loadu_si128((const __m128i*)(p-w));
+           __m128i c = _mm_loadu_si128((const __m128i*)(p-w+1));
+           SORT3_SIMD_UB(a,b,c);
+           __m128i d = _mm_loadu_si128((const __m128i*)(p-1));
+           __m128i e = _mm_loadu_si128((const __m128i*)(p));
+           __m128i f = _mm_loadu_si128((const __m128i*)(p+1));
+           SORT3_SIMD_UB(d,e,f);
+           __m128i g = _mm_loadu_si128((const __m128i*)(p+w-1));
+           __m128i h = _mm_loadu_si128((const __m128i*)(p+w));
+           __m128i i = _mm_loadu_si128((const __m128i*)(p+w+1));
+           SORT3_SIMD_UB(g,h,i);
+           e = median3(b,e,h);
+           g = _mm_max_epu8(_mm_max_epu8(a,d), g);
+           c = _mm_min_epu8(c, _mm_min_epu8(f,i));
+           return median3(c,e,g);
+       }
+       
+       void median_filter_3x3_simd(const byte* p, const int stride, const int 
n, byte* out)
+       {
+           int j = 0;
+           for (; j+15<n; j+=16, p+=16, out += 16) {
+               _mm_storeu_si128((__m128i*)out, median_3x3_simd(p,stride));
+           }
+           if (j < n) {
+               const int left = n - j;
+               if (n >= 16 && left > 8)
+                   _mm_storeu_si128((__m128i*)(out+left-16), 
median_3x3_simd(p+left-16,stride));
+               else
+                   median_filter_3x3(p, stride, left, out);
+           }
+       }
+       
+    }
+
+    void median_filter_3x3(const SubImage<byte>& I, SubImage<byte> out)
+    {
+       assert(out.size() == I.size());
+       const int s = I.row_stride();
+       const int n = I.size().x - 2;
+       for (int i=1; i+1<I.size().y; ++i)
+           median::median_filter_3x3_simd(I[i]+1, s, n, out[i]+1);
+    }
+#endif
+    
+
 };