My Weblog

Blog about programming and math

Rotating calipers [ The diameter of a convex polygon ]

Rotating calipers is great tool in computational geometry similar to sweep line . According to Wikipedia ” In computational geometry, rotating calipers is a method used to construct efficient algorithms for a number of problems”. Given code is Haskell transformation of pseudo code described at Wikipedia and computes the square of the diameter of a convex polygon . I tried to test this code on TFOSS but time limit is very tight so its hard for Haskell to get accepted.I generated some random test cases and it worked fine. A excellent paper ” Solving geometric problems with the rotating calipers ” investigates some problems of computational geometry which can be solved by this tool.

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 )  

--square of dist between two points 

distPoints :: ( Num a , Ord a  ) => Point a -> Point a -> a
distPoints ( P x1 y1 ) ( P x2 y2 ) =  ( x1 - x2 ) ^ 2 + ( y1 - y2 ) ^ 2 

--rotating caliipers 

rotCal :: ( Num a , Ord a , Floating a ) => Array Int ( Point a )  -> a -> Int -> Int -> Vector a -> Vector a -> a -> Int -> a
rotCal arr ang  pa pb ca@( V ax ay ) cb@( V bx by ) dia n
   | ang > pi = dia
   | otherwise = rotCal arr ang' pa' pb' ca' cb' dia' 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 )
	t1 = angVectors ca ( V ( x2 - x1 ) ( y2 - y1 ) )
	t2 = angVectors cb ( V ( x4 - x3 ) ( y4 - y3 ) )
	ca' = rotVector ca  $ min t1 t2
	cb' = rotVector cb  $ min t1 t2
	ang' = ang + min t1 t2
	( pa' , pb' )  = if t1 < t2 then ( mod ( pa + 1 ) n  , pb ) else ( pa , mod ( pb + 1 ) n  )
	dia' = max dia $ distPoints ( arr ! pa' ) ( arr ! pb' )

solve :: ( Num a , Ord a , Floating a ) => [ Point a ] -> a
solve [] = 0
solve [ p ] = 0
solve [ p1 , p2 ] =  distPoints p1 $ p2
solve arr =  rotCal arr' 0 pa pb ( V 1 0 ) ( V (-1) 0 ) dia   n where
	   y1 = minimumBy ( on  compare fN  ) arr
	   y2 = maximumBy ( on  compare fN  ) arr
	   pa = fromJust . findIndex (  == y1 ) $ arr
	   pb = fromJust . findIndex (  == y2 ) $ arr
	   dia = distPoints ( arr !! pa ) ( arr !! pb )
	   n = length arr
	   arr' = listArray ( 0 , n ) arr
	   fN ( P x y ) = y 

 --end of rotating caliper 

--spoj code for testing but time limit is very tight so hard to get accepted in Haskell 

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

format :: ( Num a , Ord a , Floating a ) => [ Int ] -> [ [ Point a ]]
format [] = []
format (x:xs ) =  t : format b  where
	( a , b ) = splitAt ( 2 * x ) xs
	t = helpFormat a where
	    helpFormat [] = []
	    helpFormat ( x' : y' : xs' ) = ( P ( fromIntegral x' ) ( fromIntegral  y' ) ) : helpFormat xs'

readD :: BS.ByteString -> Int
readD = fst . fromJust . BS.readInt

main = BS.interact $ BS.unlines . map  ( BS.pack . ( printf "%.0f" :: Double -> String ) .final ) . format . concat . map  ( map  readD . BS.words ) . tail  . BS.lines  

--end of spoj code
Advertisements

August 4, 2011 - Posted by | Programming | , , ,

1 Comment »

  1. i implemented this algo but getting wa in http://www.spoj.com/problems/TFOSS/
    can u tell error in this code : http://ideone.com/SRI0qi
    seems to work fine to me

    Comment by joker | December 10, 2012 | Reply


Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: