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

## 1 Comment »

### Leave a Reply

Advertisements

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 |