Julia Sets

%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
import datajoint as dj

def julia(c, size=256, center=(0.0, 0.0), zoom=1.0, iters=256):
    x, y = np.meshgrid(
        np.linspace(-1, 1, size)/zoom + center[0], 
        np.linspace(-1, 1, size)/zoom + center[1], 
    )
    z = x + 1j * y
    im = np.zeros(z.shape)
    ix = np.ones(z.shape, dtype=bool)
    for i in range(iters):
        z[ix] = z[ix] ** 2 + c
        ix = np.abs(z) < 2
        im += ix
    return im

plt.imshow(julia(-0.4+0.6j), cmap='magma')
plt.axis(False);

plt.imshow(julia(-0.4+0.6j, center=(0.34, -0.30), zoom=10000.0), cmap='magma')
plt.axis(False);

c = (
    -0.4 + 0.6j, 
    -0.74543 + 0.11301j, 
    -0.75 + 0.11j, 
    -0.1 + 0.651j,
    -0.835 - 0.2321j,
    -0.70176 - 0.3842j,
)

noise_level = 5.0

fig, ax = plt.subplots(3, 2, figsize=(7.5, 12))
for c_, a in zip(c, ax.flatten()):
    img = julia(c_, zoom=0.5) 
    img += np.random.randn(*img.shape) * noise_level
    a.imshow(img, cmap='magma')
    a.axis(False)

from skimage import data
from skimage import filters

from skimage.morphology import disk
from skimage import restoration

noise_level = 50.0
img = julia(-0.4+0.6j, size=200)
noise_img = img + np.random.randn(*img.shape) * noise_level
median_img = filters.median(noise_img, disk(3))
tv_img = restoration.denoise_tv_chambolle(noise_img, weight=20.0)
wavelet_img = restoration.denoise_wavelet(noise_img)
gaussian_img = filters.gaussian(noise_img, sigma=1.8)

fig, ax = plt.subplots(3, 2, figsize=(6, 9))
for a, (im, title) in zip(
    ax.flatten(),
    ((img, 'original'), 
     (noise_img, 'original+noise'),
     (gaussian_img, 'gaussian'),
     (median_img, 'median'), 
     (wavelet_img, 'wavelet'),
     (tv_img, 'tv'), )):
    a.imshow(im, cmap='magma', vmin=0, vmax=255)
    a.axis(False)
    a.set_title(title)

import datajoint as dj

schema = dj.Schema('julia')

@schema 
class JuliaSpec(dj.Lookup):
    definition = """
    julia_spec : smallint 
    ---
    creal : float
    cimag : float
    size=256 : smallint
    center_real=0.0 : float 
    center_imag=0.0 : float
    zoom=1.0 : float
    noise_level=50 : float
    """  

    contents =  (
        dict(julia_spec=0, creal=-0.4, cimag=0.6, noise_level=50),
        dict(julia_spec=1, creal=-0.7453, cimag=0.11301, noise_level=50),
        dict(julia_spec=2, creal=-0.75, cimag=0.11, noise_level=50),
        dict(julia_spec=3, creal=-0.1, cimag=0.651, noise_level=50),
        dict(julia_spec=4, creal=-0.835, cimag=-0.2321, noise_level=50),
        dict(julia_spec=5, creal=-0.70176, cimag=-0.3842, noise_level=50),
    )

JuliaSpec.insert1(
    dict(julia_spec=10, 
         creal=-0.4, cimag=0.6, 
         center_real=0.34, center_imag=-0.30, 
         zoom=10000.0, 
         noise_level=50.0)
)  

@schema
class JuliaImage(dj.Computed):
    definition = """
    -> JuliaSpec 
    ---
    image : longblob
    """

    def make(self, key):
        spec = (JuliaSpec & key).fetch1()
        img = julia(spec['creal'] + 1j*spec['cimag'], 
                    size=spec['size'],
                    center=(spec['center_real'], spec['center_imag']),
                    zoom=spec['zoom'],
                    )
        img += np.random.randn(*img.shape) * spec['noise_level']
        self.insert1(dict(key, image=img.astype(np.float32)))

dj.Diagram(schema)

JuliaImage.populate(display_progress=True)

JuliaImage()

plt.imshow((JuliaImage & 'julia_spec=2').fetch1('image'))
plt.axis(False);

dj.Diagram(schema)

@schema
class DenoiseMethod(dj.Lookup):
    definition = """
    denoise_method : smallint
    ---
    method : varchar(16)
    params=null : blob
    """
    contents = (
        (0, 'gaussian', dict(sigma=1.8)),
        (1, 'median', dict(radius=3)),
        (2, 'wavelet', {}),
        (3, 'tv', dict(weight=20.0))
    )

@schema 
class JuliaDenoised(dj.Computed):
    definition = """
    -> JuliaImage
    -> DenoiseMethod
    ---
    denoised_image : longblob
    """

    def make(self, key):
        img = (JuliaImage & key).fetch1('image')
        method, params = (DenoiseMethod & key).fetch1('method', 'params')

        if method == "gaussian":
            result = filters.gaussian(img, **params)
        elif method == "median":
            result = filters.median(img, disk(params['radius']))
        elif method == 'tv':
            result = restoration.denoise_tv_chambolle(img, **params)
        elif method == "wavelet":
            result = restoration.denoise_wavelet(noise_img, **params)
        else:
            raise NotImplementedError
        self.insert1(dict(key, denoised_image=result))

JuliaDenoised.populate(display_progress=True)

JuliaDenoised()

keys = JuliaDenoised.fetch('KEY')
img = ((JuliaDenoised & keys[21])).fetch1('denoised_image')
plt.imshow(img)
plt.axis(False);

dj.Diagram(schema)

@schema
class Blob(dj.Manual):
    definition = """
    id  : int
    ---
    blob : longblob
    """

Blob.insert1(dict(id=1, blob=[1, 2, 3, 'Four']))

Blob.fetch(as_dict=True)