5
\$\begingroup\$

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.

--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 , 2 * beta <= pi ]   = ( 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 ( \(_,_,_, 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 t
\$\endgroup\$
  • \$\begingroup\$ I have solved this problem. Accepted code of problem QCJ4 spoj on ideone. I hope Mihai don't need to give away his reputation points as bounty to any one :). \$\endgroup\$ – keep_learning Jul 29 '11 at 16:05
  • \$\begingroup\$ Can you post solution below? (I'll give you the bounty in this case) \$\endgroup\$ – Mihai Maruseac Aug 2 '11 at 7:30
  • \$\begingroup\$ @keep_learning: Now I got the bounty, which is not okay as I just posted some cosmetic changes, but AFAIK there is no way to give it back. If there is no other solution for that, at least I could start a 50 points bounty for a question of your choice... \$\endgroup\$ – Landei Aug 3 '11 at 10:43
  • \$\begingroup\$ @Mihai posted solution here \$\endgroup\$ – keep_learning Aug 4 '11 at 14:00
  • \$\begingroup\$ @Landel I think your code lead me to write more concise Haskell code so you can think of these bounties as taken from me :) \$\endgroup\$ – keep_learning Aug 4 '11 at 14:02
5
\$\begingroup\$

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 ) < ( 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 = ( 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 , 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   
\$\endgroup\$
3
+50
\$\begingroup\$

Not an answer, just a small improvement (not tested):

import Data.Function(on)

findAngle :: (Num a , Ord a , Floating a ) => Point a -> Point a -> Point  a  -> ( Point a , Point a , Point  a , a ) 
findAngle u v t  
  | 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 * (sqrt $ b * c))   
        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 , Point a , Point a , a ) 
findPoints u v xs 
  | 2 * theta >= pi   =     ( a , b , t , theta ) 
  | and [ 2 * alpha <= pi , 2 * beta <= pi ]   = (a , b , t , theta)  
  | otherwise = if 2 * alpha > pi then findPoints v t xs else findPoints u t xs 
     where   
       t4 (_ ,  _ , _ , t) = t 
       (a , b , t , theta) = minimumBy (compare `on` t4) . map (findAngle u v)  $ xs 
       alpha = t4 $ findAngle v t u  --angle between v u t angle subtended at u by v t
       beta = t4 $ findAngle u t v   -- angle between u v t angle subtended at v by  u t
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.