# -*- coding: utf-8 -*-
import inspect
import logging
import copy
import numpy as np
import numba
import mpmath
import fractalshades as fs
import fractalshades.settings
import fractalshades.numpy_utils.xrange as fsx
logger = logging.getLogger(__name__)
# Log attributes
# https://docs.python.org/2/library/logging.html#logrecord-attributes
[docs]
class Projection:
"""
A Projection defines a mapping acting on the screen-pixels before
iteration.
It implements the associated modifications of the derivatives, when
needed (for distance estimation or shading).
For a `xy_ratio` - zoom parameter - equal to 1, screen pixels are defined
in :math:`[0.5, 0.5] \\times [-0.5, 0.5]`
More generally, the screen is a
rectangle and pixel coordinates take values in:
.. math::
[0.5, 0.5] \\times \\left[ -\\frac{0.5}{xy\\_ratio},
\\frac{0.5}{xy\\_ratio} \\right]
The projection is defined as:
.. math::
(\\bar{x}_{pix}, \\bar{y}_{pix}) = f(x_{pix}, y_{pix})
or alternatively with complex notations:
.. math::
\\bar{z}_{pix} = f(z_{pix})
Derived classes shall implement the actual :math:`f` function
"""
@property
def init_kwargs(self):
""" Return a dict of parameters used during __init__ call"""
init_kwargs = {}
for (
p_name, param
) in inspect.signature(self.__init__).parameters.items():
# By contract, __init__ params shall be stored as instance attrib.
# ('Projection' interface)
init_kwargs[p_name] = getattr(self, p_name)
return init_kwargs
def __eq__(self, other):
"""
2 `Projections` are equal if:
- they are instances from the same subclass
- they have been created with the same init_kwargs
"""
return (
other.__class__ == self.__class__
and other.init_kwargs == self.init_kwargs
)
def adjust_to_zoom(self, fractal):
""" Some projection might impose constraints on the xy_ratio"""
# Default implementation: does nothing
self.make_impl()
def make_impl(self):
self.make_f_impl()
self.make_df_impl()
self.make_dfBS_impl()
def get_impl(self):
""" complex to complex mapping """
raise NotImplementedError("Derived classes shall implement")
def get_df_impl(self):
""" complex to complex mapping """
raise NotImplementedError("Derived classes shall implement")
def get_dfBS_impl(self):
""" complex to 4 floats mapping """
raise NotImplementedError("Derived classes shall implement")
@property
def df_scale(self):
""" Used in to Perturbation fractal as a corrective
scaling factor used in dZndC computation (and consequently distance
estimation, ...)
Projections supporting arbitray precision shall implement.
Returns:
--------
df_scale: numba compiled function pix -> scaling or None
"""
raise NotImplementedError(
f"Arbitray precision not supported by {self.__class__.__name__}"
)
@property
def bounding_box(self, xy_ratio):
""" Passed to Perturbation fractal as an enveloppe of the mapping of
the pixels (before general scaling `scale`)
Projections supporting arbitray precision shall implement.
Returns:
--------
w, h: floats, width and height of the bounding box
"""
# proj_dx = self.dx * self.projection.scale
# corner_a = lin_proj_impl_noscale(0.5 * (w + 1j * h)) * proj_dx
# corner_a is absolute the distance to center
raise NotImplementedError(
f"Arbitray precision not implemented for {self.__class__.__name__}"
)
@property
def min_local_scale(self):
""" Used in to Perturbation fractal to define min pix size. Minimum of
local |df| (before general scaling `scale`)
Projections supporting arbitray precision shall implement.
Returns:
--------
scale: mpmath.mpf or float, the scale factor
"""
# pix = proj.min_local_scale * proj.scale * (f.dx / f.nx)
raise NotImplementedError(
f"Arbitray precision not supported by {self.__class__.__name__}"
)
# Development notes: Projections and derivatives / normals
# For perturbation fractals, the derivatives are relative to a "window" scale
# dx (in Cartesian mode). This scale become dx * projection.scale (constant
# real scalar) for non-cartesian projections.
# On top of this general scaling, 2 "local" modifiers may be implemented
# - total derivative (account for variable local scaling + rotation + skew)
# (mandatory)
# - derivative "without rotation" useful if the final aim is to unwrap, for
# instance to make a videofrom a expmap (optionnal)
#==============================================================================
[docs]
class Cartesian(Projection):
[docs]
def __init__(self, expmap_seam=None):
"""
A Cartesian projection. This is simply the identity function:
.. math::
\\bar{z}_{pix} = z_{pix}
This class can be used with arbitrary-precision deep zooms.
Parameters
==========
`expmap_seam`: None | float
If not None, will ensure a smooth transition between cartesian and
exponential projections (by simulating the scaling a Expmap applies
to the derivatives). The `expmap_seam` value shall usually be set
at 1.0, if the Expmap projection `dx` is half of the Cartesian
projection `dx` (ie, the reference diameter of the Expmap matches
the reference width of the cartesian mapping).
Otherwise, it is an additionnal scaling factor.
Note: parameter only valid for arbitrary-precision implementations.
"""
self.scale = 1.0
self.expmap_seam = expmap_seam
if expmap_seam is not None:
self.make_dzndc_modifier()
def make_f_impl(self):
""" A cartesian projection just let pass-through the coordinates"""
self.f = cartesian_numba_impl
def make_df_impl(self):
self.df = None
def make_dfBS_impl(self):
self.dfBS = None
@property
def df_scale(self):
return None
def bounding_box(self, xy_ratio):
return 1., 1. / xy_ratio
@property
def min_local_scale(self):
return 1.
def make_dzndc_modifier(self):
""" Applies the scaling part of exp mapping at end of iteration.
"""
expmap_seam = self.expmap_seam
# Development notes
# Expmap at h=0 -> kdiff_proj = expmap_seam, pixabs_proj = expmap_seam / 2
# Cartesian -> kdiff_proj = 1., pixabs_proj = 1.
@numba.njit(nogil=True, fastmath=False)
def numba_impl(cpix):
return np.abs(cpix + 1.e-6) * expmap_seam
self.dzndc_modifier = numba_impl
@numba.njit(nogil=True, fastmath=True)
def cartesian_numba_impl(z):
return z
#==============================================================================
[docs]
class Expmap(Projection):
[docs]
def __init__(self, hmin, hmax, rotates_df=True, orientation="horizontal"):
"""
An exponential projection will map :math:`z_{pix}` as follows:
.. math::
\\bar{z}_{pix} = \\exp(h_{moy}) \\cdot \\exp(dh \\cdot z_{pix})
where:
.. math::
h_{moy} &= \\frac{1}{2} \\cdot (h_{min} + h_{max}) \\\\
dh &= h_{max} - h_{min}
This class can be used with arbitrary-precision deep zooms.
Parameters
==========
hmin: str or float or mpmath.mpf
scaling at the lower end of the h-axis, hmin >= 0.
hmax: str or float or mpmath.mpf
scaling at the higher end of the h-axis, hmax > hmin
rotates_df: bool
If ``True`` (default), the derivative will be scaled but also
rotated according to the mapping. If ``False``, only the scaling
will be taken into account.
A rule of thumb is this value shall be set to ``True``
for a standalone picture, and to ``False`` if used as input for a
movie making tool.
orientation: "horizontal" | "vertical"
The direction for the h axis. Defaults to "horizontal".
Notes
=====
.. note::
*Adjustment of zoom parameters:*
The `xy_ratio` of the zoom will be adjusted (at runtime) to
ensure that :math:`\\bar{y}_{pix}` extends from :math:`- \\pi`
to :math:`\\pi`. `nx` is interpreted as `nh` for all values of the
`direction` parameter ("horizontal" or "vertical").
.. note::
*Large mappings and computation by steps:*
For large exponentional mappings which may cover several orders of
magnitude, computation will be run by steps. Each step:
- reuses the same reference orbit
- updates the BLA table as needed by the scale
- updates the reference size used to avoid overflows in
derivatives
"""
if not(0 <= hmin < hmax):
raise ValueError(
"Provide hmin, hmax with: 0 <= hmin < hmax for Expmap"
)
self.use_step = False
self.numba_step_implemented = False
self.rotates_df = rotates_df
self.orientation = orientation
self.premul_1j = {"horizontal": False, "vertical": True}[orientation]
if mpmath.exp(hmax) > (1. / fs.settings.xrange_zoom_level):
# Or ~ hmax > 690... We store internally as Xrange
self.hmin = fsx.mpf_to_Xrange(hmin)
self.hmax = fsx.mpf_to_Xrange(hmax)
else:
# We store internally as float
self.hmin = float(hmin)
self.hmax = float(hmax)
self.hmoy = (hmin + hmax) * 0.5
self.dh = hmax - hmin
def nh(self, fractal):
return fractal.ny if self.premul_1j else fractal.nx
def nt(self, fractal):
return fractal.nx if self.premul_1j else fractal.ny
@property
def pix_to_ht(self):
""" Factor that transforms a pixel to ht coordinates
Note: 0 (scling 1.) is centered
"""
dh = self.dh
premul_1j = self.premul_1j
xy_ratio = self.xy_ratio
return (1j * dh * xy_ratio) if premul_1j else dh
def set_exp_zoom_step(self, exp_step_hmax, exp_step_hmin):
# property defined in ``save_db`` with exp_zoom_step set
logger.debug(
f"set_exp_zoom_step: {exp_step_hmin} {exp_step_hmax}\n"
f"for a full range of {self.hmin} {self.hmax}"
)
self.exp_step_hmax = exp_step_hmax
self.exp_step_hmin = exp_step_hmin
self.exp_step_hmoy = 0.5 * (exp_step_hmax + exp_step_hmin)
self.exp_step_dh = (exp_step_hmax - exp_step_hmin)
self.use_step = True
self.make_dzndc_modifier()
if not(self.numba_step_implemented):
self.make_df_impl()
self.make_dfBS_impl()
def del_exp_zoom_step(self):
""" property used in ``save_db`` with exp_zoom_step set """
self.use_step = False
self.make_df_impl()
self.make_dfBS_impl()
def adjust_to_zoom(self, fractal):
# We need to adjust the fractal xy_ratio in order to match hmax - hmin
# target: dh = 2. * np.pi * xy_ratio
if self.premul_1j:
xy_ratio = (np.pi * 2.) / self.dh
nx = int(fractal.nx * xy_ratio + 0.5)
else:
xy_ratio = self.dh / (np.pi * 2.)
nx = fractal.nx
fractal.xy_ratio = self.xy_ratio = xy_ratio
fractal.zoom_kwargs["xy_ratio"] = xy_ratio
fractal.nx = nx
fractal.zoom_kwargs["nx"] = nx
logger.info(
"Adjusted parameters for Expmap projection:\n"
f"zoom parameter nx: {fractal.nx}\n"
f"zoom parameter xy_ratio: {fractal.xy_ratio}"
)
self.make_impl()
self.make_dzndc_modifier()
def make_f_impl(self):
""" apply the full transform """
hmoy = self.hmoy
pix_to_ht = self.pix_to_ht
@numba.njit(nogil=True, fastmath=False)
def numba_impl(pix):
return np.exp(hmoy + pix_to_ht * pix)
self.f = numba_impl
def make_df_impl(self):
""" Rotate the derivatives of the transform, according to the options
"""
if self.use_step:
return self.make_df_impl_step()
pix_to_ht = self.pix_to_ht
if self.rotates_df:
@numba.njit(nogil=True, fastmath=False)
def numba_impl(pix):
res = np.exp(pix_to_ht * pix)
return res
else:
@numba.njit(nogil=True, fastmath=False)
def numba_impl(pix):
res = np.exp(np.real(pix_to_ht * pix))
return res
self.df = numba_impl
def make_df_impl_step(self):
self.numba_step_implemented = True # Flag to avoid re-compiling
pix_to_ht = self.pix_to_ht
if self.rotates_df:
@numba.njit(nogil=True, fastmath=False)
def numba_impl(pix):
res = np.exp(1j * np.imag(pix_to_ht * pix))
return res
else:
@numba.njit(nogil=True, fastmath=False)
def numba_impl(pix):
return 1.
self.df = numba_impl
def make_dfBS_impl(self):
""" Rotate the derivatives of the transform, according to the options
"""
if self.use_step:
return self.make_dfBS_impl_step()
pix_to_ht = self.pix_to_ht
if self.rotates_df:
@numba.njit(nogil=True, fastmath=False)
def numba_impl(pix):
ht = pix_to_ht * pix
h = np.real(ht)
t = np.imag(ht)
r = np.exp(h)
cr = np.cos(t) * r
sr = np.sin(t) * r
return -sr, -cr, cr, -sr
else:
@numba.njit(nogil=True, fastmath=False)
def numba_impl(pix):
r = np.exp(np.real(pix_to_ht * pix))
return r, 0., 0., r
self.df = numba_impl
def make_dfBS_impl_step(self):
self.numba_step_implemented = True # Flag to avoid re-compiling
pix_to_ht = self.pix_to_ht
if self.rotates_df:
@numba.njit(nogil=True, fastmath=False)
def numba_impl(z):
t = np.imag(pix_to_ht * z)
c = np.cos(t)
s = np.sin(t)
return c, -s, s, c
else:
@numba.njit(nogil=True, fastmath=False)
def numba_impl(z):
return 1., 0., 0., 1.
self.dfBS = numba_impl
def make_dzndc_modifier(self):
""" Applies the scaling part of mapping at end of iteration, so that
it is not needed at postproc stage - used for step iterations.
(Otherwise we would need to store also step information to disk for
image postprocessing stage)
"""
pix_to_ht = self.pix_to_ht
if self.use_step:
hshift = self.hmoy - self.exp_step_hmoy
else:
hshift = self.hmoy
@numba.njit(nogil=True, fastmath=False)
def numba_impl(cpix):
return np.exp(np.real(pix_to_ht * cpix) + hshift)
self.dzndc_modifier = numba_impl
@property
def scale(self):
# Returns the scaling - for perturbation implementation
if self.use_step:
hmoy_df = self.exp_step_hmoy
else:
hmoy_df = self.hmoy
return mpmath.exp(hmoy_df)
def bounding_box(self, xy_ratio):
""" This is the scaling of the bounding box for this calculation step
ie either the whole image or the current exmpap step sub-image.
"""
if self.use_step:
wh = self.exp_step_hmax
else:
wh = self.hmax
w = mpmath.exp(wh)
return w, w
@property
def min_local_scale(self):
""" This is the scaling of the minimum pixel due to the projection
for the *whole* image (accumulated by steps if activated)
"""
return mpmath.exp(self.hmin)
#==============================================================================
[docs]
class Generic_mapping(Projection):
[docs]
def __init__(self, f, df):
"""
This class allows the user to provide a custom mapping defined
on the complex plane.
The transformation between the complex plane and local pixel
coordinates (map :math:`z_{pix}`, according to the zoom level)
is managed internally.
Known limitations:
- currently only differentiable mappings of the
complex variable (aka holomorphic / meromorphic functions) are
supported.
- not suitable for arbitray precision fractals.
Parameters
==========
f: numba jitted function numba.complex128 -> numba.complex128
The complex function defining the mapping. If a tuple (f1, f2) or
(f1, f2, f3) is provided the composition will be applied (f1 o f2)
df: numba jitted function numba.complex128 -> numba.complex128
The differential of f. If a tuple (df1, df2) or (df1, df2, df3) is
provided the differential of the composition will be applied
according to the differentiation chain rule.
"""
self.f = f
self.df = df
_f = f if fs.utils.is_iterable(f) else (f,)
_df = df if fs.utils.is_iterable(df) else (df,)
if (len(_df) != len(_f)):
raise ValueError("len(f) and len(df) shall match")
self.n_func = len(_f)
self.P_f_numba = copy.copy(_f)
self.P_df_numba = copy.copy(_df)
def make_phi(self):
"""
phi : pixel to C-plane
phi_inv: C-plane to pixel
"""
z_center = self.z_center
dx = self.dx
@numba.njit(nogil=True, fastmath=False)
def numba_phi_impl(zpix):
return z_center + zpix * dx
@numba.njit(nogil=True, fastmath=False)
def numba_phi_inv_impl(z):
return (z - z_center) / dx
self.phi = numba_phi_impl
self.phi_inv = numba_phi_inv_impl
def adjust_to_zoom(self, fractal):
""" JIT compilation is delayed until we know the zoom
parameters"""
self.z_center = complex(fractal.x, fractal.y)
self.dx = fractal.dx
self.make_phi()
self.make_impl()
def make_f_impl(self):
n_func = self.n_func
phi = self.phi
phi_inv = self.phi_inv
P_f_numba = self.P_f_numba
if n_func == 1:
f0 = P_f_numba[0]
@numba.njit(numba.complex128(numba.complex128),
nogil=True, fastmath=False)
def numba_impl(z):
zz = phi(z) # pix -> P
zz = f0(zz)
zz = phi_inv(zz) # P -> pix
return zz
elif n_func == 2:
f0 = P_f_numba[0]
f1 = P_f_numba[1]
@numba.njit(numba.complex128(numba.complex128),
nogil=True, fastmath=False)
def numba_impl(z):
zz = phi(z) # pix -> P
zz = f0(zz)
zz = f1(zz)
zz = phi_inv(zz) # P -> pix
return zz
elif n_func == 3:
f0 = P_f_numba[0]
f1 = P_f_numba[1]
f2 = P_f_numba[2]
@numba.njit(numba.complex128(numba.complex128),
nogil=True, fastmath=False)
def numba_impl(z):
zz = phi(z) # pix -> P
zz = f0(zz)
zz = f1(zz)
zz = f2(zz)
zz = phi_inv(zz) # P -> pix
return zz
else:
raise NotImplementedError("Composition allowed up to 3 funcs only")
self.f = numba_impl
def make_df_impl(self):
n_func = self.n_func
phi = self.phi
P_f_numba = self.P_f_numba
P_df_numba = self.P_df_numba
if n_func == 1:
df0 = P_df_numba[0]
@numba.njit(numba.complex128(numba.complex128),
nogil=True, fastmath=False)
def numba_impl(z):
zz = phi(z) # pix -> P
dzz = df0(zz)
return dzz
elif n_func == 2:
df0 = P_df_numba[0]
f0 = P_f_numba[0]
df1 = P_df_numba[1]
@numba.njit(numba.complex128(numba.complex128),
nogil=True, fastmath=False)
def numba_impl(z):
zz = phi(z) # pix -> P
dzz = df0(zz)
zz = f0(zz)
dzz *= df1(zz)
return dzz
elif n_func == 3:
df0 = P_df_numba[0]
f0 = P_f_numba[0]
df1 = P_df_numba[1]
f1 = P_f_numba[1]
df2 = P_df_numba[2]
@numba.njit(numba.complex128(numba.complex128),
nogil=True, fastmath=False)
def numba_impl(z):
zz = phi(z) # pix -> P
dzz = df0(zz)
zz = f0(zz)
dzz *= df1(zz)
zz = f1(zz)
dzz *= df2(zz)
return dzz
else:
raise NotImplementedError("Composition allowed up to 3 funcs only")
self.df = numba_impl
def make_dfBS_impl(self):
""" Currently not implemented, will raise an assertion error if
called in LLVM """
@numba.njit(nogil=True, fastmath=False)
def numba_impl(z):
assert False
return z
self.dfBS = numba_impl