CV 0.3.6.0 → 0.3.7

raw patch · 33 files changed

+2707/−379 lines, 33 filesdep +mwc-randomdep +primitivedep −randomPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: mwc-random, primitive

Dependencies removed: random

API changes (from Hackage documentation)

Files

• CV.cabal +18/−6
• CV/Bindings/Core.hsc +11/−91
• CV/Bindings/Files.hsc +20/−0
• CV/Bindings/ImgProc.hsc +0/−1
• CV/Bindings/Matrix.hsc +6/−0
• CV/Bindings/Projection.hsc +13/−0
• CV/Bindings/Types.hsc +13/−0
• CV/Conversions.hs +5/−5
• CV/DFT.hs +40/−122
• CV/Drawing.chs +96/−56
• CV/Files.hs +28/−0
• CV/Histogram.chs +0/−1
• CV/Image.chs +55/−36
• CV/ImageOp.hs +4/−0
• CV/Matrix.hs +11/−1
• CV/Operations.hs +26/−3
• CV/Projection.hs +12/−0
• CV/Sampling.hs +25/−24
• CV/TemplateMatching.chs +6/−1
• CV/Thresholding.hs +10/−0
• CV/Transforms.chs +56/−2
• Utils/GeometryClass.hs +11/−0
• Utils/Pointer.hs +2/−1
• cbits/cvFiles.c +1963/−0
• cbits/cvFiles.h +14/−0
• cbits/cvProjection.c +150/−0
• cbits/cvProjection.h +8/−0
• cbits/cvWrapLEO.c +62/−5
• cbits/cvWrapLEO.h +9/−0
• examples/SizeFilt.hs +10/−0
• examples/SnakeTest.hs +12/−0
• examples/homography.hs +11/−13
• examples/newContours.hs +0/−11

@@ -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

@@ -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

@@ -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)

@@ -52,7 +52,6 @@             (fromIntegral r)             (cBorderType border)             (realToFrac value)- -- moments  -- | Calculates all spatial and central moments up to the 3rd order

@@ -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

@@ -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)

@@ -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

@@ -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)))))) 

@@ -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

@@ -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) -            + 

@@ -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)

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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)

@@ -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) 

@@ -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

@@ -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

@@ -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

@@ -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+

@@ -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 

@@ -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;+}

@@ -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

@@ -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");+  }+*/+

@@ -0,0 +1,8 @@+#ifndef __CVITERATORS__+#define __CVITERATORS__++#include <opencv2/core/core_c.h>++IplImage *project_polar(IplImage *src);++#endif

@@ -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 ){

@@ -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

@@ -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

@@ -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+

@@ -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)

@@ -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"