# My Weblog

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

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

```