My Weblog

Blog about programming and math

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

July 24, 2011 - Posted by | Programming | , , , ,

No comments yet.

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: