2014-11-16 15 views
13

Zaimplementowałem W3s recommended algorithm for converting SVG-path arcs from endpoint-arcs to center-arcs and back w Haskell.Dlaczego moja implementacja konwersji łuku SVG nie przechodzi testu QuickCheck?

type EndpointArc = (Double, Double, Double, Double 
        , Bool, Bool, Double, Double, Double) 

type CenterArc = (Double, Double, Double, Double 
       , Double, Double, Double) 

endpointToCenter :: EndpointArc -> CenterArc 

centerToEndpoint :: CenterArc -> EndpointArc 

See full implementation and test-code here.

Ale nie mogę dostać tę nieruchomość przekazać:

import Test.QuickCheck 
import Data.AEq ((~==)) 

instance Arbitrary EndpointArc where 
    arbitrary = do 
     ((x1,y1),(x2,y2)) <- arbitrary `suchThat` (\(u,v) -> u /= v) 
     rx    <- arbitrary `suchThat` (>0) 
     ry    <- arbitrary `suchThat` (>0) 
     phi    <- choose (0,2*pi) 
     (fA,fS)   <- arbitrary 
     return $ correctRadiiSize (x1, y1, x2, y2, fA, fS, rx, ry, phi) 

prop_conversionRetains :: EndpointArc -> Bool 
prop_conversionRetains earc = 
    let result = centerToEndpoint (endpointToCenter earc) 
    in earc ~== result 

Czasami wynika to pływający punkt (błędy, które wydają się przekraczać IEEE 754), ale czasami są Koncepcja nieliczby w wyniku.

(NaN,NaN,NaN,NaN,False,False,1.0314334509082723,2.732814841776921,1.2776112657142984) 

Który wskazuje nie ma rozwiązania, chociaż myślę, że skala Rx, Ry, jak opisano w F.6.6.2 in W3's document.

import Numeric.Matrix 

m :: [[Double]] -> Matrix Double 
m = fromList 

toTuple :: Matrix Double -> (Double, Double) 
toTuple = (\[[x],[y]] -> (x,y)) . toList 

primed :: Double -> Double -> Double -> Double -> Double 
     -> (Double, Double) 
primed x1 y1 x2 y2 phi = toTuple $ 
    m [[ cos phi, sin phi] 
     ,[-sin phi, cos phi] 
     ] 
    * m [[(x1 - x2)/2] 
     ,[(y1 - y2)/2] 
     ] 

correctRadiiSize :: EndpointArc -> EndpointArc 
correctRadiiSize (x1, y1, x2, y2, fA, fS, rx, ry, phi) = 
    let (x1',y1') = primed x1 y1 x2 y2 phi 
     lambda = (x1'^2/rx^2) + (y1'^2/ry^2) 
     (rx',ry') | lambda <= 1 = (rx, ry) 
        | otherwise = ((sqrt lambda) * rx, (sqrt lambda) * ry) 
    in (x1, y1, x2, y2, fA, fS, rx', ry', phi) 

Odpowiedz

13

OK, sam to wymyśliłem. Kluczem był oczywiście dokument W3s:

W przypadku, gdy promienie są skalowane za pomocą równania (F.6.6.3), radicand (F.6.5.2) wynosi zero i jest dokładnie jedno rozwiązanie dla środka elipsy.

F.6.5.2 w moim kod jest

(cx',cy') = (sq * rx * y1'/ry, sq * (-ry) * x1'/rx) 
       where sq = negateIf (fA == fS) $ sqrt 
         $ (rx^2 * ry^2 - rx^2 * y1'^2 - ry^2 * x1'^2) 
         /(rx^2 * y1'^2 + ry^2 * x1'^2) 

radicand, że odnosi się to

(rx^2 * ry^2 - rx^2 * y1'^2 - ry^2 * x1'^2) 
/(rx^2 * y1'^2 + ry^2 * x1'^2) 

Ale oczywiście, ponieważ pracujemy z pływaków nie jest to dokładnie zero, ale w przybliżeniu, a czasami może to być coś w rodzaju -6.99496644301622e-17, co jest negatywne! Pierwiastek kwadratowy liczby ujemnej jest liczbą zespoloną, więc obliczenia zwracają wartość NaN.

Sztuczka polegała na propagowaniu faktu, że wartości rx i ry zostały przeskalowane, aby zwrócić zero i sprawić, że sq jest zerowe, zamiast przechodzić całe obliczenia bez potrzeby, ale szybkie ustalenie jest po prostu wartością bezwzględną radicand.

(cx',cy') = (sq * rx * y1'/ry, sq * (-ry) * x1'/rx) 
       where sq = negateIf (fA == fS) $ sqrt $ abs 
         $ (rx^2 * ry^2 - rx^2 * y1'^2 - ry^2 * x1'^2) 
         /(rx^2 * y1'^2 + ry^2 * x1'^2) 

Po tym są pewne pozostałe problemy zmiennoprzecinkowe. Po pierwsze błąd przekracza co jest dozwolone przez ~== operator IEEE 754, tak zrobiłem mój własny approxEq

approxEq (x1a, y1a, x2a, y2a, fAa, fSa, rxa, rya, phia) (x1b, y1b, x2b, y2b, fAb, fSb, rxb, ryb, phib) = 
     abs (x1a - x1b ) < 0.001 
    && abs (y1a - y1b ) < 0.001 
    && abs (x2a - x2b ) < 0.001 
    && abs (y2a - y2b ) < 0.001 
    && abs (y2a - y2b ) < 0.001 
    && abs (rxa - rxb ) < 0.001 
    && abs (rya - ryb ) < 0.001 
    && abs (phia - phib) < 0.001 
    && fAa == fAb 
    && fSa == fSb 

prop_conversionRetains :: EndpointArc -> Bool 
prop_conversionRetains earc = 
    let result = centerToEndpoint (trace ("FIRST:" ++ show (endpointToCenter earc)) (endpointToCenter earc)) 
    in earc `approxEq` trace ("SECOND:" ++ show result) result 

który rozpoczyna przynosząc przypadki fA jest uzyskiwanie odwróconych. Spot liczbę Magic:

PIERWSZA: (- +5,988957688551294, -39.5430169665332,64.95929681921707,29.661347617532357,5.939852349879405, -1,2436798376040206, 3,141592653589793)

DRUGI (+4,209851895761209, -73,01839718538467, -16,18776727286379, -6,067636747681732 , Fałsz, prawda, 64.95929681921707,29.661347617532357,5.939852349879405)

*** Nie powiodło się! Falsyfikowalne (po 20 testach):
(4.209851895761204, -73.01839718538467, -16.18776781572145, -6.0676366434916655, prawda, prawda, 64.95929681921707,29.661347617532357,5.939852349879405)

Załatwione! fA = abs dtheta > pi jest w centerToEndpoint, więc jeśli jest to therabouts, to może pójść w obie strony.

Więc wyjąłem warunek Fa i zwiększyła liczbę testów w QuickCheck

approxEq (x1a, y1a, x2a, y2a, fAa, fSa, rxa, rya, phia) (x1b, y1b, x2b, y2b, fAb, fSb, rxb, ryb, phib) = 
     abs (x1a - x1b ) < 0.001 
    && abs (y1a - y1b ) < 0.001 
    && abs (x2a - x2b ) < 0.001 
    && abs (y2a - y2b ) < 0.001 
    && abs (y2a - y2b ) < 0.001 
    && abs (rxa - rxb ) < 0.001 
    && abs (rya - ryb ) < 0.001 
    && abs (phia - phib) < 0.001 
    -- && fAa == fAb 
    && fSa == fSb 

main = quickCheckWith stdArgs {maxSuccess = 50000} prop_conversionRetains 

co pokazuje, że approxEq próg nie jest jeszcze na tyle luźne.

approxEq (x1a, y1a, x2a, y2a, fAa, fSa, rxa, rya, phia) (x1b, y1b, x2b, y2b, fAb, fSb, rxb, ryb, phib) = 
     abs (x1a - x1b ) < 1 
    && abs (y1a - y1b ) < 1 
    && abs (x2a - x2b ) < 1 
    && abs (y2a - y2b ) < 1 
    && abs (y2a - y2b ) < 1 
    && abs (rxa - rxb ) < 1 
    && abs (rya - ryb ) < 1 
    && abs (phia - phib) < 1 
    -- && fAa == fAb 
    && fSa == fSb 

Które mogę w końcu przejść pomyślnie niezawodnie z dużą liczbą testów. Cóż, wszystko to po to, żeby zrobić zabawną grafikę ... Jestem pewien, że jest wystarczająco dokładna :)

Powiązane problemy