## SPOJ 6044. Minimum Diameter Circle [ Enclosing Circle Problem ]

SPOJ 6044. Minimum Diameter Circle problem is also know as enclosing circle problem. Accoriding to Wikipedia ” The smallest circle problem or minimum covering circle problem is a mathematical problem of computing the smallest circle that contains all of a given set of points in the Euclidean plane. The corresponding problem in n-dimensional space, the smallest bounding sphere problem, is to compute the smallest n-sphere that contains all of a given set of points. The smallest circle problem was initially proposed by the English mathematician James Joseph Sylvester in 1857″. Couple of algorithm explained here and implementation Pr. Chrystal algorithm. Accepted Haskell code. You can run this code on ideone.

This algorithm is due to Pr. Chrystal (1885). The Algorithm Foremost, observe that the minimal enclosing circle, MEC, is entirely determined by the Convex Hull of a given set of point. The reason is the points of the set touched by the minimal enclosing circle are always on the convex hull of the set. First step of the algorithm is to compute the convex hull, H, of the points in linear time. Clearly, we can do this since points are kept ordered by x-coordinate. Call the convex hull H and the number of convex hull vertices h. Next step is to pick any side, say S, of the convex hull. The main loop of the algorithm is as follows: Main Loop For each vertex of H other than those of S, we compute the angle subtended by S. The minimum such angle occurs at vertex v, and the value of the minimum angle is α (alpha). 1 . IF α ≥ 90 degrees THEN We are done! /* The minimal enclosing circle is determined by the diametric circle of S */ IF α < 90 degrees THEN Goto step 2 below. 2. Since α < 90 degrees, check the remaining vertices of the triangle formed by S and v. IF no vertices are obtuse THEN We are done! /* The MEC is determined by the vertices of S and the vertex v */ 3. If one of the other angles of the triangle formed by S and v is obtuse, then set S to be the side opposite the obtuse angle and repeat the main loop of the algorithm (i.e., Goto step 1 above).

import Data.List import qualified Data.Sequence as DS import Text.Printf import Data.Function ( on ) data Point a = P a a deriving ( Show , Ord , Eq ) data Turn = S | L | R deriving ( Show , Eq , Ord , Enum ) -- straight left right --start of convex hull http://en.wikipedia.org/wiki/Graham_scan {-- 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 convexHull ::( Num a , Ord a ) => [ Point a ] -> [ Point a ] convexHull xs = reverse . findHull [y,x] $ ys where (x:y:ys) = sortByangle.findMinx $ 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 --} --start of monotone hull give .04 sec lead over graham scan 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 sortPoint :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ] sortPoint xs = sortBy ( \ x y -> compPoint x y ) 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 hullComputation :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ] -> [ Point a ] hullComputation [x] ( z:ys ) = hullComputation [z,x] ys hullComputation xs [] = xs hullComputation ( y : x : xs ) ( z : ys ) | findTurn x y z == R = hullComputation ( x:xs ) ( z : ys ) | findTurn x y z == S = hullComputation ( x:xs ) ( z : ys ) | otherwise = hullComputation ( z : y : x : xs ) ys convexHull :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ] convexHull xs = final where txs = sortPoint xs ( x : y : ys ) = txs lhull = hullComputation [y,x] ys ( x': y' : xs' ) = reverse txs uhull = hullComputation [ y' , x' ] xs' final = nub $ ( reverse lhull ) ++ uhull --end of convex hull --start of finding point algorithm http://www.personal.kent.edu/~rmuhamma/Compgeometry/MyCG/CG-Applets/Center/centercli.htm Applet’s Algorithm findAngle :: ( Num a , Ord a , Floating a ) => Point a -> Point a -> Point a -> ( Point a , Point a , Point a , a ) findAngle u@(P x0 y0 ) v@(P x1 y1 ) t@(P x2 y2) | u == t || v == t = ( u , v , t , 10 * pi ) -- two points are same so set the angle more than pi | otherwise = ( u , v, t , ang ) where ang = acos $ ( b + c - a ) / ( 2.0 * ( sqrt $ b * c ) ) where sqrDist ( P x y ) ( P x' y' ) = ( x - x' ) ^ 2 + ( y - y' ) ^ 2 [ b , c , a ] = map ( uncurry sqrDist ) [ ( u , t ) , ( v , t ) , ( u , v ) ] findPoints :: ( Num a , Ord a , Floating a ) => Point a -> Point a -> [ Point a ] -> ( Point a , a ) findPoints u v xs | 2 * theta >= pi = circle2Points a b --( a , b , t , theta ) | and [ 2 * alpha <= pi , 2 * beta <= pi ] = circle3Points a b t --( a , b , t , theta ) | otherwise = if 2 * alpha > pi then findPoints v t xs else findPoints u t xs where ( a , b , t , theta ) = minimumBy ( on compare fn ) . map ( findAngle u v ) $ xs fn ( _ , _ , _ , t ) = t ( _ , _ , _ , alpha ) = findAngle v t u --angle between v u t angle subtended at u by v t ( _ , _ , _ , beta ) = findAngle u t v -- angle between u v t angle subtended at v by u t --end of finding three points algorithm --find the circle through three points http://paulbourke.net/geometry/circlefrom3/ circle2Points :: ( Num a , Ord a , Floating a ) => Point a -> Point a -> ( Point a , a ) circle2Points ( P x0 y0 ) ( P x1 y1 ) = ( P x y , r ) where x = ( x0 + x1 ) / 2.0 y = ( y0 + y1 ) / 2.0 r = sqrt $ ( x0 - x1 ) ^ 2 + ( y0 - y1 ) ^ 2 circle3Points :: ( Num a , Ord a , Floating a ) => Point a -> Point a -> Point a -> ( Point a , a ) --( center , radius ) circle3Points u@(P x1 y1 ) v@(P x2 y2 ) t@(P x3 y3 ) | x2 == x1 = circle3Points u t v | x3 == x2 = circle3Points v u t | otherwise = ( P x y , 2 * r ) where m1 = ( y2 - y1 ) / ( x2 - x1 ) m2 = ( y3 - y2 ) / ( x3 - x2 ) x = ( m1 * m2 * ( y1 - y3 ) + m2 * ( x1 + x2 ) - m1 * ( x2 + x3 ) ) / ( 2 * ( m2 - m1 ) ) y = if y2 /= y1 then ( ( x1 + x2 - 2 * x ) / ( 2 * m1 ) ) + ( ( y1 + y2 ) / 2.0 ) else ( ( x2 + x3 - 2 * x ) / ( 2 * m2 ) ) + ( ( y2 + y3 ) / 2.0 ) r = sqrt $ ( x - x1 ) ^2 + ( y - y1 ) ^ 2 --start of SPOJ code format::(Num a , Ord a ) => [[a]] -> [Point a] format xs = map (\[x0 , y0] -> P x0 y0 ) xs readInt ::( Num a , Read a ) => String -> a readInt = read solve :: ( Num a , Ord a , Floating a ) => [ Point a ] -> ( Point a , a ) solve [ P x0 y0 ] = ( P x0 y0 , 0 ) --in case of one point solve [ P x0 y0 , P x1 y1 ] = circle2Points ( P x0 y0 ) ( P x1 y1 ) -- in case of two points solve xs = findPoints x y t where t@( x : y : ys ) = convexHull xs final :: ( Num a , Ord a , Floating a ) => ( Point a , a ) -> a final ( t , w ) | w == 0 = 0 | otherwise = w main = interact $ ( printf "%.2f\n" :: Double -> String ) . final . solve . convexHull . format . map ( map readInt . words ) . tail . lines

Advertisements

No comments yet.

## Leave a Reply