content: {-#LANGUAGE ViewPatterns, ParallelListComp#-}
  module Main where
  
  import CV.Image
  
  import CV.ColourUtils
  import CV.Morphology
  import CV.Filters
  import CV.Thresholding
  import CV.Edges
  import CV.ImageMathOp
  import qualified CV.ImageMath as IM
  import CV.FunnyStatistics
  import CV.ConnectedComponents
  import CV.Transforms
  import CV.Matrix
  import CV.Fitting
  import System.Environment
  import CV.Pixelwise
  import Graphics.Gnuplot.Simple
  import Data.List
  import Data.Ord
  
  hdef d (x:xs) = x
  hdef d _      = d
  
  horMax :: Image GrayScale D32 -> [Double]
  horMax (pixels) = [fromIntegral $ snd $ maximumBy (comparing value) [(i,j) | j<-[0..height-1]] | i <- [0..width-1] ]
      where
          value (x,y) = getPixel (x,y) pixels
          (width,height) = getSize pixels
  
  getMarks :: Image GrayScale D32 -> [(Float,Float)]
  getMarks pixels = [(realToFrac i,realToFrac j) | j <- [0..height-1]
                           , i <- [0..width-1]
                           , value (i,j) > 0.5]
      where
          value (x,y) = getPixel (x,y) pixels
          value (x,y) = getPixel (x,y) pixels
  
  ntimes n op = (!! n) . iterate op
  
  smooth [x] = [x]
  smooth [x,y] = [(x+y) / 2]
  smooth (x:y:xs) = (x+y) / 2:smooth (y:xs)
  
  conv mask = map (prod mask) . tails
      where
       prod xs ys = sum (zipWith (*) xs ys)
  
  to32F = unsafeImageTo32F
  
  main = do
      Just x <- getArgs >>= loadImage . hdef (error "No file given!") >>= return . fmap (enlarge n )
      let
          oper msk = to32F . nibbly 1 0.01 . IM.abs . sobel msk s3
          up op = foldl1 (\a b -> (pyrUp a) `IM.min` b)
                   . reverse
                   . take n
                   . map op
                   . iterate pyrDown
  
          se = structuringElement (9,9) (4,4) EllipseShape
          pyrd x =  (up (oper (1,0)) x) #+ (up (oper (0,1)) x)
  
          thd =  dilate se 2 . unsafeImageTo32F . nibbly 0.9 0.001 . IM.invert 
          clean = selectSizedComponents (1e4) (1e10) . IM.moreThan 0.2 $ thd x #* (close se (pyrd x))
  
          final = to32F clean #* (IM.abs (sobel (1,0) s1 x) #+ IM.abs (sobel (0,1) s1 x))
          hm = getMarks $ to32F clean
          ell = fitEllipse $ fromList (1,length hm) $ hm
          rotated = rotate (realToFrac $ pi/2+2*pi*(angle ell/360)) final
      saveImage "Combine.png" rotated
   where n = 5
  
  nonMaxSupress w = reverse . nonMaxSupress' w . reverse . nonMaxSupress' w
  nonMaxSupress' w = map supress . map (take w) . tails
      where
          supress [] = 0
          supress [x] = x
          supress (x:xs) | all (<x) xs = x
                         | otherwise   = 0
  
  threshold t x | x <t = 0
                | otherwise = 10
  