diff --git a/CV.cabal b/CV.cabal
--- a/CV.cabal
+++ b/CV.cabal
@@ -1,5 +1,5 @@
 Name:				 CV
-Version:             0.3.6.0
+Version:             0.3.7
 Description:         OpenCV Bindings
 License:             GPL
 License-file:        LICENSE
@@ -66,15 +66,21 @@
                        cbits/cvWrapLEO.h,
                        cbits/cvWrapCore.h,
                        cbits/wrapImgProc.h,
-                       cbits/cvIterators.h
+                       cbits/cvIterators.h,
+                       cbits/cvProjection.h,
+                       cbits/cvFiles.h
     c-sources:         cbits/cvWrapLEO.c,
                        cbits/cvWrapCore.c,
                        cbits/wrapImgProc.c,
-                       cbits/cvIterators.c
+                       cbits/cvIterators.c,
+                       cbits/cvProjection.c,
+                       cbits/cvFiles.c
     install-includes:  cbits/cvWrapLEO.h,
                        cbits/cvWrapCore.h,
                        cbits/wrapImgProc.h,
-                       cbits/cvIterators.h
+                       cbits/cvIterators.h,
+                       cbits/cvProjection.h,
+                       cbits/cvFiles.h
     if flag(opencv23)
         cpp-options: -DOpenCV23
         cc-options: -DOpenCV23
@@ -100,7 +106,7 @@
                        array >= 0.2.0.0,
                        directory >= 1 && < 2,
                        mtl >= 1.1.0,
-                       random >= 1.0.0,
+                       -- random >= 1.0.0,
                        carray >= 0.1.5,
                        QuickCheck >= 2.1,
                        containers >= 0.2,
@@ -112,6 +118,8 @@
                        lazysmallcheck >= 0.5 && < 1,
                        parallel-io    >= 0.3.2 && < 0.3.3,
                        storable-tuple >= 0.0.2 && <= 1,
+                       mwc-random >= 0.12.0.0,
+                       primitive  >= 0.4.1, 
                        filepath >= 1.3.0.0 && < 1.4
 
     Exposed-modules:   CV.Image                               
@@ -152,6 +160,8 @@
                        ,CV.Tracking                           
                        ,CV.Transforms                         
                        ,CV.Video                              
+                       ,CV.Projection
+                       ,CV.Files
                        ,Utils.DrawingClass                    
                        ,Utils.Function                        
                        ,Utils.GeometryClass                   
@@ -169,7 +179,9 @@
                        CV.Bindings.Drawing,
                        CV.Bindings.Error,
                        CV.Bindings.Features,
-                       CV.Bindings.Iterators
+                       CV.Bindings.Iterators,
+                       CV.Bindings.Projection,
+                       CV.Bindings.Files
     Extensions: CPP
 
 source-repository head
diff --git a/CV/Bindings/Core.hsc b/CV/Bindings/Core.hsc
--- a/CV/Bindings/Core.hsc
+++ b/CV/Bindings/Core.hsc
@@ -11,47 +11,6 @@
 #include "cvWrapCore.h"
 #include <opencv2/core/core_c.h>
 
--- CVAPI(IplImage*) cvCreateImage(
---   CvSize size,
---   int depth,
---   int channels
--- );
-
--- Creates IPL image (header and data)
-
--- need to use wrapper functions since passing struct by value.
-
--- #ccall wrapCreateImage , CInt -> CInt -> CInt -> CInt -> IO (Ptr IplImage)
-
--- Releases IPL image header and data
--- CVAPI(void) cvReleaseImage(
---   IplImage** image
--- );
-
--- Allocates and initializes CvMat header and allocates data
--- CVAPI(CvMat*) cvCreateMat(
---  int rows,
---  int cols,
---  int type
--- );
-
--- Releases CvMat header and deallocates matrix data
--- (reference counting is used for data)
--- CVAPI(void) cvReleaseMat(
---   CvMat** mat
--- );
-
--- CVAPI(void) cvSet(
---   CvArr* arr,
---   CvScalar value,
---   const CvArr* mask CV_DEFAULT(NULL)
--- );
-
--- Sets all or "masked" elements of input array
--- to the same value
-
--- need to use wrapper functions since passing struct by value.
-
 #ccall wrapSet , Ptr <CvArr> -> Ptr <CvScalar> -> Ptr <CvArr> -> IO ()
 #ccall wrapSetAll , Ptr <CvArr> -> CDouble -> Ptr <CvArr> -> IO ()
 #ccall wrapSet1 , Ptr <CvArr> -> CDouble -> Ptr <CvArr> -> IO ()
@@ -59,50 +18,14 @@
 #ccall wrapSet3 , Ptr <CvArr> -> CDouble -> CDouble -> CDouble -> Ptr <CvArr> -> IO ()
 #ccall wrapSet4 , Ptr <CvArr> -> CDouble -> CDouble -> CDouble -> CDouble -> Ptr <CvArr> -> IO ()
 
--- CVAPI(void) cvSetZero(
---   CvArr* arr
--- );
-
 -- Clears all the array elements (sets them to 0)
-
-#ccall cvSetZero , Ptr <CvArr> -> IO ()
-
--- CVAPI(void) cvSplit(
---   const CvArr* src,
---   CvArr* dst0,
---   CvArr* dst1,
---   CvArr* dst2,
---   CvArr* dst3
--- );
-
--- Splits a multi-channel array into the set of single-channel arrays or
--- extracts particular [color] plane
-
-#ccall cvSplit , Ptr <CvArr> -> Ptr <CvArr> -> Ptr <CvArr> -> Ptr <CvArr> -> Ptr <CvArr> -> IO ()
-
--- |CVAPI(void) cvMerge(
--- |  const CvArr* src0,
--- |  const CvArr* src1,
--- |  const CvArr* src2,
--- |  const CvArr* src3,
--- |  CvArr* dst
--- |);
--- |Merges a set of single-channel arrays into the single multi-channel array
--- |or inserts one particular [color] plane to the array
-#ccall cvMerge , Ptr <CvArr> -> Ptr <CvArr> -> Ptr <CvArr> -> Ptr <CvArr> -> Ptr <CvArr> -> IO ()
+#ccall cvSetZero , Ptr BareImage -> IO ()
 
--- CVAPI(void) cvCartToPolar(
---   const CvArr* x,
---   const CvArr* y,
---   CvArr* magnitude,
---   CvArr* angle CV_DEFAULT(NULL),
---   int angle_in_degrees CV_DEFAULT(0)
--- );
+#ccall cvSplit , Ptr BareImage -> Ptr BareImage -> Ptr BareImage -> Ptr BareImage -> Ptr BareImage -> IO ()
 
--- Does cartesian->polar coordinates conversion.
--- Either of output components (magnitude or angle) is optional
+#ccall cvMerge , Ptr BareImage -> Ptr BareImage -> Ptr BareImage -> Ptr BareImage -> Ptr BareImage -> IO ()
 
-#ccall cvCartToPolar , Ptr <CvArr> -> Ptr <CvArr> -> Ptr <CvArr> -> Ptr <CvArr> -> CInt -> IO ()
+#ccall cvCartToPolar , Ptr BareImage -> Ptr BareImage -> Ptr BareImage -> Ptr BareImage -> CInt -> IO ()
 
 -- | CVAPI(void) cvPolarToCart(
 -- |   const CvArr* magnitude,
@@ -165,7 +88,7 @@
 --   const CvArr* mask CV_DEFAULT(NULL)
 -- );
 
-#ccall cvNormalize , Ptr <CvArr> -> Ptr <CvArr> -> CDouble -> CDouble -> CInt -> Ptr <CvArr> -> IO ()
+#ccall cvNormalize , Ptr BareImage -> Ptr BareImage -> CDouble -> CDouble -> CInt -> Ptr BareImage -> IO ()
 
 -- stuff related to DFT
 
@@ -178,6 +101,7 @@
 #num CV_DXT_COMPLEX_OUTPUT
 #num CV_DXT_REAL_OUTPUT
 #num CV_DXT_INV_REAL
+#num CV_DXT_INVERSE_SCALE
 
 -- CVAPI(void) cvDFT(
 --   const CvArr* src,
@@ -191,11 +115,7 @@
 --   real->ccs (forward),
 --   ccs->real (inverse)
 
-#ccall cvDFT , Ptr <CvArr> -> Ptr <CvArr> -> CInt -> CInt -> IO ()
-
--- #ccall cvWrapDFT , Ptr <CvArr> -> Ptr <CvArr> -> IO ()
-
--- #ccall cvWrapIDFT , Ptr <CvArr> -> Ptr <CvArr> -> IO ()
+#ccall cvDFT , Ptr BareImage -> Ptr BareImage -> CInt -> CInt -> IO ()
 
 -- CVAPI(void) cvMulSpectrums(
 --   const CvArr* src1,
@@ -204,11 +124,9 @@
 --   int flags
 -- );
 
-#ccall swapQuadrants, Ptr <CvArr> -> IO ()
-
 -- Multiply results of DFTs: DFT(X)*DFT(Y) or DFT(X)*conj(DFT(Y))
 
-#ccall cvMulSpectrums , Ptr <CvArr> -> Ptr <CvArr> -> Ptr <CvArr> -> CInt -> IO ()
+#ccall cvMulSpectrums , Ptr BareImage -> Ptr BareImage -> Ptr BareImage -> CInt -> IO ()
 
 -- CVAPI(int) cvGetOptimalDFTSize(
 --   int size0
@@ -227,5 +145,7 @@
 
 -- Discrete Cosine Transform
 
-#ccall cvDCT , Ptr <CvArr> -> Ptr <CvArr> -> CInt -> IO ()
+#ccall cvDCT , Ptr BareImage -> Ptr BareImage -> CInt -> IO ()
 
+-- TODO: This might be expanded:
+#opaque_t CvSeq
diff --git a/CV/Bindings/Files.hsc b/CV/Bindings/Files.hsc
new file mode 100644
--- /dev/null
+++ b/CV/Bindings/Files.hsc
@@ -0,0 +1,20 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+module CV.Bindings.Files where
+
+import Foreign.Ptr
+import Foreign.C.String
+import CV.Bindings.Types
+import CV.Image(BareImage)
+
+#strict_import
+
+#include <bindings.dsl.h>
+#include "cvFiles.h"
+
+#ccall read_from_tcr , CString ->  Ptr (Ptr CUInt) -> IO (Ptr BareImage)
+
+#ccall read_from_tcr_rectified , CString -> IO (Ptr BareImage)
+
+#ccall read_from_tcr_mag , CString -> IO (Ptr BareImage)
+
+#ccall calc_derivative_direction , Ptr BareImage -> IO (Ptr BareImage)
diff --git a/CV/Bindings/ImgProc.hsc b/CV/Bindings/ImgProc.hsc
--- a/CV/Bindings/ImgProc.hsc
+++ b/CV/Bindings/ImgProc.hsc
@@ -52,7 +52,6 @@
             (fromIntegral r)
             (cBorderType border)
             (realToFrac value)
-
 -- moments
 
 -- | Calculates all spatial and central moments up to the 3rd order
diff --git a/CV/Bindings/Matrix.hsc b/CV/Bindings/Matrix.hsc
--- a/CV/Bindings/Matrix.hsc
+++ b/CV/Bindings/Matrix.hsc
@@ -26,9 +26,15 @@
 #ccall cvReleaseMat , Ptr (Ptr <CvMat>) -> IO ()
 
 #ccall cvTranspose  , Ptr <CvMat> -> Ptr <CvMat> -> IO ()
+#ccall cvInvert     , Ptr <CvMat> -> Ptr <CvMat> -> CInt -> IO ()
 #ccall cvGEMM       , Ptr <CvMat> -> Ptr <CvMat> -> Double -> Ptr <CvMat> -> Double -> Ptr <CvMat> -> Int -> IO ()
 
 #ccall cvRodrigues2  , Ptr <CvMat> -> Ptr <CvMat> -> Ptr <CvMat> -> IO Int
+
+-- Matrix inversions
+#num CV_LU 
+#num CV_SVD 
+#num CV_SVD_SYM 
 
 #num CV_GEMM_A_T
 #num CV_GEMM_B_T
diff --git a/CV/Bindings/Projection.hsc b/CV/Bindings/Projection.hsc
new file mode 100644
--- /dev/null
+++ b/CV/Bindings/Projection.hsc
@@ -0,0 +1,13 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+module CV.Bindings.Projection where
+
+import Foreign.Ptr
+import CV.Bindings.Types
+import CV.Image(BareImage)
+
+#strict_import
+
+#include <bindings.dsl.h>
+#include "cvProjection.h"
+
+#ccall project_polar , Ptr BareImage -> IO (Ptr BareImage)
diff --git a/CV/Bindings/Types.hsc b/CV/Bindings/Types.hsc
--- a/CV/Bindings/Types.hsc
+++ b/CV/Bindings/Types.hsc
@@ -151,11 +151,19 @@
 #field height , CInt
 #stoptype
 
+instance Sized C'CvSize where
+    type Size  C'CvSize = (Int,Int)
+    getSize (C'CvSize w h)  = (fromIntegral w,fromIntegral h)
+
 #starttype CvSize2D32f
 #field width , CFloat
 #field height , CFloat
 #stoptype
 
+instance Sized C'CvSize2D32f where
+    type Size  C'CvSize2D32f = (Float,Float)
+    getSize (C'CvSize2D32f w h)  = (realToFrac w, realToFrac h)
+
 #starttype CvConnectedComp
 #field area, CDouble
 #field value, <CvScalar>
@@ -183,6 +191,10 @@
    pt (C'CvPoint2D32f x y) = (realToFrac x,realToFrac y)
    toPt (x,y) = C'CvPoint2D32f (realToFrac x) (realToFrac y)
 
+instance Sized C'CvBox2D where
+    type Size  C'CvBox2D = (Float,Float)
+    getSize (C'CvBox2D _ s _) = getSize s
+
 -- // #starttype CV_32FC2
 -- // #field x , Float
 -- // #field y , Float
@@ -255,6 +267,7 @@
 #num CV_TERMCRIT_ITER    
 #num CV_TERMCRIT_NUMBER  
 #num CV_TERMCRIT_EPS     
+
 
 
 -- Memory Storage
diff --git a/CV/Conversions.hs b/CV/Conversions.hs
--- a/CV/Conversions.hs
+++ b/CV/Conversions.hs
@@ -22,7 +22,7 @@
     ,acquireImageSlow8URGB'
     ) where
 
-import Data.Complex
+import Data.Complex as C
 
 import CV.Image
 import Data.Word
@@ -76,7 +76,7 @@
 
 
 -- |Copy the real part of an array to image
-copyComplexCArrayToImage :: CArray (Int,Int) (Complex Double) -> Image GrayScale D32
+copyComplexCArrayToImage :: CArray (Int,Int) (C.Complex Double) -> Image GrayScale D32
 copyComplexCArrayToImage carr = S $ unsafePerformIO $
                           creatingBareImage (withCArray carr (acquireImageSlowComplex' w h))
     where
@@ -99,7 +99,7 @@
      (w,h) = getSize img
 
 -- |Copy image as a real part of a complex CArray
-copyImageToComplexCArray :: Image GrayScale D32 -> CArray (Int,Int) (Complex Double)
+copyImageToComplexCArray :: Image GrayScale D32 -> CArray (Int,Int) (C.Complex Double)
 copyImageToComplexCArray (S img) = unsafePerformIO $
          withBareImage img $ \cimg ->
           createCArray ((0,0),(w-1,h-1)) (exportImageSlowComplex' cimg) --({#call exportImageSlow#} cimg)
@@ -113,7 +113,7 @@
   exportImageSlowF' :: ((Ptr (BareImage)) -> ((Ptr Float) -> (IO ())))
 
 foreign import ccall safe "CV/cvWrapLeo.h exportImageSlowComplex"
-  exportImageSlowComplex' :: ((Ptr (BareImage)) -> ((Ptr (Complex Double)) -> (IO ())))
+  exportImageSlowComplex' :: ((Ptr (BareImage)) -> ((Ptr (C.Complex Double)) -> (IO ())))
 
 foreign import ccall safe "CV/cvWrapLeo.h acquireImageSlow"
   acquireImageSlow' :: (Int -> (Int -> ((Ptr Double) -> (IO (Ptr (BareImage))))))
@@ -131,5 +131,5 @@
   acquireImageSlow8U' :: (Int -> (Int -> ((Ptr Word8) -> (IO (Ptr (BareImage))))))
 
 foreign import ccall safe "CV/cvWrapLeo.h acquireImageSlowComplex"
-  acquireImageSlowComplex' :: (Int -> (Int -> ((Ptr (Complex Double)) -> (IO (Ptr (BareImage))))))
+  acquireImageSlowComplex' :: (Int -> (Int -> ((Ptr (C.Complex Double)) -> (IO (Ptr (BareImage))))))
 
diff --git a/CV/DFT.hs b/CV/DFT.hs
--- a/CV/DFT.hs
+++ b/CV/DFT.hs
@@ -1,143 +1,61 @@
 {-# LANGUAGE ScopedTypeVariables #-}
 module CV.DFT where
 
+import CV.Bindings.Types
 import CV.Bindings.Core
 import CV.Bindings.ImgProc
-import CV.Bindings.Types
 import CV.Image
-import CV.Operations
-
-import Foreign.Ptr (castPtr,nullPtr)
+import CV.ImageMath as IM
+import CV.ImageMathOp
 import System.IO.Unsafe
-import Data.Complex
-
-type I32 = Image GrayScale D32
-type Idft32 = Image DFT D32
-data Ipolar32 = Ipolar32 (Complex D32)
-
-data Swap = SwapQuadrants | DontSwapQuadrants deriving (Eq,Show)
+import Foreign.Ptr
 
-dft = dft' SwapQuadrants
-dft' :: Swap -> Image GrayScale d -> Image DFT D32
-dft' swap i = unsafePerformIO $ do
-  --n::(Image GrayScale D32) <- create (w', h')
-  --n <- copyMakeBorder i 0 (h'-h) 0 (w'-w) BorderReplicate 0
-  z::(Image GrayScale D32) <- create (w, h)
-  d::(Image DFT D32) <- create (w, h)
-  withImage i $ \i_ptr ->
-    withImage z $ \z_ptr ->
-      withImage d $ \d_ptr -> do
-        c'cvMerge (castPtr i_ptr) (castPtr z_ptr) nullPtr nullPtr (castPtr d_ptr)
-        c'cvDFT (castPtr d_ptr) (castPtr d_ptr) c'CV_DXT_FORWARD (fromIntegral 0)
-        case swap of
-         SwapQuadrants -> c'swapQuadrants (castPtr d_ptr)
-         DontSwapQuadrants -> return ()
-        return d
+dft :: Image GrayScale d -> Image Complex D32
+dft i = unsafePerformIO $ do
+  n::(Image GrayScale D32) <- create (w', h')
+  n <- copyMakeBorder i 0 (h'-h) 0 (w'-w) BorderConstant 0
+  c::(Image Complex D32) <- create (w', h')
+  withGenImage n $ \nimg ->
+    withGenImage c $ \cimg -> do
+      c'cvMerge nimg nullPtr nullPtr nullPtr cimg
+      c'cvDFT cimg cimg c'CV_DXT_FORWARD 0
+      return c
   where
     (w,h) = getSize i
-    --w' = fromIntegral $ c'cvGetOptimalDFTSize (fromIntegral w)
-    --h' = fromIntegral $ c'cvGetOptimalDFTSize (fromIntegral h)
+    w' = fromIntegral $ c'cvGetOptimalDFTSize (fromIntegral w)
+    h' = fromIntegral $ c'cvGetOptimalDFTSize (fromIntegral h)
 
-idft = idft' SwapQuadrants
-idft' :: Swap -> Image DFT D32 -> Image GrayScale D32
-idft' swap d = unsafePerformIO $ do
+idft :: Image Complex D32 -> Image GrayScale D32
+idft c = unsafePerformIO $ do
   n::(Image GrayScale D32) <- create s
-  --z::(Image GrayScale D32) <- create s
-  withImage d $ \d_ptr ->
-    withImage n $ \n_ptr -> do
-      --withImage z $ \z_ptr -> do
-        case swap of
-         SwapQuadrants -> c'swapQuadrants (castPtr d_ptr)
-         DontSwapQuadrants -> return ()
-        --c'swapQuadrants (castPtr d_ptr)
-        c'cvDFT (castPtr d_ptr) (castPtr n_ptr) c'CV_DXT_INV_SCALE (fromIntegral 0)
-        --c'cvSplit (castPtr d_ptr) (castPtr n_ptr) (castPtr z_ptr) nullPtr nullPtr
-        return n
+  withGenImage c $ \c_ptr ->
+    withGenImage n $ \n_ptr -> do
+      c'cvDFT c_ptr c_ptr c'CV_DXT_INVERSE 0
+      c'cvSplit c_ptr n_ptr nullPtr nullPtr nullPtr
+      return n
   where
-    s = getSize d
-
-swapQuadrants img = unsafePerformIO $ do
-            e <- cloneImage img
-            withImage e (c'swapQuadrants . castPtr)
-            return e
+    s = getSize c
 
-dftSplit :: Image DFT D32 -> (Image GrayScale D32, Image GrayScale D32)
-dftSplit d = unsafePerformIO $ do
+complexSplit :: Image Complex D32 -> (Image GrayScale D32, Image GrayScale D32)
+complexSplit c = unsafePerformIO $ do
   re::(Image GrayScale D32) <- create (w, h)
   im::(Image GrayScale D32) <- create (w, h)
-  withImage d $ \d_ptr ->
-    withImage re $ \re_ptr ->
-      withImage im $ \im_ptr -> do
-        c'cvSplit (castPtr d_ptr) (castPtr re_ptr) (castPtr im_ptr) nullPtr nullPtr
+  withGenImage c $ \c_ptr ->
+    withGenImage re $ \re_ptr ->
+      withGenImage im $ \im_ptr -> do
+        c'cvSplit c_ptr re_ptr im_ptr nullPtr nullPtr
         return (re,im)
   where
-    (w,h) = getSize d
-
-dftMerge :: (Image GrayScale D32, Image GrayScale D32) -> Image DFT D32
-dftMerge (re,im) = unsafePerformIO $ do
-  d::(Image DFT D32) <- create (w, h)
-  withImage re $ \re_ptr ->
-    withImage im $ \im_ptr ->
-      withImage d $ \d_ptr -> do
-        c'cvMerge (castPtr re_ptr) (castPtr im_ptr) nullPtr nullPtr (castPtr d_ptr)
-  return d
-  where
-    (w,h) = getSize re
+    (w,h) = getSize c
 
-dftToPolar :: Image DFT D32 -> (Image GrayScale D32, Image GrayScale D32)
-dftToPolar d = unsafePerformIO $ do
+complexToMagnitude :: Image Complex D32 -> Image GrayScale D32
+complexToMagnitude c = unsafePerformIO $ do
   mag::(Image GrayScale D32) <- create (w, h)
-  ang::(Image GrayScale D32) <- create (w, h)
-  withImage re $ \re_ptr ->
-    withImage im $ \im_ptr ->
-      withImage mag $ \mag_ptr ->
-        withImage ang $ \ang_ptr ->
-          c'cvCartToPolar (castPtr re_ptr) (castPtr im_ptr) (castPtr mag_ptr) (castPtr ang_ptr) (fromIntegral 0)
-  return (mag,ang)
-  where
-    (re,im) = dftSplit d
-    (w,h) = getSize d
-
-dftFromPolar :: (Image GrayScale D32, Image GrayScale D32) -> Image DFT D32
-dftFromPolar (mag,ang) = unsafePerformIO $ do
-  re::(Image GrayScale D32) <- create (w, h)
-  im::(Image GrayScale D32) <- create (w, h)
-  withImage mag $ \mag_ptr ->
-    withImage ang $ \ang_ptr ->
-      withImage re $ \re_ptr -> do
-        withImage im $ \im_ptr -> do
-          c'cvPolarToCart (castPtr mag_ptr) (castPtr ang_ptr) (castPtr re_ptr) (castPtr im_ptr) (fromIntegral 0)
-  return $ dftMerge (re,im)
-  where
-    (w,h) = getSize mag
-
---withPolar :: Image DFT D32 ->
---    ((Image GrayScale D32, Image GrayScale D32) -> (Image GrayScale D32, Image GrayScale D32)) ->
---    Image DFT D32
-
-
-rgbSplit :: Image RGB D32 -> (Image GrayScale D32, Image GrayScale D32, Image GrayScale D32)
-rgbSplit i = unsafePerformIO $ do
-  r::(Image GrayScale D32) <- create (w, h)
-  g::(Image GrayScale D32) <- create (w, h)
-  b::(Image GrayScale D32) <- create (w, h)
-  withImage i $ \i_ptr ->
-    withImage r $ \r_ptr ->
-      withImage g $ \g_ptr ->
-        withImage b $ \b_ptr -> do
-          c'cvSplit (castPtr i_ptr) (castPtr r_ptr) (castPtr g_ptr) (castPtr b_ptr) nullPtr
-          return (r,g,b)
-  where
-    (w,h) = getSize i
-
-rgbMerge :: (Image GrayScale D32, Image GrayScale D32, Image GrayScale D32) -> Image RGB D32
-rgbMerge (r,g,b) = unsafePerformIO $ do
-  i::(Image RGB D32) <- create (w, h)
-  withImage r $ \r_ptr ->
-    withImage g $ \g_ptr ->
-      withImage b $ \b_ptr ->
-        withImage i $ \i_ptr -> do
-          c'cvMerge (castPtr r_ptr) (castPtr g_ptr) (castPtr b_ptr) nullPtr (castPtr i_ptr)
-          return i
+  withGenImage re $ \re_ptr ->
+    withGenImage im $ \im_ptr ->
+      withGenImage mag $ \mag_ptr -> do
+        c'cvCartToPolar re_ptr im_ptr mag_ptr nullPtr (fromIntegral 0)
+        return $ IM.log $ 1 |+ mag
   where
-    (w,h) = getSize r
+    (re,im) = complexSplit c
+    (w,h) = getSize c
diff --git a/CV/Drawing.chs b/CV/Drawing.chs
--- a/CV/Drawing.chs
+++ b/CV/Drawing.chs
@@ -1,4 +1,4 @@
-{-#LANGUAGE ForeignFunctionInterface, TypeFamilies, MultiParamTypeClasses, TypeSynonymInstances, 
+{-#LANGUAGE ForeignFunctionInterface, TypeFamilies, MultiParamTypeClasses, TypeSynonymInstances,
             ViewPatterns, FlexibleContexts #-}
 #include "cvWrapLEO.h"
 -- | Module for exposing opencv drawing functions. These are meant for quick and dirty marking
@@ -14,6 +14,7 @@
                 ,Drawable(..)
                 -- * Extra drawing operations
                 ,drawLinesOp
+                ,drawPolyLineOp 
                 ,drawBox2Dop 
                 -- * Floodfill operations
                 ,fillOp
@@ -37,12 +38,13 @@
 import CV.Bindings.Drawing
 import Utils.Point
 
-{#import CV.Image#}
+import CV.Image hiding (Complex)
+import qualified CV.Image
 
 import CV.ImageOp
 import Utils.GeometryClass
 import Utils.Rectangle
-import Data.Complex
+import Data.Complex as C
 
 -- | Is the shape filled or just a boundary?
 data ShapeStyle = Filled | Stroked Int
@@ -52,10 +54,10 @@
 styleToCV (Stroked w) = fromIntegral w
 
 -- TODO: The instances in here could be significantly smaller..
--- |Typeclass for images that support elementary drawing operations. 
+-- |Typeclass for images that support elementary drawing operations.
 class Drawable a b where
     -- | Type of the pixel, i.e. Float for a grayscale image and 3-tuple for RGB image.
-    type Color a b :: * 
+    type Color a b :: *
     -- | Put text of certain color to given coordinates. Good size seems to be around 0.5-1.5.
     putTextOp :: (Color a b) -> Float -> String -> (Int,Int) -> ImageOperation a b
     -- | Draw a line between two points.
@@ -66,10 +68,12 @@
     rectOp   :: (BoundingBox bb, Integral (ELBB bb)) => (Color a b) -> Int -> bb -> ImageOperation a b
     -- | Draw a filled polygon
     fillPolyOp :: (Color a b) -> [(Int,Int)] -> ImageOperation a b
+    -- | Draw an ellipse or elliptic arc
+    ellipseOp :: (Color a b) -> Int -> (Int,Int) -> (Int,Int) -> Float -> (Float,Float) -> ImageOperation a b
     ellipseBoxOp :: (Color a b) -> C'CvBox2D -> Int -> Int -> ImageOperation a b
 
 primRectOp (r,g,b) t (bounds -> Rectangle x y w h) = ImgOp $ \i -> do
-                         withGenImage i $ \img -> 
+                         withGenImage i $ \img ->
                               {#call wrapDrawRectangle#} img (fromIntegral x)
                                (fromIntegral y) (fromIntegral $ x+w) (fromIntegral $ y+h)
                                (realToFrac r) (realToFrac g)(realToFrac b)(fromIntegral t)
@@ -77,45 +81,46 @@
 -- | Primitive form of ellipse box. Not typesafe, not for end user.
 primEllipseBox :: (D32,D32,D32,D32) -> C'CvBox2D -> Int -> Int -> ImageOperation c d
 
-primEllipseBox (a,b,c,e) box thickness shift = 
-            ImgOp          $ \i -> 
-            withGenImage i $ \c_img -> 
+primEllipseBox (a,b,c,e) box thickness shift =
+            ImgOp          $ \i ->
+            withGenImage i $ \c_img ->
             with box       $ \c_box ->
-            with (C'CvScalar (rtf a) (rtf b) (rtf c) (rtf 0)) $ \c_color -> 
-             c'wrapEllipseBox c_img c_box c_color (fromIntegral thickness) 8 
+            with (C'CvScalar (rtf a) (rtf b) (rtf c) (rtf 0)) $ \c_color ->
+             c'wrapEllipseBox c_img c_box c_color (fromIntegral thickness) 8
                               (fromIntegral shift)
 
 rtf = realToFrac
 
 instance Drawable RGB D32 where
     type Color RGB D32 = (D32,D32,D32)
-    putTextOp    (r,g,b)  = primTextOp (r,g,b) 
-    lineOp       (r,g,b)  = primLineOp (r,g,b) 
+    putTextOp    (r,g,b)  = primTextOp (r,g,b)
+    lineOp       (r,g,b)  = primLineOp (r,g,b)
     circleOp     (r,g,b)  = primCircleOp (r,g,b)
-    ellipseBoxOp (r,g,b)  = primEllipseBox (r,g,b,0) 
+    ellipseBoxOp (r,g,b)  = primEllipseBox (r,g,b,0)
     rectOp       (r,g,b)  = primRectOp (r,g,b)
-    fillPolyOp   (r,g,b)   = primFillPolyOp (r,g,b)
+    fillPolyOp   (r,g,b)  = primFillPolyOp (r,g,b)
+    ellipseOp    (r,g,b)  = primEllipseOp (r,g,b)
 
 primTextOp (c1,c2,c3) size text (x,y)  = ImgOp $ \img -> do
                                    withGenImage img $ \cimg ->
                                     withCString text $ \(ctext) ->
-                                    {#call wrapDrawText#} cimg ctext (realToFrac size) 
-                                        (fromIntegral x) (fromIntegral y)   
-                                        (realToFrac c1) (realToFrac c2) (realToFrac c3) 
+                                    {#call wrapDrawText#} cimg ctext (realToFrac size)
+                                        (fromIntegral x) (fromIntegral y)
+                                        (realToFrac c1) (realToFrac c2) (realToFrac c3)
 
 primLineOp (c1,c2,c3) t (x,y) (x1,y1) = ImgOp $ \i -> do
-                         withGenImage i $ \img -> 
-                              {#call wrapDrawLine#} img (fromIntegral x) (fromIntegral y) 
-                                                        (fromIntegral x1) (fromIntegral y1) 
-                                                        (realToFrac c1) (realToFrac c2) 
-                                                        (realToFrac c3) (fromIntegral t) 
+                         withGenImage i $ \img ->
+                              {#call wrapDrawLine#} img (fromIntegral x) (fromIntegral y)
+                                                        (fromIntegral x1) (fromIntegral y1)
+                                                        (realToFrac c1) (realToFrac c2)
+                                                        (realToFrac c3) (fromIntegral t)
 
 primCircleOp (c1,c2,c3) (x,y) r s = ImgOp $ \i -> do
-                        when (r>0) $ withGenImage i $ \img -> 
-                              ({#call wrapDrawCircle#} img (fromIntegral x) (fromIntegral y) 
-                                                           (fromIntegral r) 
-                                                           (realToFrac c1) (realToFrac c2) 
-                                                           (realToFrac c3) 
+                        when (r>0) $ withGenImage i $ \img ->
+                              ({#call wrapDrawCircle#} img (fromIntegral x) (fromIntegral y)
+                                                           (fromIntegral r)
+                                                           (realToFrac c1) (realToFrac c2)
+                                                           (realToFrac c3)
                                                            $ styleToCV s)
 
 primFillPolyOp (c1,c2,c3) pts = ImgOp $ \i -> do
@@ -123,58 +128,93 @@
                                   let (xs,ys) = unzip pts
                                   xs' <- newArray $ map fromIntegral xs
                                   ys' <- newArray $ map fromIntegral  ys
-                                  {#call wrapFillPolygon#} img 
-                                       (fromIntegral $ length xs) xs' ys' 
-                                   (realToFrac c1) (realToFrac c2) (realToFrac c3) 
+                                  {#call wrapFillPolygon#} (castPtr img)
+                                       (fromIntegral $ length xs) xs' ys'
+                                   (realToFrac c1) (realToFrac c2) (realToFrac c3)
                                   free xs'
                                   free ys'
 
-instance Drawable DFT D32 where
-   type Color DFT D32 = Complex D32
-   putTextOp (r:+i)   = primTextOp (r,i,0) -- Boy does this feel silly :)
-   lineOp (r:+i)      = primLineOp (r,i,0) 
-   circleOp (r:+i)    = primCircleOp (r,i,0)
-   rectOp   (r:+i)    = primRectOp (r,i,0)
-   ellipseBoxOp (r:+i) = primEllipseBox (r,i,0,0) 
+primEllipseOp (c1,c2,c3) t (x,y) (r1,r2) a (a1,a2) =
+  ImgOp $ \i -> do
+    withGenImage i $ \img ->
+      {#call wrapDrawEllipse#} img
+        (fromIntegral x)
+        (fromIntegral y)
+        (fromIntegral r1)
+        (fromIntegral r2)
+        (realToFrac a)
+        (realToFrac a1)
+        (realToFrac a2)
+        (realToFrac c1)
+        (realToFrac c2)
+        (realToFrac c3)
+        (fromIntegral t)
+
+
+instance Drawable CV.Image.Complex D32 where
+   type Color CV.Image.Complex D32  = Complex D32
+   putTextOp (r:+i)    = primTextOp (r,i,0) -- Boy does this feel silly :)
+   lineOp (r:+i)       = primLineOp (r,i,0)
+   circleOp (r:+i)     = primCircleOp (r,i,0)
+   rectOp   (r:+i)     = primRectOp (r,i,0)
+   ellipseBoxOp (r:+i) = primEllipseBox (r,i,0,0)
    fillPolyOp (r:+i)   = primFillPolyOp (r,i,0)
+   ellipseOp (r:+i)    = primEllipseOp (r,i,0)
 
-instance Drawable GrayScale D32 where
-    type Color GrayScale D32 = D32
+instance Drawable GrayScale D8 where
+    type Color GrayScale D8 = D8
     putTextOp color = primTextOp (color,color,color) 
     lineOp c = primLineOp (c,c,c) 
     circleOp c = primCircleOp (c,c,c)
-    ellipseBoxOp c  = primEllipseBox (c,c,c,0) 
+    ellipseBoxOp c  = primEllipseBox (fromIntegral c,fromIntegral c,fromIntegral c,0) 
     rectOp c = primRectOp (c,c,c)
     fillPolyOp c = primFillPolyOp (c,c,c)
+    ellipseOp  c = primEllipseOp (c,c,c)
 
+instance Drawable GrayScale D32 where
+    type Color GrayScale D32 = D32
+    putTextOp color = primTextOp (color,color,color)
+    lineOp c = primLineOp (c,c,c)
+    circleOp c = primCircleOp (c,c,c)
+    ellipseBoxOp c  = primEllipseBox (c,c,c,0)
+    rectOp c = primRectOp (c,c,c)
+    fillPolyOp c = primFillPolyOp (c,c,c)
+    ellipseOp c   = primEllipseOp (c,c,c)
+
 -- | Flood fill a region of the image
 fillOp :: (Int,Int) -> D32 -> D32 -> D32 -> Bool -> ImageOperation GrayScale D32
-fillOp (x,y) color low high floats = 
+fillOp (x,y) color low high floats =
     ImgOp $ \i -> do
-      withImage i $ \img -> 
-        ({#call wrapFloodFill#} img (fromIntegral x) (fromIntegral y)
+      withImage i $ \img ->
+        ({#call wrapFloodFill#} (castPtr img) (fromIntegral x) (fromIntegral y)
             (realToFrac color) (realToFrac low) (realToFrac high) (toCINT $ floats))
     where
      toCINT False = 0
      toCINT True  = 1
 
 -- | Apply rectOp to an image
-rectangle :: (BoundingBox bb, Integral (ELBB bb), Drawable c d) 
+rectangle :: (BoundingBox bb, Integral (ELBB bb), Drawable c d)
              => Color c d -> Int -> bb -> Image c d
              -> IO (Image c d)
-rectangle color thickness rect i = 
+rectangle color thickness rect i =
     operate (rectOp color thickness rect) i
 
 -- | Apply fillPolyOp to an image
 fillPoly :: Drawable c d => Color c d -> [(Int, Int)] -> Image c d -> IO (Image c d)
 fillPoly c pts i = operate (fillPolyOp c pts) i
 
--- | Draw a polyline
+-- | Draw a line segments
 drawLinesOp :: Drawable c d => Color c d -> Int -> [((Int, Int), (Int, Int))] -> CV.ImageOp.ImageOperation c d
-drawLinesOp color thickness segments = 
-    foldl (#>) nonOp 
+drawLinesOp color thickness segments =
+    foldl (#>) nonOp
      $ map (\(a,b) -> lineOp color thickness a b) segments
 
+-- | Draw a polyline
+drawPolyLineOp :: Drawable c d => Color c d -> Int -> [((Int, Int))] -> CV.ImageOp.ImageOperation c d
+drawPolyLineOp color thickness segments = 
+    foldl (#>) nonOp 
+     $ map (\(a,b) -> lineOp color thickness a b) $ zip segments (tail segments)
+
 -- | Apply drawLinesOp to an image
 drawLines :: Drawable c d => Image c d -> Color c d -> Int -> [((Int, Int), (Int, Int))]
                                 -> IO (Image c d)
@@ -184,13 +224,13 @@
 -- | Draw C'CvBox2D
 drawBox2Dop :: Drawable c d => Color c d -> C'CvBox2D -> ImageOperation c d
 drawBox2Dop color (C'CvBox2D (C'CvPoint2D32f (realToFrac -> x) (realToFrac ->y))
-                             (C'CvSize2D32f  (realToFrac -> w) (realToFrac ->h)) 
-                             (degToRad -> θ)) 
-    = drawLinesOp color 1 (zip corners $ tail (cycle corners)) 
+                             (C'CvSize2D32f  (realToFrac -> w) (realToFrac ->h))
+                             (degToRad -> θ))
+    = drawLinesOp color 1 (zip corners $ tail (cycle corners))
   where
     rot (x,y) = (x * sin (-θ) - y * cos (-θ)
                 ,x * cos (-θ) + y * sin (-θ))
-    corners = map (both round . (+ (x,y)) . rot) 
+    corners = map (both round . (+ (x,y)) . rot)
               [( 0.5*h,  0.5*w)
               ,(-0.5*h,  0.5*w)
               ,(-0.5*h, -0.5*w)
@@ -206,8 +246,8 @@
 
 -- | Apply fillOp to an image
 floodfill :: (Int, Int) -> D32 -> D32 -> D32 -> Bool -> Image GrayScale D32 -> Image GrayScale D32
-floodfill (x,y) color low high floats = 
-    unsafeOperate (fillOp (x,y) color low high floats) 
+floodfill (x,y) color low high floats =
+    unsafeOperate (fillOp (x,y) color low high floats)
 
-            
+
 
diff --git a/CV/Files.hs b/CV/Files.hs
new file mode 100644
--- /dev/null
+++ b/CV/Files.hs
@@ -0,0 +1,28 @@
+module CV.Files where
+
+import CV.Bindings.Files
+import CV.Image
+
+import Foreign.C.String
+import Foreign.Ptr(nullPtr, castPtr)
+import System.IO.Unsafe
+
+readFromTcr :: String -> Image GrayScale D32
+readFromTcr p = unsafePerformIO $ creatingImage $
+  withCString p $ \path ->
+    c'read_from_tcr path nullPtr
+
+readFromTcrRectified :: String -> Image GrayScale D32
+readFromTcrRectified p = unsafePerformIO $ creatingImage $
+  withCString p $ \path ->
+    c'read_from_tcr_rectified path
+
+readFromTcrMag :: String -> Image GrayScale D32
+readFromTcrMag p = unsafePerformIO $ creatingImage $
+  withCString p $ \path ->
+    c'read_from_tcr_mag path
+
+calcDerivativeDirection :: Image GrayScale D32 -> Image GrayScale D32
+calcDerivativeDirection image = unsafePerformIO $ creatingImage $
+  withImage image $ \pimage ->
+    c'calc_derivative_direction (castPtr pimage)
diff --git a/CV/Histogram.chs b/CV/Histogram.chs
--- a/CV/Histogram.chs
+++ b/CV/Histogram.chs
@@ -39,7 +39,6 @@
 -- own number of bins.
 histogram :: [(Image GrayScale D8, Int)] -> Bool -> Maybe (Image GrayScale D8)
                -> I.Histogram
-
 histogram imageBins accumulate mask  = unsafePerformIO $
  I.creatingHistogram $ do
         hist <-  I.emptyUniformHistogramND ds
diff --git a/CV/Image.chs b/CV/Image.chs
--- a/CV/Image.chs
+++ b/CV/Image.chs
@@ -4,8 +4,9 @@
 -- * Basic types
  Image(..)
 , create
-, empty 
-, emptyCopy 
+, empty
+, emptyCopy
+, emptyCopy'
 , cloneImage
 , withClone
 , withCloneValue
@@ -14,7 +15,7 @@
 -- * Colour spaces
 , ChannelOf
 , GrayScale
-, DFT
+, Complex
 , RGB
 , RGBA
 , RGB_Channel(..)
@@ -46,6 +47,7 @@
 -- * Image information
 , ImageDepth
 , Sized(..)
+, biggerThan
 , getArea
 , getChannel
 , getImageChannels
@@ -70,12 +72,14 @@
 , tileImages
 
 -- * Conversions
-, rgbToGray 
+, rgbToGray
 , grayToRGB
-, rgbToLab 
+, rgbToLab
 , bgrToRgb
 , rgbToBgr
+, cloneTo64F
 , unsafeImageTo32F 
+, unsafeImageTo64F 
 , unsafeImageTo8Bit 
 
 -- * Low level access operations
@@ -94,6 +98,8 @@
 -- * Extended error handling
 , setCatch
 , CvException
+, CvSizeError(..)
+, CvIOError(..)
 ) where
 
 import System.Mem
@@ -108,6 +114,7 @@
 import Foreign.Ptr
 import Control.Parallel.Strategies
 import Control.DeepSeq
+
 import CV.Bindings.Error
 
 import Data.Maybe(catMaybes)
@@ -119,22 +126,22 @@
 import Foreign.Storable
 import System.IO.Unsafe
 import Data.Word
-import Data.Complex
-import Data.Complex
+import qualified Data.Complex as C
 import Control.Monad
 import Control.Exception
 import Data.Data
 import Data.Typeable
 
+import Utils.GeometryClass
 
 
 
+
 -- Colorspaces
 
 -- | Single channel grayscale image
 data GrayScale
-
-data DFT
+data Complex
 data RGB
 data RGB_Channel = Red | Green | Blue deriving (Eq,Ord,Enum)
 
@@ -301,11 +308,8 @@
 loadColorImage8 :: FilePath -> IO (Maybe (Image BGR D8))
 loadColorImage8 = unsafeloadUsing imageTo8Bit 1
 
--- | Typeclass for elements with a size, such as images and matrices.
-class Sized a where
-    type Size a :: *
-    getSize :: a -> Size a
 
+
 instance Sized BareImage where
     type Size BareImage = (Int,Int)
    -- getSize :: (Integral a, Integral b) => Image c d -> (a,b)
@@ -428,37 +432,37 @@
     CV_BayerGB2BGR_VNG =63,
     CV_BayerRG2BGR_VNG =64,
     CV_BayerGR2BGR_VNG =65,
-    
+
     CV_BayerBG2RGB_VNG =CV_BayerRG2BGR_VNG,
     CV_BayerGB2RGB_VNG =CV_BayerGR2BGR_VNG,
     CV_BayerRG2RGB_VNG =CV_BayerBG2BGR_VNG,
     CV_BayerGR2RGB_VNG =CV_BayerGB2BGR_VNG,
-    
+
     CV_BGR2HSV_FULL = 66,
     CV_RGB2HSV_FULL = 67,
     CV_BGR2HLS_FULL = 68,
     CV_RGB2HLS_FULL = 69,
-    
+
     CV_HSV2BGR_FULL = 70,
     CV_HSV2RGB_FULL = 71,
     CV_HLS2BGR_FULL = 72,
     CV_HLS2RGB_FULL = 73,
-    
+
     CV_LBGR2Lab     = 74,
     CV_LRGB2Lab     = 75,
     CV_LBGR2Luv     = 76,
     CV_LRGB2Luv     = 77,
-    
+
     CV_Lab2LBGR     = 78,
     CV_Lab2LRGB     = 79,
     CV_Luv2LBGR     = 80,
     CV_Luv2LRGB     = 81,
-    
+
     CV_BGR2YUV      = 82,
     CV_RGB2YUV      = 83,
     CV_YUV2BGR      = 84,
     CV_YUV2RGB      = 85,
-    
+
     CV_COLORCVT_MAX  =100
 };
 #endc
@@ -516,8 +520,8 @@
                                          s <- {#get IplImage->widthStep#} c_i
                                          peek (castPtr (d`plusPtr` (y*(fromIntegral s) +x*sizeOf (0::Word8))):: Ptr Word8)
 
-instance GetPixel (Image DFT D32) where
-    type P (Image DFT D32) = Complex D32
+instance GetPixel (Image Complex D32) where
+    type P (Image Complex D32) = C.Complex D32
     {-#INLINE getPixel#-}
     getPixel (x,y) i = unsafePerformIO $
                         withGenImage i $ \c_i -> do
@@ -527,7 +531,7 @@
                                              fs = sizeOf (undefined :: Float)
                                          re <- peek (castPtr (d`plusPtr` (y*cs + x*2*fs)))
                                          im <- peek (castPtr (d`plusPtr` (y*cs +(x*2+1)*fs)))
-                                         return (re:+im)
+                                         return (re C.:+ im)
 
 -- #define UGETC(img,color,x,y) (((uint8_t *)((img)->imageData + (y)*(img)->widthStep))[(x)*3+(color)])
 instance GetPixel (Image RGB D32) where
@@ -633,7 +637,7 @@
 
 instance CreateImage (Image GrayScale D32) where
     create (w,h) = creatingImage $ {#call wrapCreateImage32F#} (fromIntegral w) (fromIntegral h) 1
-instance CreateImage (Image DFT D32) where
+instance CreateImage (Image Complex D32) where
     create (w,h) = creatingImage $ {#call wrapCreateImage32F#} (fromIntegral w) (fromIntegral h) 2
 instance CreateImage (Image LAB D32) where
     create (w,h) = creatingImage $ {#call wrapCreateImage32F#} (fromIntegral w) (fromIntegral h) 3
@@ -670,12 +674,15 @@
 emptyCopy :: (CreateImage (Image a b)) => Image a b -> (Image a b)
 emptyCopy img = unsafePerformIO $ create (getSize img)
 
+emptyCopy' :: (CreateImage (Image a b)) => Image a b -> IO (Image a b)
+emptyCopy' img = create (getSize img)
+
 -- | Save image. This will convert the image to 8 bit one before saving
 class Save a where
-    save :: FilePath -> a -> IO () 
+    save :: FilePath -> a -> IO ()
 
 instance Save (Image BGR D32) where
-    save filename image = primitiveSave filename (unS . unsafeImageTo8Bit $ image) 
+    save filename image = primitiveSave filename (unS . unsafeImageTo8Bit $ image)
 
 instance Save (Image RGB D32) where
     save filename image = primitiveSave filename (swapRB . unS . unsafeImageTo8Bit $ image)
@@ -687,10 +694,10 @@
     save filename image = primitiveSave filename (unS $ image)
 
 instance Save (Image GrayScale D32) where
-    save filename image = primitiveSave filename (unS . unsafeImageTo8Bit $ image) 
-     
+    save filename image = primitiveSave filename (unS . unsafeImageTo8Bit $ image)
+
 primitiveSave :: FilePath -> BareImage -> IO ()
-primitiveSave filename fpi = do 
+primitiveSave filename fpi = do
        exists <- doesDirectoryExist (takeDirectory filename)
        when (not exists) $ throw (CvIOError $ "Directory does not exist: " ++ (takeDirectory filename))
        withCString  filename $ \name  ->
@@ -802,6 +809,16 @@
                 r <- fun result
                 return r
 
+cloneTo64F :: Image c d -> IO (Image c D64)
+cloneTo64F img = withGenImage img $ \image ->
+                creatingImage
+                 ({#call ensure64F #} image)
+
+unsafeImageTo64F :: Image c d -> Image c D64
+unsafeImageTo64F img = unsafePerformIO $ withGenImage img $ \image ->
+                creatingImage
+                 ({#call ensure64F #} image)
+
 unsafeImageTo32F :: Image c d -> Image c D32
 unsafeImageTo32F img = unsafePerformIO $ withGenImage img $ \image ->
                 creatingImage
@@ -925,13 +942,13 @@
                                          let cs = fromIntegral s
                                              fs = sizeOf (undefined :: Float)
                                          poke (castPtr (d`plusPtr` (y*cs +x*3*fs)))     b
-                                         poke (castPtr (d`plusPtr` (y*cs +(x*3+1)*fs))) g 
+                                         poke (castPtr (d`plusPtr` (y*cs +(x*3+1)*fs))) g
                                          poke (castPtr (d`plusPtr` (y*cs +(x*3+2)*fs))) r
 
-instance SetPixel (Image DFT D32) where
-    type SP (Image DFT D32) = Complex D32
+instance SetPixel (Image Complex D32) where
+    type SP (Image Complex D32) = C.Complex D32
     {-#INLINE setPixel#-}
-    setPixel (x,y) (re:+im) image = withGenImage image $ \c_i -> do
+    setPixel (x,y) (re C.:+ im) image = withGenImage image $ \c_i -> do
                              d <- {#get IplImage->imageData#} c_i
                              s <- {#get IplImage->widthStep#} c_i
                              let cs = fromIntegral s
@@ -958,13 +975,13 @@
 --  between images. Images are assumed to be the same size (determined by the first image)
 montage :: (CreateImage (Image GrayScale D32)) => (Int,Int) -> Int -> [Image GrayScale D32] -> Image GrayScale D32
 montage (u',v') space' imgs
-    | u'*v' /= (length imgs) = error ("Montage mismatch: "++show (u,v, length imgs))
+    | u'*v' < (length imgs) = error ("Montage mismatch: "++show (u,v, length imgs))
     | otherwise              = resultPic
     where
      space = fromIntegral space'
      (u,v) = (fromIntegral u', fromIntegral v')
      (rw,rh) = (u*xstep,v*ystep)
-     (w,h) = getSize (head imgs)
+     (w,h) = foldl (\(mx,my) (x,y) -> (max mx x, max my y)) (0,0) $ map getSize imgs
      (xstep,ystep) = (fromIntegral space + w,fromIntegral space + h)
      edge = space`div`2
      resultPic = unsafePerformIO $ do
@@ -978,16 +995,18 @@
      deriving (Show, Typeable)
 
 data CvIOError = CvIOError String deriving (Show,Typeable)
+data CvSizeError = CvSizeError String deriving (Show,Typeable)
 
 instance Exception CvException
 instance Exception CvIOError
+instance Exception CvSizeError
 
 setCatch = do
    let catch i cstr1 cstr2 cstr3 j = do
          func <- peekCString cstr1
          msg  <- peekCString cstr2
          file <- peekCString cstr3
-         throw (CvException (fromIntegral i) func msg file (fromIntegral j)) 
+         throw (CvException (fromIntegral i) func msg file (fromIntegral j))
          return 0
    cb <- mk'CvErrorCallback catch
    c'cvRedirectError cb nullPtr nullPtr
diff --git a/CV/ImageOp.hs b/CV/ImageOp.hs
--- a/CV/ImageOp.hs
+++ b/CV/ImageOp.hs
@@ -20,6 +20,10 @@
 -- |Apply image operation to a Copy of an image
 img <# op = unsafeOperate op img
 
+blitOp img pos = ImgOp $ \i -> blit img i pos
+
+setPixelOp pos v = ImgOp $ \i -> setPixel pos v i
+
 -- motivating example:
 -- >>> hop i = stretchHistogram $ i #- gaussian (5,5)
 -- allocates two extra images
diff --git a/CV/Matrix.hs b/CV/Matrix.hs
--- a/CV/Matrix.hs
+++ b/CV/Matrix.hs
@@ -5,7 +5,7 @@
     (
     Exists(..),
     Matrix, emptyMatrix, fromFunction, fromList,toList,toRows,toCols,get,put,withMatPtr
-    , transpose, mxm, rodrigues2, identity
+    , transpose, mxm, invert,  rodrigues2, identity
     )where
 
 
@@ -135,6 +135,16 @@
                  withForeignPtr f_m $ \c_m ->
                   withForeignPtr f_c $ \c_c -> c'cvRodrigues2 c_m c_c nullPtr
                  return res
+
+
+-- | Matrix inversion
+invert :: (Exists (Matrix a), Args (Matrix a) ~ Size (Matrix a)) => Matrix a -> Matrix a
+invert m@(Matrix f_m) = unsafePerformIO $ do
+                 res@(Matrix f_c) <- create (getSize m)
+                 withForeignPtr f_m $ \c_m ->
+                  withForeignPtr f_c $ \c_c -> c'cvInvert c_m c_c (fromIntegral . fromEnum $ c'CV_LU)
+                 return res
+
 
 
 -- | Ordinary matrix multiplication
diff --git a/CV/Operations.hs b/CV/Operations.hs
--- a/CV/Operations.hs
+++ b/CV/Operations.hs
@@ -1,14 +1,18 @@
+{-# LANGUAGE ScopedTypeVariables #-}
 module CV.Operations
 ( clear
 , set
+, expand
 , NormType(..)
 , normalize
 , unitNormalize
 , unitStretch
 , logNormalize
+, cartToPolar
 ) where
 
 import CV.Bindings.Core
+import CV.Bindings.ImgProc
 import CV.Bindings.Types
 import CV.Image
 import CV.ImageMath as IM
@@ -29,6 +33,12 @@
     c'wrapSetAll (castPtr i_ptr) (realToFrac v) nullPtr
   return i
 
+
+expand :: (Int,Int,Int,Int) -> Image d c -> Image d c
+expand (top,bottom,left,right) i = unsafePerformIO $
+  copyMakeBorder i top bottom left right BorderReplicate 0
+
+
 data NormType =
   NormC |
   NormL1 |
@@ -64,9 +74,9 @@
 normalize a b t src =
   unsafePerformIO $ do
     withCloneValue src $ \clone ->
-      withImage src $ \si ->
-        withImage clone $ \ci -> do
-          c'cvNormalize (castPtr si) (castPtr ci) (realToFrac a) (realToFrac b) (cNormType t) nullPtr
+      withGenImage src $ \si ->
+        withGenImage clone $ \ci -> do
+          c'cvNormalize si ci (realToFrac a) (realToFrac b) (cNormType t) nullPtr
           return clone
 
 unitNormalize i
@@ -78,3 +88,16 @@
 unitStretch i = normalize 0 1 NormMinMax i
 
 logNormalize = unitNormalize . IM.log . (1 |+)
+
+cartToPolar :: (Image GrayScale D32, Image GrayScale D32) -> (Image GrayScale D32, Image GrayScale D32)
+cartToPolar (x,y) = unsafePerformIO $ do
+  r::(Image GrayScale D32) <- create (w, h)
+  a::(Image GrayScale D32) <- create (w, h)
+  withImage x $ \px ->
+    withImage y $ \py ->
+      withImage r $ \pr ->
+        withImage a $ \pa -> do
+          c'cvCartToPolar (castPtr px) (castPtr py) (castPtr pr) (castPtr pa) (fromIntegral 0)
+          return (r,a)
+  where
+    (w,h) = getSize x
diff --git a/CV/Projection.hs b/CV/Projection.hs
new file mode 100644
--- /dev/null
+++ b/CV/Projection.hs
@@ -0,0 +1,12 @@
+module CV.Projection where
+
+import CV.Bindings.Projection
+import CV.Image
+
+import Foreign.Ptr
+import System.IO.Unsafe
+
+projectPolar :: Image c d -> Image GrayScale D32
+projectPolar i = unsafePerformIO $ creatingImage $
+  withImage i $ \i_ptr ->
+    c'project_polar (castPtr i_ptr)
diff --git a/CV/Sampling.hs b/CV/Sampling.hs
--- a/CV/Sampling.hs
+++ b/CV/Sampling.hs
@@ -1,8 +1,9 @@
 module CV.Sampling where
 
 import CV.Image
-import System.Random
+import Control.Monad.Primitive
 import Control.Monad
+import System.Random.MWC
 
 import Foreign.C.Types
 import qualified CV.ImageMath as IM
@@ -63,11 +64,11 @@
                                          adjust (x,y) = (x-w`div`2
                                                         ,y-h`div`2)
 
--- Make a random selections in IO monad
-randomSelect lst = randomRIO (0,length lst -1) >>= \x ->
-                              return (lst !! x)
+---- Make a random selections in IO monad
+--randomSelect lst = randomRIO (0,length lst -1) >>= \x ->
+--                              return (lst !! x)
                               
-select k lst = sequence $ replicate k (randomSelect lst)
+-- select k lst = sequence $ replicate k (randomSelect lst)
 
 -- Discard coords around image borders. Useful for safely picking patches
 discardAroundEdges (iw,ih) (vb,hb) coords = filter inRange coords
@@ -89,28 +90,28 @@
     where coords = getCoordsFromMarks marks
 
 
----- Get some random image patches
---randomPatches size count image = do
---    coords <- replicateM count $ randomCoord (w,h)
---    return $ getPatches size coords image
--- where
---    (pwidth,pheight) = size
---    (iwidth,iheight) = getSize image
---    (w,h) = (iwidth - pwidth , iheight-pheight) 
+-- Get some random image patches
+randomPatches size count image gen = do
+    coords <- replicateM count $ randomCoord (w,h) gen
+    return $ getPatches size coords image
+ where
+    (pwidth,pheight) = size
+    (iwidth,iheight) = getSize image
+    (w,h) = (iwidth - pwidth , iheight-pheight) 
 
----- Get some random pixels from image
---randomPixels count image = do
---   coords <- replicateM count $ randomCoord size
---   return $ map (flip getPixel $ image) $ coords 
--- where
---  size = getSize image
+-- Get some random pixels from image
+-- randomPixels count image = do
+--    coords <- replicateM count $ randomCoord size
+--    return $ map (flip getPixel $ image) $ coords 
+--  where
+--   size = getSize image
 
 ---- Get some random coords from image
 --randomCoords :: MonadRandom m => Int -> (Int,Int) -> m [(Int,Int)]
 --randomCoords count area = replicateM count $ randomCoord area
 
---randomCoord :: MonadRandom m => (Int,Int) -> m (Int,Int)
---randomCoord (w,h) = do
---            x <- (getRandomR (0::Int,fromIntegral $ w-1))
---            y <- (getRandomR (0::Int,fromIntegral $ h-1))
---            return (x,y) 
+randomCoord :: PrimMonad m => (Int,Int) -> Gen (PrimState m) -> m (Int,Int)
+randomCoord (w,h) g = do
+            x <- uniformR (0::Int,fromIntegral $ w-1) g
+            y <- uniformR (0::Int,fromIntegral $ h-1) g
+            return (x,y) 
diff --git a/CV/TemplateMatching.chs b/CV/TemplateMatching.chs
--- a/CV/TemplateMatching.chs
+++ b/CV/TemplateMatching.chs
@@ -16,6 +16,7 @@
 
 {#import CV.Image#}
 import System.IO.Unsafe
+import Control.Exception
 
 getTemplateMap image template = unsafePerformIO $
 	   withImage image $ \cvimg ->
@@ -38,6 +39,7 @@
 
 simpleTemplateMatch :: MatchType -> Image GrayScale D32 -> Image GrayScale D32 -> ((Int,Int),Double)
 simpleTemplateMatch mt image template 
+ | image `biggerThan` template
 	= unsafePerformIO $ do
 	   withImage image $ \cvimg ->
 	    withImage template $ \cvtemp ->
@@ -49,9 +51,11 @@
 			y <- peek ptrinty;
 			v <- peek ptrdblval;
 		    return ((fromIntegral x,fromIntegral y),realToFrac v); }
+ | otherwise = throw $ CvSizeError "simpleTemplateMatch: template is bigger than the image"
 
 matchTemplate :: MatchType-> Image GrayScale D32 -> Image GrayScale D32 -> Image GrayScale D32 
-matchTemplate mt image template = unsafePerformIO $ do
+matchTemplate mt image template 
+ | image `biggerThan` template = unsafePerformIO $ do
      let isize = getSize image
          tsize = getSize template
          size  = isize - tsize + (1,1) 
@@ -61,6 +65,7 @@
        withGenImage res $ \cresult -> 
         {#call cvMatchTemplate#} cimg ctempl cresult (fromIntegral . fromEnum $ mt)
      return res
+ | otherwise = throw $ CvSizeError $ "MatchTemplate: template larger than the image"
 
 
 -- | Perform subpixel template matching using intensity interpolation
diff --git a/CV/Thresholding.hs b/CV/Thresholding.hs
--- a/CV/Thresholding.hs
+++ b/CV/Thresholding.hs
@@ -4,6 +4,7 @@
   -- * Interfaces to OpenCV functions
   ThresholdType(..)
 , threshold
+, thresholdInPlace
 , thresholdOtsu
 , AdaptiveType(..)
 , adaptiveThreshold
@@ -24,6 +25,7 @@
 import System.IO.Unsafe
 import CV.Sampling
 import Utils.List
+import Unsafe.Coerce
 import Data.List
 import CV.Histogram
 import CV.Bindings.ImgProc
@@ -95,6 +97,14 @@
           c'cvThreshold (castPtr pimage) (castPtr presult) (realToFrac tval)
             (realToFrac (maxval image)) (cThresholdType ttype)
           return result
+
+-- TODO: Convert into imageOperation
+thresholdInPlace :: (MaxVal d) => ThresholdType -> Double -> Image GrayScale d -> IO (Image GrayScale D8)
+thresholdInPlace ttype tval image = do
+      withImage image $ \pimage ->
+          c'cvThreshold (castPtr pimage) (castPtr pimage) (realToFrac tval)
+            (realToFrac (maxval image)) (cThresholdType ttype)
+      return (unsafeCoerce image)
 
 -- | Thresholds a grayscale image using the otsu method according to the
 --   selected type. Threshold value is selected automatically, and only 8-bit
diff --git a/CV/Transforms.chs b/CV/Transforms.chs
--- a/CV/Transforms.chs
+++ b/CV/Transforms.chs
@@ -10,7 +10,8 @@
 import System.IO.Unsafe
 {#import CV.Image#}
 import CV.ImageMathOp
-import CV.Matrix as M
+import qualified CV.Matrix as M
+import CV.Matrix (Matrix,withMatPtr)
 
 -- |Since DCT is valid only for even sized images, we provide a
 -- function to crop images to even sizes.
@@ -116,7 +117,16 @@
     = unsafePerformIO $ 
        withImage img $ \cimg -> creatingImage $ {#call wrapPerspective#} cimg a1 a2 a3 a4 a5 a6 a7 a8 a9
 
+perspectiveTransform' :: (CreateImage (Image c d)) => Matrix Float -> Image c d -> (Int,Int)-> Image c d
+perspectiveTransform' mat img size
+    = unsafePerformIO $ do
+       r <- create  size
+       withImage img $ \c_img ->
+         withMatPtr mat $ \c_mat ->
+         withImage r $ \c_r -> {#call wrapWarpPerspective#} (castPtr c_img) (castPtr c_r) (castPtr c_mat)
+       return r
 
+
 -- |Find a homography between two sets of points in. The resulting 3x3 matrix is returned as a list.
 getHomography srcPts dstPts = 
     unsafePerformIO $ withArray src $ \c_src ->
@@ -230,6 +240,21 @@
      result = pyrDown img 
      (w,h)  = getSize result 
 
+-- | Enlargen the image so that its size is a power of two.
+minEnlarge :: Image GrayScale D32 -> Image GrayScale D32
+minEnlarge i = enlargeShadow (min (ceiling (logBase 2 (f w))) (ceiling (logBase 2 (f h)))) i
+    where 
+     f = fromIntegral
+     (w,h) = getSize i
+
+-- | Calculate an infinite gaussian pyramid of an image while keeping track of
+--   various corner cases and gotchas.
+gaussianPyramid :: Image GrayScale D32 -> [Image GrayScale D32]
+gaussianPyramid = iterate pyrDown' . minEnlarge
+    where 
+     pyrDown' i = let (w,h) = getSize i
+                  in if (w`div`2) <=1 || (h`div`2) <= 1 then i else pyrDown i
+
 -- |Calculate the laplacian pyramid of an image up to the nth level.
 --  Notice that the image size must be divisible by 2^n or opencv 
 --  will abort (TODO!)
@@ -247,7 +272,8 @@
   --   safeAdd x y = sameSizePad y x #+ y  
 
 -- TODO: Could have wider type
--- |Enlarge image so, that it's size is divisible by 2^n 
+-- |Enlargen the image so that its size is divisible by 2^n. Fill the area
+--  outside the image with black.
 enlarge :: Int -> Image GrayScale D32 -> Image GrayScale D32
 enlarge n img =  unsafePerformIO $ do
                    i <- I.create (w2,h2)
@@ -259,6 +285,20 @@
      pad x = x + (np - x `mod` np)
      np = 2^n
 
+-- | Enlargen the image so that its size is is divisible by 2^n. Replicate
+--   the border of the image.
+enlargeShadow :: Int -> Image GrayScale D32 -> Image GrayScale D32
+enlargeShadow n img =  unsafePerformIO $ do
+                   i <- create (w2,h2)
+                   withImage img $ \c_img -> 
+                    withImage i  $ \c_i   -> {#call blitShadow#} c_i c_img 
+                   return i
+    where
+     (w,h) = getSize img
+     (w2,h2) = (pad w, pad h)
+     pad x = x + (np - x `mod` np)
+     np = 2^n
+
 #c
 enum DistanceType {
      C =  CV_DIST_C
@@ -267,6 +307,15 @@
 };
 #endc
 {#enum DistanceType {}#}
+#ifdef OpenCV24
+#c
+enum LabelType {
+     CCOMP = CV_DIST_LABEL_CCOMP
+    ,PIXEL = CV_DIST_LABEL_PIXEL
+};
+#endc
+{#enum LabelType {}#}
+#endif
 
 -- |Mask sizes accepted by distanceTransform
 data MaskSize = M3 | M5 deriving (Eq,Ord,Enum,Show)
@@ -281,6 +330,11 @@
                                   (fromIntegral . fromEnum $ dtype) 
                                   (fromIntegral . fromEnum $ maskSize)
                                    nullPtr nullPtr
+#ifdef OpenCV24
+                                  (fromIntegral . fromEnum $ CCOMP)
+#endif
+
     return result
+
     -- TODO: Add handling for labels
     -- TODO: Add handling for custom masks
diff --git a/Utils/GeometryClass.hs b/Utils/GeometryClass.hs
--- a/Utils/GeometryClass.hs
+++ b/Utils/GeometryClass.hs
@@ -60,3 +60,14 @@
    type ELS a :: *
    startEnd :: a -> ((ELS a, ELS a),(ELS a, ELS a))
 
+-- | Typeclass for elements with a size, such as images and matrices.
+class Sized a where
+    type Size a :: *
+    getSize :: a -> Size a
+
+biggerThan :: (Sized a, Sized b, Size a~(Int,Int), Size b ~Size a) => a -> b -> Bool
+biggerThan a b = w1>=w2 && h1>=h2
+    where
+     (w1,h1) = getSize a
+     (w2,h2) = getSize b
+
diff --git a/Utils/Pointer.hs b/Utils/Pointer.hs
--- a/Utils/Pointer.hs
+++ b/Utils/Pointer.hs
@@ -1,6 +1,7 @@
 module Utils.Pointer where 
 import Foreign.C.Types
-import Foreign.ForeignPtr
+import Foreign.ForeignPtr hiding (unsafeForeignPtrToPtr)
+import Foreign.ForeignPtr.Unsafe
 import Foreign.Ptr
 import Foreign.Marshal.Array
 
diff --git a/cbits/cvFiles.c b/cbits/cvFiles.c
new file mode 100644
--- /dev/null
+++ b/cbits/cvFiles.c
@@ -0,0 +1,1963 @@
+#include "cvFiles.h"
+#include <opencv2/imgproc/imgproc_c.h>
+#include <opencv2/highgui/highgui_c.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+
+typedef unsigned int uint_t;
+
+const float epsilon = 0.001;
+
+int *mask_idx;
+int mask_size;
+
+float calc_sum(float *pos, uint_t *count)
+{
+  int *p, *e;
+  float s, v;
+  uint_t c;
+
+  c = 0;
+  s = 0;
+  e = &mask_idx[mask_size];
+  for (p = mask_idx; p < e; p++) {
+    v = *(pos + *p);
+    if (v > epsilon) {
+      s += v;
+      c++;
+    }
+  }
+  *count = c;
+  return s;
+}
+
+float calc_sum2(float *pos, uint_t *count)
+{
+  int *p, *e;
+  float s, v;
+  uint_t c;
+
+  c = 0;
+  s = 0;
+  e = &mask_idx[mask_size];
+  for (p = mask_idx; p < e; p++) {
+    v = *(pos + *p);
+    if (v > epsilon) {
+      s += v*v;
+      c++;
+    }
+  }
+  *count = c;
+  return s;
+}
+
+float calc_sum9(float *pos, size_t stride)
+{
+  float *pos1, *pos2, *pos3, *pos4, *pos5;
+  pos1 = pos - 1 * stride - 1;
+  pos2 = pos - 1;
+  pos3 = pos + stride - 1;
+  return
+    *pos1 + *(pos1+1) + *(pos1+2) +
+    *pos2 + *(pos2+1) + *(pos2+2) +
+    *pos3 + *(pos3+1) + *(pos3+2);
+}
+
+float calc_dx(float *pos)
+{
+  float d, p1, p2, n1, n2;
+
+  p1 = *(pos + 1);
+  p2 = *(pos + 2);
+  n1 = *(pos - 1);
+  n2 = *(pos - 2);
+
+  d = 0;
+  if (p1 > epsilon && n1 > epsilon) {
+    d = p1 - n1;
+  }
+  if (p2 > epsilon && n2 > epsilon) {
+    d += (2*p2 - 2*n2);
+  }
+
+  return d;
+}
+
+float calc_dy(float *pos, uint_t stride)
+{
+  float d, p1, p2, n1, n2;
+
+  p1 = *(pos + stride);
+  p2 = *(pos + 2 * stride);
+  n1 = *(pos - stride);
+  n2 = *(pos - 2 * stride);
+
+  d = 0;
+  if (p1 > epsilon && n1 > epsilon) {
+    d = p1 - n1;
+  }
+  if (p2 > epsilon && n2 > epsilon) {
+    d += (2*p2 - 2*n2);
+  }
+
+  return d;
+}
+
+#define C_PI 3.14159265358979323846
+const float pi = C_PI;
+const float dir_diff_small = C_PI / 8;
+const float dir_diff_large = C_PI / 4;
+
+float direction_difference(float a, float b)
+{
+  float diff = a - b;
+  if (diff < 0) diff = -diff;
+  if (diff > pi) diff = 2*pi - diff;
+  return diff;
+}
+/*
+bool is_same_direction(float a, float b, float diff)
+{
+  float lower_bound, upper_bound;
+
+  lower_bound = a - diff;
+  upper_bound = b + diff;
+  
+  if (lower_bound < -pi) {
+    return ((((lower_bound + 2*pi) < b) || -pi < b) && b < upper_bound)
+  }
+  else
+  if (upper_bound > pi) {
+    return (lower_bound < b && (b < pi || b < (upper_bound - 2*pi))); 
+  }
+  else {
+    return (lower_bound < b && b < upper_bound);
+  }
+  return false;
+}
+*/
+
+float calc_directionality_1(float *pos, uint_t stride)
+{
+  float result, dir_c, dir_a, dir_b, dir_diff_s1, dir_diff_s2, value, diff_s, diff_c;
+  
+  result = 0;
+  dir_c = *pos;
+  if (direction_difference(dir_c, pi) < dir_diff_small || direction_difference(dir_c, 0) < dir_diff_small) {
+    dir_a = *(pos - stride);
+    dir_b = *(pos + stride);
+    
+    value = *(pos - stride - 1);
+    dir_diff_s1 = direction_difference(dir_c, value);
+    value = *(pos - 1);
+    dir_diff_s1 += direction_difference(dir_c, value);
+    value = (*pos + stride - 1);
+    dir_diff_s1 += direction_difference(dir_c, value);
+    dir_diff_s1 /= 3;
+    
+    value = *(pos - stride + 1);
+    dir_diff_s2 = direction_difference(dir_c, value);
+    value = *(pos + 1);
+    dir_diff_s2 += direction_difference(dir_c, value);
+    value = *(pos + stride + 1);
+    dir_diff_s2 += direction_difference(dir_c, value);
+    dir_diff_s2 /= 3;
+    
+    if (dir_diff_s1 > dir_diff_s2) {
+      diff_s = dir_diff_s1;
+    }
+    else {
+      diff_s = dir_diff_s2;
+    }
+    
+    if (diff_s > dir_diff_small) {
+      diff_c = direction_difference(dir_c, dir_a);
+      if (diff_c < dir_diff_small) {
+        result += 1.0;
+      }
+      else
+      if (diff_c < diff_s) {
+        result += 0.5;
+      }
+      diff_c = direction_difference(dir_c, dir_b);
+      if (diff_c < dir_diff_small) {
+        result += 1.0;
+      }
+      else
+      if (diff_c < diff_s) {
+        result += 0.5;
+      }
+    }
+  }
+  return result;
+}
+
+float calc_directionality_2(float *pos, uint_t stride)
+{
+  float result, dir_c, dir_a, dir_b, dir_diff_s1, dir_diff_s2, value, diff_s, diff_c;
+  
+  result = 0;
+  dir_c = *pos;
+  if (direction_difference(dir_c, 3*pi/4) < dir_diff_small || direction_difference(dir_c, -pi/4) < dir_diff_small) {
+    dir_a = *(pos - stride + 1);
+    dir_b = *(pos + stride - 1);
+    
+    value = *(pos - stride - 1);
+    dir_diff_s1 = direction_difference(dir_c, value);
+    value = *(pos - stride);
+    dir_diff_s1 += direction_difference(dir_c, value);
+    value = (*pos - 1);
+    dir_diff_s1 += direction_difference(dir_c, value);
+    dir_diff_s1 /= 3;
+    
+    value = *(pos + stride + 1);
+    dir_diff_s2 = direction_difference(dir_c, value);
+    value = *(pos + stride);
+    dir_diff_s2 += direction_difference(dir_c, value);
+    value = *(pos + 1);
+    dir_diff_s2 += direction_difference(dir_c, value);
+    dir_diff_s2 /= 3;
+    
+    if (dir_diff_s1 > dir_diff_s2) {
+      diff_s = dir_diff_s1;
+    }
+    else {
+      diff_s = dir_diff_s2;
+    }
+    
+    if (diff_s > dir_diff_small) {
+      diff_c = direction_difference(dir_c, dir_a);
+      if (diff_c < dir_diff_small) {
+        result += 1.0;
+      }
+      else
+      if (diff_c < diff_s) {
+        result += 0.5;
+      }
+      diff_c = direction_difference(dir_c, dir_b);
+      if (diff_c < dir_diff_small) {
+        result += 1.0;
+      }
+      else
+      if (diff_c < diff_s) {
+        result += 0.5;
+      }
+    }
+  }
+  return result;
+}
+
+float calc_directionality_3(float *pos, uint_t stride)
+{
+  float result, dir_c, dir_a, dir_b, dir_diff_s1, dir_diff_s2, value, diff_s, diff_c;
+  
+  result = 0;
+  dir_c = *pos;
+  if (direction_difference(dir_c, pi/2) < dir_diff_small || direction_difference(dir_c, -pi/2) < dir_diff_small) {
+    dir_a = *(pos + 1);
+    dir_b = *(pos - 1);
+    
+    value = *(pos - stride - 1);
+    dir_diff_s1 = direction_difference(dir_c, value);
+    value = *(pos - stride);
+    dir_diff_s1 += direction_difference(dir_c, value);
+    value = (*pos - stride + 1);
+    dir_diff_s1 += direction_difference(dir_c, value);
+    dir_diff_s1 /= 3;
+    
+    value = *(pos + stride - 1);
+    dir_diff_s2 = direction_difference(dir_c, value);
+    value = *(pos + stride);
+    dir_diff_s2 += direction_difference(dir_c, value);
+    value = *(pos + stride + 1);
+    dir_diff_s2 += direction_difference(dir_c, value);
+    dir_diff_s2 /= 3;
+    
+    if (dir_diff_s1 > dir_diff_s2) {
+      diff_s = dir_diff_s1;
+    }
+    else {
+      diff_s = dir_diff_s2;
+    }
+    
+    if (diff_s > dir_diff_small) {
+      diff_c = direction_difference(dir_c, dir_a);
+      if (diff_c < dir_diff_small) {
+        result += 1.0;
+      }
+      else
+      if (diff_c < diff_s) {
+        result += 0.5;
+      }
+      diff_c = direction_difference(dir_c, dir_b);
+      if (diff_c < dir_diff_small) {
+        result += 1.0;
+      }
+      else
+      if (diff_c < diff_s) {
+        result += 0.5;
+      }
+    }
+  }
+  return result;
+}
+
+float calc_directionality_4(float *pos, uint_t stride)
+{
+  float result, dir_c, dir_a, dir_b, dir_diff_s1, dir_diff_s2, value, diff_s, diff_c;
+  
+  result = 0;
+  dir_c = *pos;
+  if (direction_difference(dir_c, pi/4) < dir_diff_small || direction_difference(dir_c, -3*pi/4) < dir_diff_small) {
+    dir_a = *(pos + stride + 1);
+    dir_b = *(pos - stride - 1);
+    
+    value = *(pos - stride);
+    dir_diff_s1 = direction_difference(dir_c, value);
+    value = *(pos - stride + 1);
+    dir_diff_s1 += direction_difference(dir_c, value);
+    value = (*pos + 1);
+    dir_diff_s1 += direction_difference(dir_c, value);
+    dir_diff_s1 /= 3;
+    
+    value = *(pos - 1);
+    dir_diff_s2 = direction_difference(dir_c, value);
+    value = *(pos + stride - 1);
+    dir_diff_s2 += direction_difference(dir_c, value);
+    value = *(pos + stride);
+    dir_diff_s2 += direction_difference(dir_c, value);
+    dir_diff_s2 /= 3;
+    
+    if (dir_diff_s1 > dir_diff_s2) {
+      diff_s = dir_diff_s1;
+    }
+    else {
+      diff_s = dir_diff_s2;
+    }
+    
+    if (diff_s > dir_diff_small) {
+      diff_c = direction_difference(dir_c, dir_a);
+      if (diff_c < dir_diff_small) {
+        result += 1.0;
+      }
+      else
+      if (diff_c < diff_s) {
+        result += 0.5;
+      }
+      diff_c = direction_difference(dir_c, dir_b);
+      if (diff_c < dir_diff_small) {
+        result += 1.0;
+      }
+      else
+      if (diff_c < diff_s) {
+        result += 0.5;
+      }
+    }
+  }
+  return result;
+}
+
+IplImage *read_from_tcr(const char *path, uint_t **timestamp_array)
+{
+  IplImage *dst;
+  float *dst_data, *dst_pos;
+  uint_t dst_stride, *timestamps, timestamp;
+  size_t width, height, file_size, read_size, line_size, timestamp_size;
+  size_t line_count, timestamp_count;
+  CvSize size;
+  FILE *input_file;
+
+  width = 2592;
+  line_size = width * sizeof(float);
+  timestamp_size = sizeof(uint_t);
+  dst = NULL;
+
+  input_file = fopen(path, "rb");
+  if (input_file == NULL) {
+    printf("Error: unable to open file %s\n", path);
+    return NULL;
+  }
+
+  /* seek file end to determine image size */
+  fseek(input_file , 0 , SEEK_END);
+  file_size = ftell(input_file);
+  rewind(input_file);
+  height = (size_t)(file_size / (line_size + timestamp_size));
+
+  size.width = width;
+  size.height = height;
+  dst = cvCreateImage(size, IPL_DEPTH_32F, 1);
+  if (dst == NULL) {
+    printf("Error: unable to create image\n");
+    return NULL;
+  }
+
+  dst_data = (float*)dst->imageData;
+  dst_stride = (uint_t)(dst->widthStep / sizeof(float));
+
+  if (timestamp_array != NULL) {
+    timestamps = (uint_t *)malloc(height * timestamp_size);
+    if (timestamps == NULL) {
+      printf("Error: unable to allocate timestamp array\n");
+      return NULL;
+    }
+    *timestamp_array = timestamps;
+  }
+  else {
+    timestamps = NULL;
+  }
+
+  line_count = 0;
+  timestamp_count = 0;
+
+  dst_pos = dst_data;
+  while (1) {
+    /* read line timestamp to use for rectification */
+    read_size = fread((void*)&timestamp, timestamp_size, 1, input_file);
+    if (read_size < 1) {
+      if (line_count < height) {
+        printf("Error: unexpected end of file after line %d\n", line_count);
+        return dst;
+      }
+      break;
+    }
+    else {
+      if (timestamps != NULL) {
+        timestamps[timestamp_count++] = timestamp;
+      }
+    }
+
+    /* read data line */
+    read_size = fread((void*)dst_pos, sizeof(float), width, input_file);
+    if (read_size < width) {
+      printf("Error unexpected end of file after line %d\n", line_count);
+      return dst;
+    }
+    else {
+      line_count += 1;
+      dst_pos += dst_stride;
+    }
+  }
+  fclose(input_file);
+
+  return dst;
+}
+
+IplImage *rectify_tcr(IplImage *src, uint_t *timestamps, uint_t startx, uint_t endx, uint_t starty, uint_t endy)
+{
+  IplImage *dst;
+  float *src_data, *dst_data, *src_pos, *dst_pos, tratio;
+  uint_t src_stride, dst_stride, tstart, tend, tdiff, *timestamps_accum;
+  size_t height, new_width, new_height, row_size;
+  int i, j, result;
+  CvSize size;
+
+  new_width = endx - startx;
+  new_height = new_width;
+
+  tstart = timestamps[starty];
+  tend = timestamps[endy];
+  tdiff = tend - tstart;
+  tratio = (float)tdiff / (float)new_height;
+
+  /* create a lookup table for rectified image time stamps */
+  timestamps_accum = (uint_t *)malloc(new_height * sizeof(uint_t));
+  if (timestamps_accum == NULL) {
+    printf("Error: rectify_tcr: unable to allocate cumulative timestamp array\n");
+    return NULL;
+  }
+
+  timestamps_accum[0] = tstart;
+  for (i = 1; i < new_height; i++) {
+    timestamps_accum[i] = tstart + (uint_t)(tratio * i);
+  }
+
+  size.width = new_width;
+  size.height = new_height;
+  dst = cvCreateImage(size, IPL_DEPTH_32F, 1);
+  if (dst == NULL) {
+    printf("Error: rectify_tcr: unable to create image\n");
+    return NULL;
+  }
+  height = src->height;
+  src_data = (float*)src->imageData;
+  src_stride = (uint_t)(src->widthStep / sizeof(float));
+  dst_data = (float*)dst->imageData;
+  dst_stride = (uint_t)(dst->widthStep / sizeof(float));
+
+  src_pos = src_data + starty * src_stride + startx;
+  j = starty;
+  row_size = new_width * sizeof(float);
+  dst_pos = dst_data;
+  for (i = 0; i < new_height; i++, dst_pos += dst_stride) {
+    if (timestamps[j] < timestamps_accum[i]) {
+      j++;
+      if (j == height) {
+        printf("Error: rectify_tcr: row overflow\n");
+        break;
+      }
+      src_pos += src_stride;
+    }
+    memcpy(dst_pos, src_pos, row_size);
+  }
+
+  free(timestamps_accum);
+  return dst;
+}
+
+IplImage *derivative_direction(IplImage *src, float ignore_val, float ignore_eps)
+{
+  IplImage *dst, *tr, *ta;
+  float *src_data, *src_pos, *dst_data, *dst_pos, *tr_data, *tr_pos, *ta_data, *ta_pos;
+  float dx, dy, r, a, asum, amean;
+  uint_t src_stride, dst_stride, tr_stride, ta_stride, width, height;
+  CvSize size;
+  int x, y;
+
+  width = src->width;
+  height = src->height;
+
+  size.width = width;
+  size.height = height;
+  tr = cvCreateImage(size, IPL_DEPTH_32F, 1);
+  ta = cvCreateImage(size, IPL_DEPTH_32F, 1);
+  dst = cvCreateImage(size, IPL_DEPTH_32F, 1);
+  cvSet(tr, cvScalarAll(0), NULL);
+  cvSet(ta, cvScalarAll(0), NULL);
+  cvSet(dst, cvScalarAll(0), NULL);
+
+  src_data = (float*)src->imageData;
+  src_stride = (uint_t)(src->widthStep / sizeof(float));
+  dst_data = (float*)dst->imageData;
+  dst_stride = (uint_t)(dst->widthStep / sizeof(float));
+  tr_data = (float*)tr->imageData;
+  tr_stride = (uint_t)(tr->widthStep / sizeof(float));
+  ta_data = (float*)ta->imageData;
+  ta_stride = (uint_t)(ta->widthStep / sizeof(float));
+
+  for (y = 2; y < height - 2; y++) {
+    src_pos = src_data + y * src_stride + 2;
+    tr_pos = tr_data + y * tr_stride + 2;
+    ta_pos = ta_data + y * ta_stride + 2;
+    for (x = 2; x < width - 2; x++, src_pos++, tr_pos++, ta_pos++) {
+      dx = calc_dx(src_pos);
+      dy = calc_dy(src_pos, src_stride);
+      r = sqrt(dx*dx+dy*dy);
+      if (dx == 0 && dy == 0) {
+        a = 0;
+      }
+      else {
+        a = atan2(dy,dx);
+      }
+      *tr_pos = r;
+      *ta_pos = a;
+    }
+  }
+
+  for (y = 2; y < height - 2; y++) {
+    dst_pos = dst_data + y * dst_stride + 2;
+    tr_pos = tr_data + y * tr_stride + 2;
+    ta_pos = ta_data + y * ta_stride + 2;
+    for (x = 2; x < width - 2; x++, dst_pos++, tr_pos++, ta_pos++) {
+      a = calc_directionality_1(ta_pos, ta_stride) + 
+          0.5 * calc_directionality_2(ta_pos, ta_stride) +
+          calc_directionality_3(ta_pos, ta_stride) +
+          0.5 * calc_directionality_4(ta_pos, ta_stride);
+      /*
+      asum = calc_sum9(ta_pos, ta_stride);
+      amean = asum / 9;
+      */
+      /*r = sqrt(dx*dx+dy*dy);*/
+      /*
+      if (abs(a - amean) > 0.3) {
+        *dst_pos = a;
+      }
+      else {
+        *dst_pos = 0;
+      }
+      */
+      *dst_pos = a;
+    }
+  }
+
+  cvReleaseImage(&tr);
+  cvReleaseImage(&ta);
+  return dst;
+}
+
+IplImage *mean_fill_holes(IplImage *src, int r, float ignore_val, float ignore_eps)
+{
+  IplImage *dst;
+  float *src_data, *src_pos, *dst_data, *dst_pos, value, sum;
+  uint_t src_stride, dst_stride, width, height, sum_count, min_count;
+  CvSize size;
+  int x, y, p, d;/*, mask_size, *mask_idx;*/
+
+  width = src->width;
+  height = src->height;
+
+  size.width = width;
+  size.height = height;
+  dst = cvCreateImage(size, IPL_DEPTH_32F, 1);
+  cvSet(dst, cvScalarAll(ignore_val), NULL);
+
+  src_data = (float*)src->imageData;
+  src_stride = (uint_t)(src->widthStep / sizeof(float));
+  dst_data = (float*)dst->imageData;
+  dst_stride = (uint_t)(dst->widthStep / sizeof(float));
+
+  d = (2*r+1);
+  /* create the mask offset table for the averaging filter mask */
+  mask_size = d*d;
+  mask_idx = malloc(mask_size * sizeof(int));
+  for (y = -r, p = 0; y <= r; y++) {
+    for (x = -r; x <= r; x++) {
+      mask_idx[p++] = y * src_stride + x;
+    }
+  }
+
+  min_count = (uint_t)(mask_size / 2);
+
+  for (y = r; y < height - r; y++) {
+    src_pos = src_data + y * src_stride + r;
+    dst_pos = dst_data + y * dst_stride + r;
+    for (x = r; x < width - r; x++, src_pos++, dst_pos++) {
+      value = *src_pos;
+      if ((ignore_val - ignore_eps) < value && value < (ignore_val + ignore_eps)) {
+        sum = calc_sum(src_pos, &sum_count);
+        if (sum_count > min_count) {
+          value = (sum/sum_count);
+        }
+      }
+      *dst_pos = value;
+    }
+  }
+
+  free(mask_idx);
+  mask_idx = NULL;
+  return dst;
+}
+
+IplImage *filter_mean(IplImage *src, int r, float ignore_val, float ignore_eps)
+{
+  IplImage *dst;
+  float *src_data, *src_pos, *dst_data, *dst_pos, value, sum;
+  uint_t src_stride, dst_stride, width, height, sum_count, min_count;
+  CvSize size;
+  int x, y, p, d;/*, mask_size, *mask_idx;*/
+
+  width = src->width;
+  height = src->height;
+
+  size.width = width;
+  size.height = height;
+  dst = cvCreateImage(size, IPL_DEPTH_32F, 1);
+  cvSet(dst, cvScalarAll(ignore_val), NULL);
+
+  src_data = (float*)src->imageData;
+  src_stride = (uint_t)(src->widthStep / sizeof(float));
+  dst_data = (float*)dst->imageData;
+  dst_stride = (uint_t)(dst->widthStep / sizeof(float));
+
+  d = (2*r+1);
+  /* create the mask offset table for the averaging filter mask */
+  mask_size = d*d;
+  mask_idx = malloc(mask_size * sizeof(int));
+  for (y = -r, p = 0; y <= r; y++) {
+    for (x = -r; x <= r; x++) {
+      mask_idx[p++] = y * src_stride + x;
+    }
+  }
+
+  min_count = (uint_t)(mask_size / 2);
+
+  for (y = r; y < height - r; y++) {
+    src_pos = src_data + y * src_stride + r;
+    dst_pos = dst_data + y * dst_stride + r;
+    for (x = r; x < width - r; x++, src_pos++, dst_pos++) {
+      value = *src_pos;
+      sum = calc_sum(src_pos, &sum_count);
+      if (sum_count > min_count) {
+        *dst_pos = (sum / sum_count);
+      }
+      else {
+        *dst_pos = ignore_val;
+      }
+    }
+  }
+
+  free(mask_idx);
+  mask_idx = NULL;
+  return dst;
+}
+
+IplImage *image_diff(IplImage *src1, IplImage *src2, float ignore_val, float ignore_eps)
+{
+  IplImage *dst;
+  float *src1_data, *src2_data, *dst_data, *src1_pos, *src2_pos, *dst_pos, value, diff;
+  uint_t src1_stride, src2_stride, dst_stride, width, height;
+  CvSize size;
+  int x, y;
+
+  width = src1->width;
+  height = src1->height;
+
+  size.width = width;
+  size.height = height;
+  dst = cvCreateImage(size, IPL_DEPTH_32F, 1);
+  cvSet(dst, cvScalarAll(ignore_val), NULL);
+
+  src1_data = (float*)src1->imageData;
+  src1_stride = (uint_t)(src1->widthStep / sizeof(float));
+  src2_data = (float*)src2->imageData;
+  src2_stride = (uint_t)(src2->widthStep / sizeof(float));
+  dst_data = (float*)dst->imageData;
+  dst_stride = (uint_t)(dst->widthStep / sizeof(float));
+
+  for (y = 0; y < height; y++) {
+    src1_pos = src1_data + y * src1_stride;
+    src2_pos = src2_data + y * src2_stride;
+    dst_pos = dst_data + y * dst_stride;
+    for (x = 0; x < width; x++, src1_pos++, src2_pos++, dst_pos++) {
+      value = *src1_pos;
+      if (value > ignore_val + ignore_eps || value < ignore_val - ignore_eps) {
+        diff = abs(value - *src2_pos);
+        *dst_pos = diff;
+      }
+      else {
+        *dst_pos = 0;
+      }
+    }
+  }
+
+  return dst;
+}
+
+IplImage *abs_diff_mean(IplImage *src, float ignore_val, float ignore_eps)
+{
+  IplImage *dst;
+  float *src_data, *dst_data, *src_pos, *dst_pos;
+  float value, sum, count, avg;
+  uint_t src_stride, dst_stride, width, height;
+  CvSize size;
+  int x, y;
+
+  width = src->width;
+  height = src->height;
+
+  size.width = width;
+  size.height = height;
+  dst = cvCreateImage(size, IPL_DEPTH_32F, 1);
+  cvSet(dst, cvScalarAll(0), NULL);
+
+  src_data = (float*)src->imageData;
+  src_stride = (uint_t)(src->widthStep / sizeof(float));
+  dst_data = (float*)dst->imageData;
+  dst_stride = (uint_t)(dst->widthStep / sizeof(float));
+
+  sum = 0;
+  count = 0;
+  for (y = 0; y < height; y++) {
+    src_pos = src_data + y * src_stride;
+    for (x = 0; x < width; x++, src_pos++) {
+      value = *src_pos;
+      if (value > (ignore_val + ignore_eps) || value < (ignore_val - ignore_eps)) {
+        sum += value;
+        count += 1;
+      }
+    }
+  }
+
+  avg = sum / count;
+
+  for (y = 0; y < height; y++) {
+    src_pos = src_data + y * src_stride;
+    dst_pos = dst_data + y * dst_stride;
+    for (x = 0; x < width; x++, src_pos++, dst_pos++) {
+      value = *src_pos;
+      if (value > ignore_val + ignore_eps || value < ignore_val - ignore_eps) {
+        *dst_pos = abs(value - avg);
+      }
+    }
+  }
+
+  return dst;
+}
+
+IplImage *box_mean_abs_diff(IplImage *src, int r, float ignore_val, float ignore_eps)
+{
+  IplImage *dst;
+  float *src_data, *src_pos, *dst_data, *dst_pos, value, sum;
+  uint_t src_stride, dst_stride, width, height, sum_count, min_count;
+  CvSize size;
+  int x, y, p, d;/*, mask_size, *mask_idx;*/
+
+  width = src->width;
+  height = src->height;
+
+  size.width = width;
+  size.height = height;
+  dst = cvCreateImage(size, IPL_DEPTH_32F, 1);
+  cvSet(dst, cvScalarAll(ignore_val), NULL);
+
+  src_data = (float*)src->imageData;
+  src_stride = (uint_t)(src->widthStep / sizeof(float));
+  dst_data = (float*)dst->imageData;
+  dst_stride = (uint_t)(dst->widthStep / sizeof(float));
+
+  d = (2*r+1);
+  /* create the mask offset table for the averaging filter mask */
+  mask_size = d*d;
+  mask_idx = malloc(mask_size * sizeof(int));
+  for (y = -r, p = 0; y <= r; y++) {
+    for (x = -r; x <= r; x++) {
+      mask_idx[p++] = y * src_stride + x;
+    }
+  }
+
+  min_count = (uint_t)(mask_size / 2);
+
+  for (y = r; y < height - r; y++) {
+    src_pos = src_data + y * src_stride + r;
+    dst_pos = dst_data + y * dst_stride + r;
+    for (x = r; x < width - r; x++, src_pos++, dst_pos++) {
+      value = *src_pos;
+      sum = calc_sum(src_pos, &sum_count);
+      if (sum_count > min_count) {
+        if (value > ignore_val + ignore_eps || value < ignore_val - ignore_eps) {
+          value = abs((sum/sum_count) - value);
+        }
+        else {
+          value = ignore_val;
+        }
+        *dst_pos = value;
+      }
+      else {
+        *dst_pos = ignore_val;
+      }
+    }
+  }
+
+  free(mask_idx);
+  mask_idx = NULL;
+  return dst;
+}
+
+void stretch_histogram_avg_sdv(IplImage *img, float ignore_val, float ignore_eps)
+{
+  float *img_data, *img_pos;
+  uint_t img_stride, x, y, width, height;
+  float value, min, max, sum1, sum2, count, avg, sdv;
+  float bound_1, bound_2, stretch_1, stretch_2, stretch_3, stretch_4;
+
+  width = img->width;
+  height = img->height;
+  img_data = (float*)img->imageData;
+  img_stride = (uint_t)(img->widthStep / sizeof(float));
+
+  min = 2000000000;
+  max = 0;
+  count = 0;
+  sum1 = 0;
+  sum2 = 0;
+  for (y = 0; y < height; y++) {
+    img_pos = img_data + y * img_stride;
+    for (x = 0; x < width; x++, img_pos++) {
+      value = *img_pos;
+      if (value > (ignore_val + ignore_eps) || value < (ignore_val - ignore_eps)) {
+        if (value < min) min = value; else if (value > max) max = value;
+        sum1 += value;
+        sum2 += value*value;
+        count += 1;
+      }
+    }
+  }
+
+  avg = sum1 / count;
+  sdv = sqrt(sum2 / count - avg*avg);
+  printf("min=%f max=%f avg=%f sdv=%f\n", min, max, avg, sdv);
+
+  if ((avg - 2 * sdv) > min) {
+    min = avg - 2 * sdv;
+  }
+  if ((avg + 2 * sdv) < max) {
+    max = avg + 2 * sdv;
+  }
+  if (avg - sdv < min) {
+    bound_1 = min + 0.32 * (avg-min);
+  }
+  else {
+    bound_1 = avg - sdv;
+  }
+  if (avg + sdv > max) {
+    bound_2 = max - 0.32 * (max-avg);
+  }
+  else {
+    bound_2 = avg + sdv;
+  }
+  stretch_1 = 0.25 / (bound_1 - min);
+  stretch_2 = 0.25 / (avg - bound_1);
+  stretch_3 = 0.25 / (bound_2 - avg);
+  stretch_4 = 0.25 / (max - bound_2);
+
+  /* normalize values and create the histogram */
+  printf("b1=%f b2=%f b3=%f b4=%f b5=%f s1=%f s2=%f s3=%f s4=%f\n", min, bound_1, avg, bound_2, max, stretch_1, stretch_2, stretch_3, stretch_4);
+
+  for (y = 0; y < height; y++) {
+    img_pos = img_data + y * img_stride;
+    for (x = 0; x < width; x++, img_pos++) {
+      value = *img_pos;
+      if ((ignore_val - ignore_eps) < value && value < (ignore_val + ignore_eps)) {
+        value = 0;
+      }
+      else
+      if (value < min) {
+        value = 0;
+      }
+      else
+      if (value < bound_1) {
+        value = stretch_1 * (value - min);
+      }
+      else
+      if (value < avg) {
+        value = 0.25 + stretch_2 * (value - bound_1);
+      }
+      else
+      if (value < bound_2) {
+        value = 0.5 + stretch_3 * (value - avg);
+      }
+      else
+      if (value < max) {
+        value = 0.75 + stretch_4 * (value - bound_2);
+      }
+      else {
+        value = 1;
+      }
+
+      *img_pos = value;
+    }
+  }
+}
+
+void equalize_histogram(IplImage *img, float ignore_val, float ignore_eps)
+{
+  float *img_data, *img_pos;
+  uint_t img_stride, i, x, y, width, height, hidx;
+  float value, count, min, max, stretch;
+  float hist[256];
+
+  width = img->width;
+  height = img->height;
+  img_data = (float*)img->imageData;
+  img_stride = (uint_t)(img->widthStep / sizeof(float));
+
+  min = 2000000000;
+  max = 0;
+  for (y = 0; y < height; y++) {
+    img_pos = img_data + y * img_stride;
+    for (x = 0; x < width; x++, img_pos++) {
+      value = *img_pos;
+      if (value > (ignore_val + ignore_eps) || value < (ignore_val - ignore_eps)) {
+        if (value < min) min = value; else if (value > max) max = value;
+      }
+    }
+  }
+
+  stretch = 1.0 / (max - min);
+  printf("min=%f max=%f stretch=%f\n", min, max, stretch);
+
+  /* normalize values and create the histogram */
+  for (i = 0; i < 256; i++) {
+    hist[i] = 0;
+  }
+
+  count = 0;
+  for (y = 0; y < height; y++) {
+    img_pos = img_data + y * img_stride;
+    for (x = 0; x < width; x++, img_pos++) {
+      value = *img_pos;
+      if (value > ignore_val + ignore_eps || value < ignore_val - ignore_eps) {
+        value = stretch * (value - min);
+        if (value < 0) value = 0; else if (value > 1) value = 1;
+        hidx = (int)(255 * value);
+        if (hidx < 0) hidx = 0; else if (hidx > 255) hidx = 255;
+        hist[hidx] += 1;
+        count += 1;
+      }
+      else {
+        value = 0;
+      }
+      *img_pos = value;
+    }
+  }
+
+  /* transform the histogram into a cumulative distribution */
+  printf("%f ", hist[0]);
+  hist[0] /= count;
+  for (i = 1; i < 256; i++) {
+    /*printf("(%.2f,", hist[i]);*/
+    hist[i] = (hist[i] / count) + hist[i-1];
+    /*printf("%.2f) ", hist[i]);*/
+  }
+  /*printf("\n");*/
+
+  for (y = 0; y < height; y++) {
+    img_pos = img_data + y * img_stride;
+    for (x = 0; x < width; x++, img_pos++) {
+      value = *img_pos;
+      hidx = (int)(255 * value);
+      if (hidx < 0) hidx = 0; else if (hidx > 255) hidx = 255;
+      value = (hist[hidx] * value);
+      if (value < 0) value = 0; else if (value > 1) value = 1;
+      *img_pos = value;
+    }
+  }
+}
+
+void stretch_histogram(IplImage *img, float ignore_val, float ignore_eps)
+{
+  float *img_data, *img_pos;
+  uint_t img_stride, x, y, width, height, count;
+  float value, min, max, stretch;
+
+  width = img->width;
+  height = img->height;
+  img_data = (float*)img->imageData;
+  img_stride = (uint_t)(img->widthStep / sizeof(float));
+
+  min = 2000000000;
+  max = 0;
+  for (y = 0; y < height; y++) {
+    img_pos = img_data + y * img_stride;
+    for (x = 0; x < width; x++, img_pos++) {
+      value = *img_pos;
+      if (value > (ignore_val + ignore_eps) || value < (ignore_val - ignore_eps)) {
+        if (value < min) min = value; else if (value > max) max = value;
+      }
+    }
+  }
+
+  stretch = 1.0 / (max - min);
+  printf("min=%f max=%f stretch=%f\n", min, max, stretch);
+
+  for (y = 0; y < height; y++) {
+    img_pos = img_data + y * img_stride;
+    for (x = 0; x < width; x++, img_pos++) {
+      value = *img_pos;
+      if (value > (ignore_val + ignore_eps) || value < (ignore_val - ignore_eps)) {
+        value = stretch * (value - min);
+        if (value < 0) value = 0; else if (value > 1) value = 1;
+      }
+      else {
+        value = 0;
+      }
+      *img_pos = value;
+    }
+  }
+}
+
+IplImage *to_8bit(IplImage *src)
+{
+  IplImage *dst;
+  float *src_data, *src_pos, value;
+  unsigned char *dst_data, *dst_pos;
+  uint_t x, y, width, height, src_stride, dst_stride;
+  int temp;
+  CvSize size;
+
+  width = src->width;
+  height = src->height;
+
+  size.width = width;
+  size.height = height;
+  dst = cvCreateImage(size, IPL_DEPTH_8U, 1);
+
+  src_data = (float*)src->imageData;
+  src_stride = (uint_t)(src->widthStep / sizeof(float));
+  dst_data = (unsigned char*)dst->imageData;
+  dst_stride = (uint_t)(dst->widthStep / sizeof(unsigned char));
+
+  for (y = 0; y < height; y++) {
+    src_pos = src_data + y * src_stride;
+    dst_pos = dst_data + y * dst_stride;
+    for (x = 0; x < width; x++, src_pos++, dst_pos++) {
+      value = *src_pos;
+      temp = (int)(255 * value);
+      if (temp < 0) temp = 0;
+      else if (temp > 255) temp = 255;
+      *dst_pos = (unsigned char)temp;
+    }
+  }
+
+  return dst;
+}
+
+IplImage *read_from_tcr_mag(const char *path)
+{
+  IplImage *src, *tmp1, *tmp2, *dst;
+  uint_t *timestamps;
+
+  printf("read\n");
+  src = read_from_tcr(path, &timestamps);
+  printf("filter\n");
+  tmp1 = filter_mean(src, 3, 0, 0.001);
+  cvReleaseImage(&src);
+  printf("derivate\n");
+  tmp2 = derivative_direction(tmp1, 0, 0.001);
+  cvReleaseImage(&tmp1);
+  printf("equalize\n");
+  equalize_histogram(tmp2, 0, 0.001);
+  printf("convert\n");
+  dst = to_8bit(tmp2);
+  cvReleaseImage(&tmp2);
+  cvSaveImage("tcr_mag.png", dst, 0);
+  return dst;
+}
+
+IplImage *calc_derivative_direction(IplImage *src)
+{
+  IplImage *dst;
+  dst = derivative_direction(src, 0, 0.001);
+  equalize_histogram(dst, 0, 0.001);
+  return dst;
+}
+
+IplImage *read_from_tcr_rectified(const char *path)
+{
+  IplImage *src, *tmp1, *tmp2, *dst;
+  float *tmp_data, *tmp_pos, value;
+  uint_t *timestamps, tmp_stride, width, height, x, y, r;
+  uint_t minx, maxx, startx, endx, miny, maxy, starty, endy;
+
+  src = read_from_tcr(path, &timestamps);
+  if (src == NULL || timestamps == NULL) {
+    return NULL;
+  }
+
+  /*src = mean_fill_holes(tmp, 3, 0, 0.001);*/
+  /*cvReleaseImage(&tmp);*/
+  /*tmp1 = filter_mean(src, 5, 0, 0.001);*/
+  /*cvSaveImage("tcr_mean.png", tmp, 0);*/
+  /*tmp2 = image_diff(src, tmp1, 0, 0.001);*/
+  /*tmp2 = abs_diff_mean(tmp1, 0, 0.001);*/
+  /*cvSaveImage("tcr_diff.png", dst, 0);*/
+  tmp1 = box_mean_abs_diff(src, 5, 0, 0.001);
+  cvReleaseImage(&src);
+  /*cvReleaseImage(&tmp1);*/
+  stretch_histogram_avg_sdv(tmp1, 0, 0.001);
+  /*equalize_histogram(tmp2, 0, 0.001);*/
+  /*stretch_histogram(tmp2, 0, 0.001);*/
+
+  width = tmp1->width;
+  height = tmp1->height;
+  tmp_data = (float*)tmp1->imageData;
+  tmp_stride = (uint_t)(tmp1->widthStep / sizeof(float));
+
+  minx = 2000000000;
+  maxx = 0;
+  miny = 2000000000;
+  maxy = 0;
+  for (y = 0; y < height; y++) {
+    tmp_pos = tmp_data + y * tmp_stride;
+    for (x = 0; x < width; x++, tmp_pos++) {
+      value = *tmp_pos;
+      if (value > 0.001) {
+        if (x < minx) minx = x; else if (x > maxx) maxx = x;
+        if (y < miny) miny = y; else if (y > maxy) maxy = y;
+      }
+    }
+  }
+  printf("minx=%d maxx=%d miny=%d maxy=%d\n", minx, maxx, miny, maxy);
+
+  startx = 0;
+  if (minx > 2) startx = minx - 2;
+  endx = width - 1;
+  if (maxx < endx - 2) endx = maxx + 2;
+  starty = 0;
+  if (miny >  2) starty = miny - 2;
+  endy = height - 1;
+  if (maxy < endy - 2) endy = maxy + 2;
+
+  tmp2 = rectify_tcr(tmp1, timestamps, startx, endx, starty, endy);
+  cvReleaseImage(&tmp1);
+  dst = to_8bit(tmp2);
+  cvReleaseImage(&tmp2);
+  cvSaveImage("tcr_eq.png", dst, 0);
+  free(timestamps);
+
+  return dst;
+}
+
+IplImage *read_from_tcr_rectified_old(const char *path)
+{
+  IplImage *src, *tmp, *tmp2, *tmp3, *tmp4, *dst;
+  float *src_data, *tmp_data, *dst_data, *src_pos, *tmp_pos, *dst_pos;
+  int src_stride, tmp_stride, dst_stride;
+  unsigned int *timestamps, *timestamps_accum;
+  unsigned int timestamp, timestamp_count, tstart, tend, tdiff;
+  float tratio;
+  CvSize size;
+  size_t width, height, read_size;
+  int line_width, timestamp_width, line_count, file_size;
+  FILE *input_file;
+  /* histogram for grayscale equalization */
+  float hist[256];
+
+  width = 2592;
+  line_width = width * sizeof(float);
+  timestamp_width = sizeof(unsigned int);
+  src = NULL;
+  tmp = NULL;
+  dst = NULL;
+
+  input_file = fopen(path, "rb");
+  if (input_file != NULL) {
+    /* seek file end to determine image size */
+    fseek (input_file , 0 , SEEK_END);
+    file_size = ftell(input_file);
+    rewind (input_file);
+    height = (size_t)(file_size / (line_width + timestamp_width));
+    /*printf("File size: %d, rows: %d\n", file_size, height);*/
+
+    size.width = width;
+    size.height = height;
+    src = cvCreateImage(size, IPL_DEPTH_32F, 1);
+    tmp = cvCreateImage(size, IPL_DEPTH_32F, 1);
+    tmp2 = cvCreateImage(size, IPL_DEPTH_32F, 1);
+    tmp3 = cvCreateImage(size, IPL_DEPTH_32F, 1);
+    tmp4 = cvCreateImage(size, IPL_DEPTH_32F, 1);
+    src_data = (float*)src->imageData;
+    tmp_data = (float*)tmp->imageData;
+    src_stride = (int)(src->widthStep / sizeof(float));
+    tmp_stride = (int)(tmp->widthStep / sizeof(float));
+
+    cvSet(tmp, cvScalarAll(0), NULL);
+
+    uint_t r = 2;
+    uint_t d = (2*r+1);
+    /* create the mask offset table for the averaging filter mask */
+    mask_size = d*d;
+    mask_idx = malloc(mask_size * sizeof(int));
+    for (int i = -r, pos = 0; i <= r; i++) {
+      for (int j = -r; j <= r; j++) {
+        mask_idx[pos++] = i * src_stride + j;
+        /*printf("%d ", mask_idx[pos-1]);*/
+      }
+    }
+    /*printf("\n");*/
+
+    timestamps = (unsigned int *)malloc(height * timestamp_width);
+    if (src_data != NULL) {
+      line_count = 0;
+      timestamp_count = 0;
+
+      src_pos = src_data;
+      while (1) {
+        /* read line timestamp to use for rectification */
+        read_size = fread((void*)&timestamp, timestamp_width, 1, input_file);
+        if (read_size < 1) {
+          break;
+        }
+        else {
+          timestamps[timestamp_count++] = timestamp;
+        }
+
+        /* read data line */
+        read_size = fread((void*)src_pos, sizeof(float), width, input_file);
+        if (read_size < width) {
+          break;
+        }
+        else {
+          line_count += 1;
+          src_pos += src_stride;
+        }
+      }
+      fclose(input_file);
+
+      /*cvSaveImage("tcr_raw.png", src, 0);*/
+      /*
+      cvSmooth(src,tmp,CV_GAUSSIAN,5,5,0,0);
+      cvLaplace(tmp,tmp2,5);
+      */
+      /*
+      cvSmooth(src,tmp,CV_BLUR,5,5,0,0);
+      cvSmooth(src,tmp2,CV_GAUSSIAN,3,3,0,0);
+      cvAbsDiff(tmp,src,tmp3);
+      cvSaveImage("tcr_diff.png", tmp3, 0);
+      cvSaveImage("tcr_smooth2.png", tmp2, 0);
+      */
+      /*
+      cvSaveImage("tcr_laplacian5.png", tmp2, 0);
+      */
+      {
+        unsigned int x, y, minx, maxx, miny, maxy, count;
+        unsigned int new_width, new_height, startx, endx, starty, endy;
+        float value, sum, sum1, sum2, min, max, ratio, stretch, avg, sdv;
+        float sum_count, min_count, stretch_1, stretch_2;
+
+        min = 2000000000;
+        max = 0;
+        minx = 2000000000;
+        maxx = 0;
+        miny = 2000000000;
+        maxy = 0;
+        count = 0;
+        sum1 = 0;
+        sum2 = 0;
+        min_count = mask_size / 2.0;
+        for (y = r; y < height - r; y++) {
+          src_pos = src_data + y * src_stride + r;
+          tmp_pos = tmp_data + y * tmp_stride + r;
+          for (x = r; x < width - r; x++, src_pos++, tmp_pos++) {
+            value = *src_pos;
+            sum = calc_sum(src_pos, &sum_count);
+            if (sum_count > min_count) {
+              if (value > epsilon) {
+                value = abs((sum/sum_count) - value);
+                if (value < min) min = value; else if (value > max) max = value;
+                if (x < minx) minx = x; else if (x > maxx) maxx = x;
+                if (y < miny) miny = y; else if (y > maxy) maxy = y;
+                sum1 += value;
+                sum2 += value*value;
+                count += 1;
+              }
+              else {
+                value = 0;//sum/count;
+              }
+              *tmp_pos = value;
+            }
+            else {
+              *tmp_pos = 0;
+            }
+          }
+        }
+        printf("\n");
+        avg = sum1 / count;
+        sdv = sqrt(sum2 / count - avg*avg);
+        printf("min=%f max=%f avg=%f sdv=%f minx=%d maxx=%d miny=%d maxy=%d\n", min, max, avg, sdv, minx, maxx, miny, maxy);
+        /*cvSaveImage("tcr_filtered.png", tmp, 0);*/
+
+        startx = r;
+        if (minx > startx + 2) startx = minx - 2;
+        endx = width - r - 1;
+        if (maxx < endx - 2) endx = maxx + 2;
+        starty = r;
+        if (miny > starty + 2) starty = miny - 2;
+        endy = height - r - 1;
+        if (maxy < endy - 2) endy = maxy + 2;
+
+        new_width = endx - startx;
+        /*new_height = endy - starty;*/
+        new_height = new_width;
+
+        tstart = timestamps[starty];
+        tend = timestamps[endy];
+        tdiff = tend - tstart;
+        ratio = (float)tdiff / (float)new_height;
+        printf("start=%d end=%d diff=%d ratio=%f\n", tstart, tend, tdiff, ratio);
+
+        /* create a lookup table for rectified image time stamps */
+        timestamps_accum = (unsigned int *)malloc(new_height * timestamp_width);
+        timestamps_accum[0] = tstart;
+        for (int i = 1; i < new_height; i++) {
+          timestamps_accum[i] = tstart + (unsigned int)(ratio * i);
+        }
+
+        size.width = new_width;
+        size.height = new_height;
+        dst = cvCreateImage(size, IPL_DEPTH_32F, 1);
+        dst_data = (float*)dst->imageData;
+        dst_stride = (int)(dst->widthStep / sizeof(float));
+
+        dst_pos = dst_data;
+        int j = starty;
+        int row_size = new_width * sizeof(float);
+        tmp_pos = tmp_data + starty * width + startx;
+        for (int i = 0; i < new_height; i++, dst_pos += dst_stride) {
+          if (timestamps[j] < timestamps_accum[i]) {
+            j += 1;
+            if (j == height) {
+              printf("row overflow\n");
+              break;
+            }
+            tmp_pos += tmp_stride;
+          }
+          memcpy(dst_pos, tmp_pos, row_size);
+        }
+
+        free(timestamps_accum);
+        /*cvSaveImage("tcr_rectified.png", dst, 0);*/
+
+        int hidx;
+        /* create a cumulative gaussian distribution table for normalizing the histogram */
+        float norm[256];
+        float nsum = 0;
+        for (int i = 0; i < 256; i++) {
+          norm[i] = exp(-(pow(((float)i/255.0) - 0.5, 2.0) / (2*0.25*0.25)));
+          nsum += norm[i];
+        }
+        norm[0] /= nsum;
+        for (int i = 1; i < 256; i++) {
+          norm[i] = (norm[i] / nsum) + norm[i-1];
+          /*printf("%f ", norm[i]);*/
+        }
+        /*printf("\n");*/
+        for (int i = 0; i < 256; i++) {
+          hist[i] = 0;
+        }
+        stretch = 1.0 / ((avg+2*sdv) - 0);
+        stretch_1 = 0.5 / (avg - 0);
+        stretch_2 = 0.5 / (max - avg);
+
+        /* normalize values and create the histogram */
+        printf("stretch_1=%f stretch_2=%f\n", stretch_1, stretch_2);
+        count = 0;
+        for (y = 0; y < new_height; y++) {
+          dst_pos = dst_data + y * dst_stride;
+          for (x = 0; x < new_width; x++, dst_pos++) {
+            value = *dst_pos;
+            if (value > epsilon) {
+              if (value < (avg+epsilon)) {
+                value = stretch_1 * (value - 0);
+              }
+              else {
+                value = 0.5 + stretch_2 * (value - avg);
+              }
+              /*value = stretch * (value - 0);*/
+              hidx = (int)(255.0 * value);
+              if (hidx < 1) hidx = 1;
+              else if (hidx > 255) hidx = 255;
+              hist[hidx] += 1;
+              count += 1;
+              *dst_pos = value;
+            }
+            else {
+              *dst_pos = 0;
+            }
+          }
+        }
+        /*cvSaveImage("tcr_normalized.png", dst, 0);*/
+
+        for (int i = 1; i < 256; i++) {
+          printf("%.3f ", hist[i]);
+        }
+        printf("\n");
+        /* transform the histogram into a cumulative distribution */
+        hist[0] = 0;
+        for (int i = 1; i < 256; i++) {
+          hist[i] = (hist[i] / count) + hist[i-1];
+          printf("%.3f ", hist[i]);
+        }
+        printf("\n");
+
+        for (y = 0; y < new_height; y++) {
+          dst_pos = dst_data + y * dst_stride;
+          for (x = 0; x < new_width; x++, dst_pos++) {
+            value = *dst_pos;
+            if (value < epsilon) {
+              value = 0;
+            }
+            else {
+              hidx = (int)(255.0 * value);
+              if (hidx < 1) hidx = 1;
+              else if (hidx > 255) hidx = 255;
+              value = (hist[hidx] * value); /*norm[hidx] * */
+            }
+            *dst_pos = value;
+          }
+        }
+        /*cvSaveImage("tcr_stretched.png", dst, 0);*/
+      }
+    }
+    else {
+      printf("Failed to allocate memory for image");
+    }
+    free(mask_idx);
+    free(timestamps);
+    cvReleaseImage(&tmp);
+    cvReleaseImage(&tmp2);
+    cvReleaseImage(&tmp3);
+    cvReleaseImage(&tmp4);
+    cvReleaseImage(&src);
+  }
+  else {
+    printf("Failed to open file %s\n", path);
+  }
+
+  return dst;
+}
+
+IplImage *read_from_tcr_projected(const char *path)
+{
+
+}
+
+IplImage *read_from_tcr_old(const char *path)
+{
+  IplImage *tmp, *dst;
+  CvSize size;
+  float *tmp_data, *dst_data, *tmp_pos, *dst_pos;
+  /*unsigned char *b_data, *b_pos;*/
+  float hist[256];
+  size_t width = 2592;
+  size_t height;
+  int tmp_stride, dst_stride;/*, b_stride;*/
+  int line_width = width * sizeof(float);
+  int timestamp_width = sizeof(unsigned int);
+
+  tmp = NULL;
+  dst = NULL;
+  FILE *input_file = fopen(path, "rb");
+  size_t readsize;
+  int line_count;
+  int file_size;
+  if (input_file != NULL) {
+    unsigned int *timestamps, *timestamps_accum;
+    unsigned int timestamp, prev_timestamp, timestamp_diff;
+    unsigned int total_timestamp_diff;
+    unsigned int total_timestamp_diff_2;
+    unsigned int timestamp_count;
+    float timestamp_avg, timestamp_sdv;
+
+    fseek (input_file , 0 , SEEK_END);
+    file_size = ftell (input_file);
+    rewind (input_file);
+    height = (size_t)(file_size / (line_width + timestamp_width));
+    printf("File size: %d, rows: %d\n", file_size, height);
+    size.width = width;
+    size.height = height;
+    tmp = cvCreateImage(size, IPL_DEPTH_32F, 1);
+    tmp_data = (float*)tmp->imageData;
+    tmp_stride = (int)(tmp->widthStep / sizeof(float));
+    timestamps = (unsigned int *)malloc(height * timestamp_width);
+    timestamps_accum = (unsigned int *)malloc(height * timestamp_width);
+    if (tmp_data != NULL) {
+      line_count = 0;
+      prev_timestamp = 0;
+      total_timestamp_diff = 0;
+      total_timestamp_diff_2 = 0;
+      timestamp_count = 1;
+      timestamps[0] = 0;
+      tmp_pos = tmp_data;
+
+      while (1) {
+        readsize = fread((void*)&timestamp, timestamp_width, 1, input_file);
+        /*printf("Read %d elements for timestamp\n", readsize);*/
+        if (readsize < 1) {
+          break;
+        }
+        else {
+          if (prev_timestamp > 0) {
+            timestamp_diff = (timestamp - prev_timestamp);
+            /*printf("%d ", timestamp_diff);*/
+            total_timestamp_diff += timestamp_diff;
+            total_timestamp_diff_2 += (timestamp_diff * timestamp_diff);
+            timestamps[timestamp_count] = timestamp_diff;
+            timestamp_count += 1;
+          }
+          prev_timestamp = timestamp;
+        }
+        readsize = fread((void*)tmp_pos, sizeof(float), width, input_file);
+        /*printf("Read %d elements for line\n", readsize);*/
+        if (readsize < width) {
+          break;
+        }
+        else {
+          line_count += 1;
+          tmp_pos += tmp_stride;
+        }
+      }
+
+      fclose(input_file);
+      timestamp_avg = total_timestamp_diff / timestamp_count;
+      timestamp_sdv = sqrt(total_timestamp_diff_2 / timestamp_count - timestamp_avg*timestamp_avg);
+      printf("Read %d lines, avg diff %f, diff sdv %f\n", line_count, timestamp_avg, timestamp_sdv);
+
+      timestamps_accum[0] = 0;
+      for (int i = i; i < height; i++) {
+        timestamps_accum[i] = timestamps[i] / total_timestamp_diff + timestamps_accum[i-1];
+      }
+
+      {
+        unsigned int pos, x, y, px, py, minx, maxx, miny, maxy, sumx, sumy, count;
+        unsigned int new_width, new_height;
+        float value, min, max, sum, sum2, avg, sdv, ratio, stretch, stretch_1, stretch_2;
+        float a, r, rmin, rmax, cx, cy;
+        unsigned char bval;
+
+        min = 2000000000;
+        max = 0;
+        sum = 0;
+        sum2 = 0;
+        minx = 2000000000;
+        maxx = 0;
+        miny = 2000000000;
+        maxy = 0;
+        sumx = 0;
+        sumy = 0;
+        count = 0;
+
+        for (y = 0; y < height; y++) {
+          tmp_pos = tmp_data + y * tmp_stride;
+          for (x = 0; x < width; x++, tmp_pos++) {
+            value = *tmp_pos;
+            if (value > 0.0001) {
+              if (value < min) min = value;
+              if (value > max) max = value;
+              if (x < minx) minx = x;
+              if (x > maxx) maxx = x;
+              if (y < miny) miny = y;
+              if (y > maxy) maxy = y;
+              sum += value;
+              sum2 += (value * value);
+              sumx += x;
+              sumy += y;
+              count += 1;
+            }
+          }
+        }
+
+        for (int i = 0; i < 256; i++) {
+          hist[i] = 0;
+        }
+        int hidx;
+
+        avg = sum / count;
+        sdv = sqrt((sum2 / count) - (avg * avg));
+        cx = (float)sumx / count;
+        cy = (float)sumy / count;
+        ratio = (float)(maxy - miny) / (float)(maxx - minx);
+
+        printf("amin=%f amax=%f aavg=%f asdv=%f\n", min, max, avg, sdv);
+        printf("cx=%f cy=%f ratio=%f\n", cx, cy, ratio);
+
+        stretch = 1.0 / (max - min);
+        stretch_1 = 0.2 / (avg - min);
+        stretch_2 = 0.8 / (max - avg);
+
+        rmin = 2000000000;
+        rmax = 0;
+        count = 0;
+        for (y = 0; y < height; y++) {
+          tmp_pos = tmp_data + y * tmp_stride;
+          for (x = 0; x < width; x++, tmp_pos++) {
+            value = *tmp_pos;
+            if (value > 0.0001) {
+              px = floor(x - cx);
+              py = floor((y / ratio) - (cy / ratio));
+              r = sqrt(py*py + px*px);
+              if (r < rmin) rmin = r;
+              if (r > rmax) rmax = r;
+              hidx = (int)(255 * stretch * (value - min));
+              if (hidx < 1) hidx = 1;
+              else if (hidx > 255) hidx = 255;
+              hist[hidx] += 1;
+              count += 1;
+            }
+          }
+        }
+
+        hist[0] = 0;
+        for (int i = 1; i < 256; i++) {
+          hist[i] = (hist[i] / count) + hist[i-1];
+          printf("%f ", hist[i]);
+        }
+
+        new_height = ceil(rmax - rmin);
+        new_width = ceil(2 * M_PI * rmax);
+        printf("w=%d h=%d rmin=%f rmax=%f\n", new_width, new_height, rmin, rmax);
+        size.width = new_width;
+        size.height = new_height;
+
+        dst = cvCreateImage(size, IPL_DEPTH_32F, 1);
+        dst_data = (float*)dst->imageData;
+        dst_stride = (int)(dst->widthStep / sizeof(float));
+
+        for (y = 0; y < new_height; y++) {
+          dst_pos = dst_data + y * dst_stride;
+          r = (float)(rmax - y);
+          for (x = 0; x < new_width; x++, dst_pos++) {
+            a = ((float)x / (float)new_width) * 2 * M_PI;
+            px = floor(cx + (r * cos(a)));
+            py = floor(((cy / ratio) + (r * sin(a))) * ratio);
+            if (px < 0 || py < 0 || px >= width || py >= height) {
+              /*printf("x=%d y=%d\n", px, py);*/
+            }
+            else {
+              value = *(tmp_data + py * tmp_stride + px);
+              if (value < 0.0001) {
+                value = 0;
+              }
+              if (value < avg) {
+                value = stretch_1 * (value - min);
+                hidx = (int)(255 * value);
+                if (hidx < 1) hidx = 1;
+                else if (hidx > 255) hidx = 255;
+                value = (hist[hidx] * value);
+              }
+              else {
+                value = 0.2 + stretch_2 * (value - avg);
+                hidx = (int)(255 * value);
+                if (hidx < 1) hidx = 1;
+                else if (hidx > 255) hidx = 255;
+                value = (hist[hidx] * value);
+              }
+              *dst_pos = value;
+            }
+          }
+        }
+      }
+    }
+    else {
+      printf("Failed to allocate memory for image");
+    }
+    cvReleaseImage(&tmp);
+  }
+  else {
+    printf("Failed to open file %s\n", path);
+  }
+  return dst;
+}
+
+
+IplImage *read_from_tcr_sdv(const char *path)
+{
+  IplImage *dst;
+  CvSize size;
+  float *data, *data2, *data_pos, *data2_pos;
+  size_t width = 2592;
+  size_t height;
+  int line_width = width * sizeof(float);
+  int timestamp_width = sizeof(unsigned int);
+
+  dst = NULL;
+  FILE *input_file = fopen(path, "rb");
+  size_t readsize;
+  int line_count;
+  int file_size;
+  if (input_file != NULL) {
+    unsigned int timestamp, prev_timestamp, timestamp_diff;
+    unsigned int total_timestamp_diff;
+    unsigned int total_timestamp_diff_2;
+    unsigned int timestamp_count;
+
+    fseek (input_file , 0 , SEEK_END);
+    file_size = ftell (input_file);
+    rewind (input_file);
+    height = (size_t)(file_size / (line_width + timestamp_width));
+    printf("File size: %d, rows: %d\n", file_size, height);
+
+    data = (float *)malloc(height * width * sizeof(float));
+    data2 = (float *)malloc(height * width * sizeof(float));
+    memset(data2, 0, height * width * sizeof(float));
+    if (data != NULL) {
+      line_count = 0;
+      prev_timestamp = 0;
+      total_timestamp_diff = 0;
+      total_timestamp_diff_2 = 0;
+      timestamp_count = 0;
+      data_pos = data;
+
+      while (1) {
+        readsize = fread((void*)&timestamp, timestamp_width, 1, input_file);
+        /*printf("Read %d elements for timestamp\n", readsize);*/
+        if (readsize < 1) {
+          break;
+        }
+        else {
+          if (prev_timestamp > 0) {
+            timestamp_diff = (timestamp - prev_timestamp);
+            total_timestamp_diff += timestamp_diff;
+            total_timestamp_diff_2 += (timestamp_diff * timestamp_diff);
+            timestamp_count += 1;
+          }
+          prev_timestamp = timestamp;
+        }
+        readsize = fread((void*)data_pos, sizeof(float), width, input_file);
+        /*printf("Read %d elements for line\n", readsize);*/
+        if (readsize < width) {
+          break;
+        }
+        else {
+          line_count += 1;
+          data_pos += width;
+        }
+      }
+      size.width = width;
+      size.height = line_count;
+      dst = cvCreateImageHeader(size, IPL_DEPTH_32F, 1);
+      cvSetData(dst, (void*)data2, line_width);
+      fclose(input_file);
+      printf("Read %d lines, average diff %d\n", line_count, (int)(total_timestamp_diff / timestamp_count));
+
+      {
+        int size, pos, x, y, count, c, c2;
+        float amin, amax, smin, smax, asum, asum2, ssum, ssum2, aavg, asdv, savg, ssdv;
+        float value,v1,v2,v3,v4,v5,v6,v7,v8,v9,vsum,vsum2,vavg,vavg2,vsdv,vsdv2;
+        float stretch_1, stretch_2;
+        amin = 2000000000;
+        smin = 2000000000;
+        amax = 0;
+        smax = 0;
+        asum = 0;
+        asum2 = 0;
+        ssum = 0;
+        ssum2 = 0;
+        count = 0;
+        vsum = 0;
+        vsum2 = 0;
+        vavg = 0;
+        vavg2 = 0;
+        vsdv = 0;
+        vsdv2 = 0;
+        for (y = 2; y < height - 2; y++) {
+          for (x = 2; x < width - 2; x++) {
+            data_pos = &data[y * width + x];
+            if (*data_pos > 0.0001) {
+              /*
+              v1 = *(data_pos - width - 1);
+              v2 = *(data_pos - width);
+              v3 = *(data_pos - width + 1);
+              v4 = *(data_pos - 1);
+              v5 = *data_pos;
+              v6 = *(data_pos + 1);
+              v7 = *(data_pos + width - 1);
+              v8 = *(data_pos + width);
+              v9 = *(data_pos + width + 1);
+              */
+              vsum = calc_sum(data_pos, &c);
+              /*v1+v2+v3+v4+v5+v6+v7+v8+v9;*/
+              vsum2 = calc_sum2(data_pos, &c2);
+              /*(v1*v1)+(v2*v2)+(v3*v3)+(v4*v4)+(v5*v5)+(v6*v6)+(v7*v7)+(v8*v8)+(v9*v9);*/
+              vavg = vsum / c;
+              vavg2 = vavg * vavg;
+              if (vavg < amin) amin = vavg;
+              else if (vavg > amax) amax = vavg;
+              asum += vavg;
+              asum2 += vavg2;
+
+              vsum2 /= c2;
+              if (vsum2 < vavg2) {
+                vsdv = 0;
+              }
+              else {
+                vsdv = sqrt(vsum2 - vavg2);
+              }
+              if (vsdv < smin) smin = vsdv;
+              else if (vsdv > smax) smax = vsdv;
+              ssum += vsdv;
+              ssum2 += (vsdv*vsdv);
+              count += 1;
+            }
+          }
+        }
+        /*
+        for (pos = 0; pos < size; pos++) {
+          value = data[pos];
+          v1 = data[pos-
+          if (value < 0.0001) {
+            count_0 += 1;
+          }
+          else {
+            if (value < min) min = value;
+            if (value > max) max = value;
+            sum_b += value;
+            sum2_b += (value * value);
+            count_b += 1;
+          }
+        }
+        */
+        aavg = asum / count;
+        asdv = sqrt((asum2 / count) - (aavg * aavg));
+        savg = ssum / count;
+        ssdv = sqrt((ssum2 / count) - (savg * savg));
+        /*
+        if (min < (avg - 2 * sdv)) min = (avg - 2 * sdv);
+        if (max > (avg + 2 * sdv)) max = (avg + 2 * sdv);
+        */
+        /*
+        stretch_1 = 0.1 / (avg - min);
+        stretch_2 = 0.9 / (max - avg);
+        */
+        printf("amin=%f amax=%f aavg=%f asdv=%f\n", amin, amax, aavg, asdv);
+        printf("smin=%f smax=%f savg=%f ssdv=%f\n", smin, smax, savg, ssdv);
+
+        smin = 2000000000;
+        smax = 0;
+        ssum = 0;
+        ssum2 = 0;
+        count = 0;
+        vsum = 0;
+        vsum2 = 0;
+        for (y = 2; y < height-2; y++) {
+          for (x = 2; x < width-2; x++) {
+            data_pos = &data[y * width + x];
+            data2_pos = &data2[y * width + x];
+            if (*data_pos > 0.0001) {
+              /*
+              v1 = *(data_pos - width - 1);
+              v2 = *(data_pos - width);
+              v3 = *(data_pos - width + 1);
+              v4 = *(data_pos - 1);
+              v5 = *data_pos;
+              v6 = *(data_pos + 1);
+              v7 = *(data_pos + width - 1);
+              v8 = *(data_pos + width);
+              v9 = *(data_pos + width + 1);
+              */
+              vsum = calc_sum(data_pos, &c);
+              /*v1+v2+v3+v4+v5+v6+v7+v8+v9;*/
+              vsum2 = calc_sum2(data_pos, &c2);
+              /*(v1*v1)+(v2*v2)+(v3*v3)+(v4*v4)+(v5*v5)+(v6*v6)+(v7*v7)+(v8*v8)+(v9*v9);*/
+              vavg = vsum/c;
+              vavg2 = vavg*vavg;
+              vsum2 /= c2;
+              if (vsum2 < vavg2) {
+                vsdv = 0;
+              }
+              else {
+                vsdv = sqrt(vsum2 - vavg2);
+              }
+              /*if (vsdv > savg) {*/
+                value = vsdv / ssdv; /* - savg */
+                *data2_pos = value;
+                if (value > 1) value = 1;
+                if (value < smin) smin = value;
+                if (value > smax) smax = value;
+                ssum += value;
+                ssum2 += value*value;
+                count += 1;
+              /*}*/
+            }
+          }
+        }
+        savg = ssum / count;
+        ssdv = sqrt((ssum2 / count) - (savg*savg));
+        printf("vmin=%f vmax=%f vavg=%f vsdv=%f\n", smin, smax, savg, ssdv);
+
+        size = width * height;
+        stretch_1 = 1.0 / (smax - smin);
+        for (pos = 0; pos < size; pos++) {
+          value = data2[pos];
+          if (value < 0.0001) {
+            value = 0.5;
+          }
+          else {
+            value = stretch_1 * (value - smin);
+          }
+          data2[pos] = value;
+        }
+      }
+    }
+    else {
+      printf("Failed to allocate memory for image");
+    }
+    free(data);
+  }
+  else {
+    printf("Failed to open file %s\n", path);
+  }
+  return dst;
+}
diff --git a/cbits/cvFiles.h b/cbits/cvFiles.h
new file mode 100644
--- /dev/null
+++ b/cbits/cvFiles.h
@@ -0,0 +1,14 @@
+#ifndef __CVFILES__
+#define __CVFILES__
+
+#include <opencv2/core/core_c.h>
+
+#ifndef M_PI
+#define M_PI           3.14159265358979323846
+#endif
+
+IplImage *read_from_tcr(const char *path, unsigned int **timestamp_array);
+IplImage *read_from_tcr_rectified(const char *path);
+IplImage *read_from_tcr_mag(const char *path);
+
+#endif
diff --git a/cbits/cvProjection.c b/cbits/cvProjection.c
new file mode 100644
--- /dev/null
+++ b/cbits/cvProjection.c
@@ -0,0 +1,150 @@
+#include "cvProjection.h"
+
+#include <stdio.h>
+#include <math.h>
+
+#ifndef M_PI
+#define M_PI           3.14159265358979323846
+#endif
+
+float calc_sum_9(float *pos, size_t width)
+{
+  float *pos1, *pos2, *pos3, *pos4, *pos5;
+  pos1 = pos - 1 * width - 1;
+  pos2 = pos - 1;
+  pos3 = pos + width - 1;
+  return 
+    *pos1 + *(pos1+1) + *(pos1+2) +
+    *pos2 + *(pos2+1) + *(pos2+2) +
+    *pos3 + *(pos3+1) + *(pos3+2); 
+}
+
+IplImage *project_polar(IplImage *src)
+{
+  IplImage *dst;
+  unsigned int w, h, x, y, px, py, minx, maxx, miny, maxy, sumx, sumy, count;
+  unsigned int new_width, new_height, src_stride, dst_stride;
+  float value, a, r, rmin, rmax, cx, cy, ratio;
+  float *src_data, *dst_data, *src_pos, *dst_pos;
+  
+  CvSize size = cvGetSize(src);
+  src_data = (float*)src->imageData;
+  src_stride = (int)(src->widthStep / sizeof(float));
+  w = size.width;
+  h = size.height;
+  
+  minx = 2000000000;
+  maxx = 0;
+  miny = 2000000000;
+  maxy = 0;
+  sumx = 0;
+  sumy = 0;
+  count = 0;
+  
+  /* find the extents of object */
+  for (y = 0; y < h; y++) {
+    src_pos = src_data + y * src_stride;
+    for (x = 0; x < w; x++, src_pos++) {
+      value = *src_pos;
+      if (value > 0.0001) {
+        if (x < minx) minx = x;
+        if (x > maxx) maxx = x;
+        if (y < miny) miny = y;
+        if (y > maxy) maxy = y;
+        sumx += x;
+        sumy += y;
+        count += 1;
+      }
+    }
+  }
+  
+  cx = (float)sumx / count;
+  cy = (float)sumy / count;
+  ratio = (float)(maxy - miny) / (float)(maxx - minx);
+  printf("w=%d h=%d cx=%f cy=%f ratio=%f\n", w, h, cx, cy, ratio);
+  
+  rmin = 2000000000;
+  rmax = 0;
+  for (y = 0; y < h; y++) {
+    src_pos = src_data + y * src_stride;
+    for (x = 0; x < w; x++, src_pos++) {
+      value = *src_pos;
+      if (value > 0.0001) {
+        px = floor(x - cx);
+        py = floor((y / ratio) - (cy / ratio));
+        r = sqrt(py*py + px*px);
+        if (r < rmin) rmin = r;
+        if (r > rmax) rmax = r;
+      }
+    }
+  }
+
+  new_height = ceil(rmax - rmin);
+  new_width = ceil(2 * M_PI * rmax);
+  printf("w=%d h=%d rmin=%f rmax=%f\n", new_width, new_height, rmin, rmax);
+
+  size.width = new_width;
+  size.height = new_height;
+  dst = cvCreateImage(size, IPL_DEPTH_32F, 1);
+  dst_data = (float*)dst->imageData;
+  dst_stride = (int)(dst->widthStep / sizeof(float));
+  
+  for (y = 0; y < new_height; y++) {
+    dst_pos = dst_data + y * dst_stride;
+    r = (float)(rmax - y);
+    for (x = 0; x < new_width; x++, dst_pos++) {
+      a = ((float)x / (float)new_width) * 2 * M_PI;
+      px = floor(cx + (r * cos(a)));
+      py = floor(((cy / ratio) + (r * sin(a))) * ratio);
+      if (px < 0 || py < 0 || px >= w || py >= h) {
+        /*printf(".");*/
+      }
+      else {
+        src_pos = src_data + py * src_stride + px;
+        value = *src_pos;
+        if (value < 0.0001) {
+          if (px > 0 && py > 0 && px < w-1 && py < h-1) {
+            value = calc_sum_9(src_pos, src_stride) / 9;
+            if (value < 0.0001) {
+              value = 0.5;
+            }
+          }
+          else {
+            value = 0.5;
+          }
+        }
+        *dst_pos = value;
+      }
+    }
+  }
+  return dst;
+}
+
+/*
+  
+  int w = (int)(2 * M_PI * r);
+  tmp = ensure32F(src);
+  dst = wrapCreateImage32F(w,h,1);
+  if (r < cx && r < cy && h <= r) {
+    float *src_data = (float *)tmp->imageData;
+    float *dst_data = (float *)dst->imageData;
+    int src_stride = (int)(tmp->widthStep / sizeof(float));
+    int dst_stride = (int)(dst->widthStep / sizeof(float));
+    int px, py;
+    double angle;
+    double radius;
+    for (int i = 0; i < w; i++) {
+      angle = ((double)i / (double)w) * (2 * M_PI);
+      for (int j = 0; j < h; j++) {
+        radius = (double)(r-j);
+        px = cx + floor(radius * cos(angle));
+        py = cy + floor(radius * sin(angle));
+        *(dst_data + j * dst_stride + i) = *(src_data + py * src_stride + px);
+      }
+    }
+  }
+  else {
+    printf("image too small\n");
+  }
+*/
+
diff --git a/cbits/cvProjection.h b/cbits/cvProjection.h
new file mode 100644
--- /dev/null
+++ b/cbits/cvProjection.h
@@ -0,0 +1,8 @@
+#ifndef __CVITERATORS__
+#define __CVITERATORS__
+
+#include <opencv2/core/core_c.h>
+
+IplImage *project_polar(IplImage *src);
+
+#endif
diff --git a/cbits/cvWrapLEO.c b/cbits/cvWrapLEO.c
--- a/cbits/cvWrapLEO.c
+++ b/cbits/cvWrapLEO.c
@@ -280,6 +280,36 @@
 
 // Return image that is IPL_DEPTH_32F version of
 // given src
+IplImage* ensure64F(const IplImage *src)
+{
+ CvSize size;
+ IplImage *result;
+ int channels = src->nChannels;
+ int dstDepth = IPL_DEPTH_64F;
+ size = cvGetSize(src);
+ result = cvCreateImage(size,dstDepth,channels);
+
+ switch(src->depth) {
+  case IPL_DEPTH_32F:
+  case IPL_DEPTH_64F:
+   cvConvertScale(src,result,1,0); // Scale the values to [0,255]
+   return result;
+  case IPL_DEPTH_8U:
+  case IPL_DEPTH_8S:
+   cvConvertScale(src,result,1.0/255.0,0);
+   return result;
+  case IPL_DEPTH_16S:
+   cvConvertScale(src,result,1.0/65535.0,0);
+   return result;
+  case IPL_DEPTH_32S:
+   cvConvertScale(src,result,1.0/4294967295.0,0);
+   return result;
+  default:
+   printf("Cannot convert to floating image");
+   abort();
+ }
+}
+
 IplImage* ensure32F(const IplImage *src)
 {
  CvSize size;
@@ -307,7 +337,6 @@
   default:
    printf("Cannot convert to floating image");
    abort();
-
  }
 }
 
@@ -375,7 +404,10 @@
  free(pts);
 }
 
-
+void wrapDrawEllipse(IplImage *img, int x, int y, int r1, int r2, float a, float a1, float a2, float r, float g, float b, int thickness)
+{
+  cvEllipse(img, cvPoint(x,y),cvSize(r1,r2),a,a1,a2,CV_RGB(r,g,b),thickness,8,0);
+}
 
 int getImageWidth(IplImage *img)
 {
@@ -454,6 +486,30 @@
 // printf("Done!\n"); fflush(stdout);
 }
 
+// Assuming a is the bigger image
+void blitShadow(IplImage *a, IplImage *b)
+{
+ CvSize sa = cvGetSize(a);
+ CvSize sb = cvGetSize(b);
+
+ for ( int i=0; i<sb.width; i++ )
+  for ( int j=0; j<sb.height; j++ )
+   FGET(a,i,j) = FGET(b,i,j);
+
+ for ( int i=sb.width; i<sa.width; i++ )
+  for ( int j=0; j<sb.height; j++ )
+   FGET(a,i,j) = FGET(b,sb.width-1,j);
+
+ for ( int i=0; i<sb.width; i++ )
+  for ( int j=sb.height; j<sa.height; j++ )
+   FGET(a,i,j) = FGET(b,i,sb.height-1);
+
+ for ( int i=sb.width; i<sa.width; i++ )
+  for ( int j=sb.height; j<sa.height; j++ )
+   FGET(a,i,j) = FGET(b,sb.width-1,sb.height-1);
+
+}
+
 IplImage* makeEvenDown(IplImage *src)
 {
  CvSize size = cvGetSize(src);
@@ -1780,7 +1836,7 @@
 //@-node:aleator.20051207074905:LBP
 //@+node:aleator.20051109102750:Selective Average
 // Assuming grayscale image calculate local selective average of point x y
-double calcSelectiveAvg(IplImage *img,double t
+inline double calcSelectiveAvg(IplImage *img,double t
                                    ,int x, int y
                                    ,int wwidth, int wheight)
 {
@@ -2128,7 +2184,7 @@
 //@+node:aleator.20060413093124.1:sizeFilter
 IplImage* sizeFilter(IplImage *src, double minSize, double maxSize)
 {
-    IplImage* dst = cvCreateImage( cvGetSize(src), IPL_DEPTH_32F, 1 );
+    IplImage* dst = cvCreateImage( cvGetSize(src), IPL_DEPTH_8U, 1 );
     CvMemStorage* storage = cvCreateMemStorage(0);
     CvSeq* contour = 0;
 
@@ -2139,7 +2195,7 @@
     {
         double area=fabs(cvContourArea(contour,CV_WHOLE_SEQ,0));
         if (area <=minSize || area >= maxSize) continue;
-        CvScalar color = cvScalar(1,1,1,1);
+        CvScalar color = cvScalar(255,255,255,255);
         cvDrawContours( dst, contour, color, color, -1, CV_FILLED, 8,
             cvPoint(0,0));
     }
@@ -2384,6 +2440,7 @@
     index += seq->elem_size;
     CV_NEXT_SEQ_ELEM( seq->elem_size, reader );
 }}
+
 
 #ifndef OpenCV24
 void wrapExtractMSER( CvArr* _img, CvArr* _mask, CvSeq** contours, CvMemStorage* storage, CvMSERParams *params ){
diff --git a/cbits/cvWrapLEO.h b/cbits/cvWrapLEO.h
--- a/cbits/cvWrapLEO.h
+++ b/cbits/cvWrapLEO.h
@@ -57,6 +57,8 @@
 
 void wrapFillPolygon(IplImage *img, int pc, int *xs, int *ys, float r, float g, float b);
 
+void wrapDrawEllipse(IplImage *img, int x, int y, int r1, int r2, float a, float a1, float a2, float r, float g, float b, int thickness);
+
 void wrapMatMul(int w, int h, double *mat
                , double *vec, double *t);
 
@@ -320,6 +322,13 @@
 
 void wrapBoundingRect(CvArr* points, int update, CvRect *r)
 { *r = cvBoundingRect(points, update); }
+
+void wrapWarpPerspective(const CvArr* src, CvArr* dst, const CvMat* mapMatrix )
+{ cvWarpPerspective( src,  dst,  mapMatrix, CV_INTER_LINEAR+CV_WARP_FILL_OUTLIERS, cvScalarAll(0)); }
+
+void blitShadow(IplImage *a, IplImage *b);
+
+IplImage* ensure64F(const IplImage *src);
 
 #endif
 //@-node:aleator.20050908101148.2:@thin cvWrapLEO.h
diff --git a/examples/SizeFilt.hs b/examples/SizeFilt.hs
new file mode 100644
--- /dev/null
+++ b/examples/SizeFilt.hs
@@ -0,0 +1,10 @@
+import CV.Image
+import CV.ConnectedComponents
+import qualified CV.ImageMath as IM
+
+main = do
+    Just x <- loadImage "meanshiftsample.png"
+    let y  = selectSizedComponents 10 100000000 $ IM.moreThan 0.5  x
+
+    saveImage "size0.png" $ IM.moreThan 0.5 x
+    saveImage "size.png" y
diff --git a/examples/SnakeTest.hs b/examples/SnakeTest.hs
new file mode 100644
--- /dev/null
+++ b/examples/SnakeTest.hs
@@ -0,0 +1,12 @@
+module Main where
+import CV.Image
+import CV.Tracking
+import CV.Bindings.Types
+
+main = do
+    Just x <- loadImage "meanshiftsample2.png"
+    let initial = [(round $ 100+80*cos x, round $ 100+80*sin x) | x <- [0,0.1..2*pi]]
+    mapM_ print initial
+    new <- snake (unsafeImageTo8Bit x) initial 0.1 0.1 0.24 (3,3) (EPS 0.1) False
+    mapM_ print new
+
diff --git a/examples/homography.hs b/examples/homography.hs
--- a/examples/homography.hs
+++ b/examples/homography.hs
@@ -4,17 +4,15 @@
 import CV.Transforms
 
 main = do
-  let source = M.fromList (5,2) [1,1
-                                ,1,2
-                                ,2,1
-                                ,5,5
-                                ,2,4]
-      target = M.fromList (5,2) [1,1
-                                ,1,2
-                                ,2,1
-                                ,5,5
-                                ,1900,4]
+  let mat = M.fromList (5,2) . concatMap (\(x,y) -> [x,y])
+      source = [(1,1)
+               ,(1,2)
+               ,(2,1)
+               ,(5,5)
+               ,(2,4)]
+      tr = map (\(x,y) -> (x+100,y+100)) source
+      sc = map (\(x,y) -> (x*100,y*100)) source
 
-  print (getHomography' source target Default 0.1)
-  print (getHomography' source target Ransac 0.001)
-  print (getHomography' source target LMeds 0.001)
+  print (getHomography' (mat source) (mat $ tr ) Ransac 0.1)
+  print (getHomography' (mat source) (mat $ sc ) Ransac 0.1)
+  print (getHomography' (mat source) (mat source) LMeds 0.1)
diff --git a/examples/newContours.hs b/examples/newContours.hs
deleted file mode 100644
--- a/examples/newContours.hs
+++ /dev/null
@@ -1,11 +0,0 @@
-module Main where
-import CV.Image
-import CV.ConnectedComponents
-import CV.Bindings.Types
-
-
-main = do
-   Just x <- loadImage "Circle.png" >>= return . fmap unsafeImageTo8Bit
-   cs <- contours x (C'CvPoint 0 0) 0 0 0
-   print cs
-   print "done"