I run into precision problems very easily
Thanks. Indeed with 24 bit float you can't do much. Another example:
But! You can use the derivative of the inflections (before Julia iterations) to calculate the required precision needed
Something like
max(16, ceil(16 - log2(abs(derivative))))
Using this precision (different for each pixel) gives a clean image:
Here's a full working Haskell program using the 'rounded' package for MPFR bindings:
{-
Simple inflection mapping with dynamic precision. First compute the inflections
with low constant precision, then use the computed derivative to decide how much
precision to use for each pixel in the final render.
Usage:
ghc -O2 Main.hs
./Main > output.pgm
Runtime is about 5m30s on one core of a 3GHz amd64 desktop running Debian.
-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}
module Main (main) where
-- MPFR binding, not on Hackage yet, get it from either of these:
-- https://code.mathr.co.uk/rounded/shortlog/refs/heads/claude5
-- https://github.com/ekmett/rounded/pull/27
import Numeric.Rounded
import Data.Complex (Complex((:+)), magnitude)
import Data.List (foldl')
import Data.Monoid ((<>))
import Data.Proxy (Proxy)
import Data.Word (Word8)
import qualified Data.ByteString as BS
import qualified Data.ByteString.Char8 as BSC
-- arbitrary precision types (p bits of precision)
type R (p :: k) = Rounded 'TowardNearest p
type C (p :: k) = Complex (R p)
-- dual numbers for automatic differentiation
data AD p = AD !(C p) !(C p)
deriving (Eq, Read, Show)
instance Precision p => Num (AD p) where
negate (AD a b) = AD (negate a) (negate b)
AD a b + AD c d = AD (a + c) (b + d)
AD a b - AD c d = AD (a - c) (b - d)
AD a b * AD c d = AD (a * c) (a * d + b * c)
fromInteger = constant . fromInteger
constant :: Precision p => C p -> AD p
constant p = AD p 0
variable :: Precision p => C p -> C p -> AD p
variable d p = AD p d
getVariable :: Precision p => AD p -> C p
getVariable (AD x _) = x
getDerivative :: Precision p => AD p -> C p
getDerivative (AD _ y) = y
-- inflection mapping
inflect :: Precision p => C p -> AD p -> AD p
inflect p c =
let f = constant p
d = c - f
in d * d + f
inflects :: Precision p => [C p] -> AD p -> AD p
inflects ps z = foldl' (flip inflect) z ps
-- estimate required precision from derivative
precisionRequired :: Precision p => AD p -> Int
precisionRequired (AD _ d) = max 16 . ceiling $ 16 - logBase 2 (magnitude d)
-- generate a grid of coordinates
grid :: Precision p => Int -> Int -> C p -> R p -> [[C p]]
grid w h c r =
[ [ c + d * (x i :+ y j) | i <- take w [0..] ] | j <- take h [0..] ]
where
d = 2 * r :+ 0
aspect = fromInt w / fromInt h
x i = ((fromInt i + 0.5) / fromInt w - 0.5) * aspect
y j = (fromInt h - fromInt j + 0.5) / fromInt h - 0.5
clamp :: Ord a => a -> a -> a -> a
clamp mi ma = min ma . max mi
-- simple greyscale colouring
grey :: Double -> Word8
grey = fromIntegral . clamp 0 (255 :: Int) . floor . (255 *) . tanh . clamp 0 8
-- distance estimator
de :: Precision p => AD p -> Double
de (AD z' dz')
| isInfinite mdz || isNaN mdz || isInfinite mz || isNaN mz = 1e10
| otherwise = 2 * mz * log mz / mdz
where
z = fmap toDouble z'
dz = fmap toDouble dz'
mz = magnitude z
mdz = magnitude dz
-- Julia iteration loop
calculate :: Precision p => AD p -> AD p -> AD p
calculate c
= head . (++ [AD 0 0])
. dropWhile ((< escapeRadius) . magnitude . getVariable)
. take maxIters
. iterate (z -> z * z + c)
-- calculate and colour a point
render :: Precision p => [C p] -> AD p -> Proxy p -> Word8
render ps z0 _proxy = grey . de . calculate c . inflects ps $ z0
where
c = constant . fmap precRound . last $ ps
-- compute a pixel
pixel :: Precision p => [C p] -> C p -> Word8
pixel ps z0 = reifyPrecision prec render'
where
-- compute precision at base precision of input
-- prec = 24 -- uncomment this and comment the next line to see the effect
prec = precisionRequired . inflects ps . variable delta0 $ z0
delta0 = fmap precRound delta
-- round to computed necessary precision for rendering
render' :: forall p . Precision p => Proxy p -> Word8
render' = render ps' (variable delta' z0')
where
z0' = fmap precRound z0
delta' = fmap precRound delta
ps' = map (fmap precRound) ps
-- pixel size
delta :: C 24
delta = radius * 2 / fromInt height :+ 0
main' :: [C 24] -> BS.ByteString
main' ps
= BSC.pack ("P5
" ++ show width ++ " " ++ show height ++ "
255
")
<> BS.pack (concatMap (map (pixel ps)) (grid width height (head ps) radius))
main :: IO ()
main = BS.putStr (main' example)
-- image parameters
width :: Int
width = 1280
height :: Int
height = 702
-- iteration parameters
escapeRadius :: Precision p => R p
escapeRadius = 65536
maxIters :: Int
maxIters = 1000
-- view size
radius :: R 24
radius = 1.1
-- example inflections
example :: [C 24]
example =
[ 0.29166666 :+ 1.20277774
, 0.15833333 :+ 1.21944439
, 0.06388889 :+ 1.20277774
, (-0.05833333) :+ 1.21388888
, (-0.22500001) :+ 1.17499995
, (-0.34166667) :+ 1.16388893
, (-0.48055553) :+ 1.19166672
, (-0.59722221) :+ 1.13611114
, (-0.69166666) :+ 1.06388891
, (-0.80833334) :+ 0.96388888
, (-0.76944447) :+ 0.13055556
]
You can download the program here:
https://mathr.co.uk/misc/2017-03-21_inflection-mapping-dynamic-precision.hsIt's not really practical for live exploring though, 24bit float on GPU is pretty much instant, while this variable precision software floating point on CPU takes several minutes at the same resolution - probably it would be beneficial to just use double precision, or double-double if required for some pixels...