Hello all , I am trying to understand rotating calipers [
http://en.wikipedia.org/wiki/Rotating_calipers ] but i am not sure if
understood this algorithm correctly . I tried to use the almost same
algorithm given on wiki but with four calipers to solve the problem
[ http://cgm.cs.mcgill.ca/~orm/rotcal.html ]. My approach is find
xminP, xmaxP, yminP ymaxP and their corresponding calipers will be ( 0
* i - j ) , ( o * i + j ) , ( i + 0 * j ) and ( -i + 0 * j ). I
implemented the algorithm in Haskell but its not working . I am not
sure if i have followed the wiki algorithm correctly and could some
one please tell me what is wrong with implementation. It would be
great if some one can explain this algorithm in   pseudo  code which
explains the rotating caliper and their implementation details . In
case of indentation , see here [ http://hpaste.org/49907 ] .
Thank you
Mukesh Tiwari

import Data.List
import Data.Array
import Data.Maybe
import Data.Function
import Text.Printf
import qualified Data.ByteString.Char8 as BS

data Point a = P a a deriving ( Show , Ord , Eq )
data Vector a = V a a deriving ( Show , Ord , Eq )
data Turn = S | L | R deriving ( Show , Eq , Ord , Enum  )

--start of convex hull

compPoint :: ( Num  a , Ord a ) => Point a -> Point a -> Ordering
compPoint ( P x1 y1 ) ( P x2 y2 )
  | compare x1 x2 == EQ = compare y1 y2
  | otherwise = compare x1 x2

findMinx :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ]
findMinx xs = sortBy ( \x  y  -> compPoint  x y  ) xs

compAngle ::(Num a , Ord a ) => Point a -> Point a -> Point a ->
Ordering
compAngle ( P x1 y1 ) ( P x2 y2 ) ( P x0 y0 ) = compare ( (  y1 - y0 )
* ( x2 - x0 )  ) ( ( y2 - y0) * ( x1 - x0 ) )

sortByangle :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ]
sortByangle (z:xs) = z : sortBy ( \x y -> compAngle x y z ) xs

findTurn :: ( Num a , Ord a , Eq a ) => Point a -> Point a -> Point a -
> Turn
findTurn ( P x0 y0 ) ( P x1 y1 ) ( P x2 y2 )
 | ( y1 - y0 ) * ( x2- x0 ) < ( y2 - y0 ) * ( x1 - x0 ) = L
 | ( y1 - y0 ) * ( x2- x0 ) == ( y2 - y0 ) * ( x1 - x0 ) = S
 | otherwise = R

findHull :: ( Num a , Ord a  )  => [ Point a ] ->   [ Point a ] ->
[ Point a ]
findHull [x]  ( z : ys )  = findHull [ z , x ]  ys  --incase of second
point  on line from x to z
findHull xs  [] = xs
findHull ( y : x : xs )  ( z : ys )
  | findTurn x y z == R = findHull (  x : xs )   ( z:ys )
  | findTurn x y z == S = findHull (  x : xs )   ( z:ys )
  | otherwise = findHull ( z : y : x : xs  )   ys


convexHull ::( Num a , Ord a )  => [ Point a ] -> [ Point a ]
convexHull xs = reverse . findHull [ y , x ]  $ ys where
        ( x : y : ys ) = sortByangle . findMinx $ xs


--end of convex hull

--start of rotating caliper part http://en.wikipedia.org/wiki/Rotating_calipers
--dot product for getting angle

angVectors :: ( Num a , Ord a , Floating a ) => Vector a -> Vector a -
> a
angVectors ( V ax ay ) ( V bx by ) = theta where
    dot = ax * bx + ay * by
    a = sqrt $ ax ^ 2 + ay ^ 2
    b = sqrt $ bx ^ 2 + by ^ 2
    theta = acos $ dot / a / b

--rotate the vector x y by angle t

rotVector :: ( Num a , Ord a , Floating a ) => Vector a -> a -> Vector
a
rotVector ( V x y ) t = V ( x * cos t - y * sin t ) ( x * sin t + y *
cos t )

--dist between two parallel vectors

distVec :: ( Num a , Ord a , Floating a ) => Vector a -> Vector a ->
a
distVec ( V x1 y1 ) ( V x2 y2 ) = sqrt $ ( x1 - x2 ) ^ 2 + ( y1 - y2 )
^ 2
--rotating caliipers

rotCal :: ( Num a , Ord a , Floating a ) => Array Int ( Point a )  ->
a -> [ Int ] -> [ Vector a ] -> a -> Int -> a
rotCal arr ang  [ pa , pb , qa , qb] [ cpa , cpb , cqa , cqb ] area  n
   | 2 * ang > pi = area
   | otherwise = rotCal arr ang' [ pa' , pb' , qa' , qb' ] [ cpa' ,
cpb' , cqa' , cqb' ] area' n where
        P x1 y1 = arr ! pa
        P x2 y2 = arr ! ( mod ( pa + 1 ) n )
        P x3 y3 = arr ! pb
        P x4 y4 = arr ! ( mod ( pb + 1 ) n )

        P x5 y5 = arr ! qa
        P x6 y6 = arr ! ( mod ( qa + 1 ) n )
        P x7 y7 = arr ! qb
        P x8 y8 = arr ! ( mod ( qb + 1 ) n )

        t1 = angVectors cpa ( V ( x2 - x1 ) ( y2 - y1 ) )
        t2 = angVectors cpb ( V ( x4 - x3 ) ( y4 - y3 ) )
        t3 = angVectors cqa ( V ( x6 - x5 ) ( y6 - y5 ) )
        t4 = angVectors cqb ( V ( x8 - x7 ) ( y8 - y7 ) )
        t = minimum [ t1 , t2 , t3 , t4 ]

        cpa' = rotVector cpa  t
        cpb' = rotVector cpb  t
        cqa' = rotVector cqa  t
        cqb' = rotVector cqb  t

        ang' = ang + t
        ( pa' , pb' , qa' , qb' ) = fN [ t1 , t2 , t3 , t4 ] t where
                fN [ t1 , t2 , t3 , t4 ] t
                   | t == t1 = ( mod ( pa + 1 ) n , pb , qa , qb )
                   | t == t2 = ( pa , mod ( pb + 1 ) n , qa , qb )
                   | t == t3 = ( pa , pb , mod ( qa + 1 ) n , qb )
                   | otherwise = ( pa , pb , qa , mod ( qb + 1 ) n )

        width = distVec cpa' cpb'
        length = distVec cqa' cqb'
        area' = min area $ length * width

solve :: ( Num a , Ord a , Floating a ) => [ Point a ] -> a
solve [] = 0
solve [ p ] = 0
solve [ p1 , p2 ] =  0
solve [ p1 , p2 , p3 ] = 0
solve arr =  rotCal arr' 0 [ pa , pb , qa , qb ] [ cpa , cpb , cqa ,
cqb ] area   n where
           y1 = minimumBy ( on  compare fN1  ) arr
           y2 = maximumBy ( on  compare fN1  ) arr
           x1 = minimumBy ( on  compare fN2  ) arr
           x2 = maximumBy ( on  compare fN2  ) arr
           pa = fromJust . findIndex (  == y1 ) $ arr
           pb = fromJust . findIndex (  == y2 ) $ arr
           qa = fromJust . findIndex (  == x1 ) $ arr
           qb = fromJust . findIndex (  == x2 ) $ arr
           cpa = V 1 0
           cpb = V ( -1 ) 0
           cqa = V 0 ( -1 )
           cqb = V 0 1
           area = 1e9
           n = length arr
           arr' = listArray ( 0 , n ) arr
           fN1 ( P x y ) = y
           fN2 ( P x y ) = x

--end of rotating caliper



final :: ( Num a , Ord a , Floating a ) => [ Point a ] -> a
final [] = 0
final [ p ] = 0
final [ p1 , p2 ] =  0
final [ p1 , p2 , p3 ] = 0
final arr = solve . convexHull $ arr

_______________________________________________
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

Reply via email to