patternMinor
Enclosing Circle Problem implementation in Haskell
Viewed 0 times
circleproblemhaskellimplementationenclosing
Problem
I implemented the algorithm by Pr. Chrystal described here in Haskell so could someone please tell me: Do i have implemented this algorithm correctly?
Initial calling of findPoints takes the first and second point of convex hull and the list of points on convex hull.
My second question is: I am trying to solve the problem on sphere online judge [link of problem in code below, because I can not post more than 2 links] using this algorithm but I'm getting the wrong answer.
The code for the problems is on ideone. Why am I getting the wrong answer.
Initial calling of findPoints takes the first and second point of convex hull and the list of points on convex hull.
My second question is: I am trying to solve the problem on sphere online judge [link of problem in code below, because I can not post more than 2 links] using this algorithm but I'm getting the wrong answer.
The code for the problems is on ideone. Why am I getting the wrong answer.
--http://www.spoj.pl/problems/QCJ4/
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 * sb * sc ) ) where
b = ( x0 - x2 ) ^ 2 + ( y0 - y2 ) ^ 2
c = ( x1 - x2 ) ^ 2 + ( y1 - y2 ) ^ 2
a = ( x0 - x1 ) ^ 2 + ( y0 - y1 ) ^ 2
sb = sqrt b
sc = sqrt c
findPoints :: ( Num a , Ord a , Floating a ) => Point a -> Point a -> [ Point a ] -> ( Point a , Point a , Point a , a )
findPoints u v xs
| 2 * theta >= pi = ( a , b , t , theta )
| and [ 2 * alpha pi then findPoints v t xs else findPoints u t xs
where
( a , b , t , theta ) = minimumBy ( \(_,_,_, t1 ) ( _ , _ , _ ,t2 ) -> compare t1 t2 ) . map ( findAngle u v ) $ xs
( _ , _ , _ , 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 tSolution
As per request by Mihai , here is the accepted code of spoj problem . The only error in previous code was " The minimal enclosing circle is determined by the diametric circle of S ". Instead of finding a circle passing through two points , i was calculating circle from three points.
```
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 ) [ 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 ) [ 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 = ( init lhull ) ++ ( init 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 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
| ot
```
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 ) [ 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 ) [ 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 = ( init lhull ) ++ ( init 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 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
| ot
Code Snippets
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'Context
StackExchange Code Review Q#3603, answer score: 5
Revisions (0)
No revisions yet.