We study applications of 3D printing to the broad goal of understanding how mathematical objects vary continuously in families.  To do so, we model the varying parameter as the vertical axis of a 3D print, introducing the notion of a static animation: a 3D printed object each of whose layers is a member of the continuously deforming family.  We survey examples and draw connections to algebraic geometry, complex dynamics, chaos theory, and more.  We also include a detailed tutorial (with accompanying code and files) so that the reader can create static animations of their own.

§ INTRODUCTION

Across mathematical disciplines there is considerable interest in understanding how a mathematical object can deform according to a shift in a single parameter, producing a varying family of related mathematical objects.  If these objects have small enough dimension, one can vary the parameter over time and display how the object changes as an animation on a screen. With a 3D printer we can replace the time dimension with the vertical $z$-axis, and stack each successively deformed layer on top of the last. This turns the entire parametrized family of objects into a single solid which we can hold in our hands, giving an entire new physical dimension to the deformation.

The study of continuously varying families is ubiquitous in mathematics.  To name a few examples:  In algebraic geometry, one studies degenerations of curves, surfaces, and higher dimensional algebraic varieties, putting them in continuously varying flat families which interpolate between one space and the next <cit.>.  In complex dynamics, one asks about limiting behavior when iterating a holomorphic function $f$, and compares this to the limiting behavior of a slight perturbation of $f$, to dazzling effects [17].  In topology, the question of whether one space can continuously deform to another forms the notion of a homotopy, defining the fundamental invariant in the rich and expansive field of homotopy theory <cit.>.  In the study of chaotic dynamical systems, one explores how slight changes in initial conditions can radically affect longer term outcomes [24].  In fact, one would be hard pressed to find a field of mathematics which does not contain a multitude of problems that aim to understand how a mathematical object or system varies when perturbed.

In this article we propose a general framework which can be used to explore situations like those described in the previous paragraph.  Namely, we describe ways to unify an entire family of continuously varying objects into a single mathematical object or deformation space.  When the members of this family are 1 or 2 dimensional, we argue that 3D printing this deformation space provides new perspectives on the illustration of continuously deforming geometry, beyond that which an animation on a screen can achieve.  We call these 3D models static animations.

There are two main components to this article, one more theoretical and one more interactive.  The theoretical component precisely formulates the notion of a static animation and describes connections to various mathematical areas, showcasing examples and suggesting avenues for exploration.  The interactive component provides detailed instructions for the reader to create static animations of their own, and is accompanied by code and files in a public github repository [5].

The structure of the article is therefore as follows.  In Section <ref> we give a rigorous definition of the notion of a static animation, and illustrate the concept with an example of a varying family of polar flowers.  In Section <ref> we introduce the author's first 3D printed deformation space, which illustrates a family of Pythagorean tree fractals.  In Section <ref> we enter the world of complex dynamics, showcasing various prints of deformations of Julia sets, including some our own work, as well as work Bernat Espigulé, Caroline Davis, and Bethany Mussman.

We then arrive at the interactive portion in Section <ref>, describing a procedure to create static animations that was developed with Bernat Espigulé.  This workflow can be understood as a sort of reverse tomography.  Tomography is a way to visualize the internal structure of 3D object by producing a collection 2D pictures of cross sections at various depths.  For example, an MRI is a type of tomography which can be applied to the human body.  The process we describe goes in the other direction, starting with a collection of 2D images and stitching them together into a 3D object.

The upshot of this particular workflow is that the user can create static animations using mostly 2D modelling methods$-$minimizing the need for more difficult 3D modeling software.  The main idea is to feed a stack of 2D images to an open source medical and molecular visualization software called Chimera [21],[25] which will stitch them together into a 3D model.[As far as the author is aware, this functionality was initially developed to reconstruct 3D models from medical tomographical images such as MRI images.]  After broadly describing the workflow, we also include a detailed step-by-step tutorial allowing the reader to follow along and create a static animation of their own.  Accompanying the tutorial we have our own code and files from each step of the process, including a file of the finalized 3D model.  It is the hope of the author that after working through this tutorial, any reader able to make 2D animations on a screen will also be able to make 3D static animations.

The idea of viewing a continuously deforming family of a geometric objects as a geometric object in its own right is certainly not new, and in fact it was the algebrogeometric perspective of viewing so-called flat morphisms as parameter spaces which motivated the author's exploration into this field.  We provide more details and an example of this connection in the appendix.

§ ACKNOWLEDGMENTS

The author's initial exploration into static animations was as the fabrication manager of the Washington Experimental Math Lab at the University of Washington$-$supported by NSF grant DMS-1559860$-$and the author thanks Jayadev Athreya for the technology and support in exploring math from an experimental and illustrative perspective during this time.  The author thanks Bernat Espigulé for his generosity and collaboration in creating static animations of Julia sets, Caroline Davis for extremely enlightening conversations, and Edmund Harriss for taking many of the photographs of their models which are included in this manuscript.  We also thank the anonymous referee for many helpful comments and contributions, including some of the python code in our tutorial section, as well as the description of our process as a reverse tomography.  Much of the development of the ideas and models were made during the 2019 ICERM semester titled Illustrating Mathematics$-$supported by NSF grant DMS-1439786$-$and the author thanks all the participants for their collaboration and support during that semester.  The author also thanks the participants of the 2021 Illustrating Mathematics Summer School at PCMI, and the 2022 Joint Mathematics Meetings mini-course 3D Printing: Challenges and Applications.  Many of the 3D models were made with UCSF Chimera, developed by the Resource for Biocomputing, Visualization, and Informatics at the University of California, San Francisco, with support from NIH P41-GM103311.  The author is currently supported by NSF grant DMS-1646385 as part of the Research Training Group in Arithmetic Geometry at UC Berkeley.

§ STATIC ANIMATIONS

In what follows we give a precise mathematical formulation of what we mean when referring to a static animation.

Suppose we are given some interval $I\subseteq\mathbb{R}$, and for each $t\in I$ we are given a subset $B_t\subseteq\mathbb{R}^2$.  The static animation associated to the data $(I,(B_t)_{t\in I})$ is the following subset of $\bR^3$.

\[\cB = \{(b,t):t\in I,b\in B_t\}\subseteq\bR^3.\]

The sets $B_t$ are called the frames of the static animation.

This gives us some subset of $\bR^3$, which we could hope to 3D print.[That said, if the set $B_t$ varies too wildly with $t$, we will probably be unable to do so.  To have any hope, we'd like $B_t$ to vary continuously in $t$.  Loosely speaking, this means that if $t\approx t'$ then $B_t\approx B_{t'}$.  This can be made precise in a number of ways, depending on the context.]  Observe that a static animation comes equipped with a projection map $\pi:\cB\to I$ given by the rule $\pi(b,t) = t$.  For any $t_0\in I$ we can recover the frame $B_{t_0}$ as the the fiber of $\pi$ over $t_0$:

\[\pi^{-1}(t_0) = \{(b,t_0)\in\cB\} \cong B_{t_0}.\]

We illustrate the notion of a static animation and its associated projection map with an example of a continuously varying family of curves, interpolating between the polar flowers $B_{t_1},B_{t_2},$ and $B_{t_3}$, below.

[domain=0:2*pi,scale=1,samples=500] plot (deg():2 + cos(5 * (r)));
at (-3,2) $B_{t_1}$;
[domain=0:2*pi,scale=1,samples=500] plot (deg():2 + .2*cos(5 * (r) + 100));
at (-2,2) $B_{t_3}$;

[domain=0:2*pi,scale=1,samples=500] plot (deg():2 + .6*cos(5 * (r ) + 50));
at (-2.5,2) $B_{t_2}$;

We start by aligning these curves vertically, with the polar curve $B_t$ curve living on the plane $z=t$.

[hide axis]
3[variable=t,domain=0:360,samples=65] ((2 + cos(5 * (t) ))*sin(t),(2 + cos(5 * (t)))*cos(t), 0);
3[variable=t,domain=0:360,samples=65] ((2 + .6*cos(5 * (t) + 50))*sin(t),(2 + .6*cos(5 * (t)+50))*cos(t), .4);
3[variable=t,domain=0:360,samples=65] ((2 + .2*cos(5 * (t)+100))*sin(t),(2 + .2*cos(5 * (t)+100))*cos(t), .8);

at (-.5,4.25) $B_{t_3}$;
at (-.5,2.75) $B_{t_2}$;
at (-.5,1.25) $B_{t_1}$;

Then we interpolate by filling in the rest of the $t$ values to get a surface, $\cB$.  We will explicitly describe this surface in Section <ref> below.

[hide axis]
3[surf,samples=50,samples y =30, domain=0:360,domain y=0.2:1] ((2 + y*cos(5*(x)+  125*(1-y)))*sin(x),(2 + y*cos(5*(x)+125*(1-y)))*cos(x), 1-y);

at (-.5,2.5) $\cB$;

The projection $\pi:\cB\to\bR$ consists of simply plucking out the $z$-coordinate of a point on this surface.  For $t\in I$, the fiber over $t$ of this projection, $\pi^{-1}(t)$, is therefore the horizontal cross section given by intersecting $\cB$ with the plane $z=t$.  By construction, this is precisely the polar flower $B_t$, floating at height $t$.

We can now cap this surface to make a solid which can be 3D printed as a deformation space illustrating how to continuously interpolate between polar flowers.

[This is a screenshot of the STL file opened and rendered in Ultimaker-Cura [4]]

We give a step by step tutorial for the creation of this object in Section <ref>.  An STL file of the 3D model is also available at our public github repository [5], and can be downloaded to print.

Of course, one could also illustrate the continuous deformation of polar flowers as an animation on a screen.[We made such an animation.  Find it here: <http://www.gabrieldorfsmanhopkins.com/digital_fabrication/polarFlowerDeformation/flowerAnimation.gif>]  Indeed, the horizontal cross sections of the solid described above correspond precisely to the frames of this animation.  This is why we call these objects static animations, and why we call their horizontal cross sections frames.

§.§ A Mathematical Description of The Deformation of Polar Flowers

We continue with a more precise formulation of the example above.  Adopting the same notation, $B_{t_1}$ is the polar flower given by the equation

\[r = 2 + cos(5\theta)\]

in polar coordinates.  Similarly, $B_{t_3}$ is given by

\[r = 2 + \frac{1}{5}\cos(5\theta + 2\pi/5).\]

In particular, both consist of waves on the circle of radius 2, although they differ in amplitude and a shift of where the wave begins.  These correspond to two varying parameters in the equation: the amplitude is controlled by the coefficient in front of the cosine, and the shift corresponds to adding a constant term to the argument of the cosine.  These can be adjusted simultaneously to interpolate continuously  between $B_{t_1}$ and $B_{t_3}$, giving a family of curves whose general element $B_t$ is given by the equation

\[r = 2 + t\cos(5\theta + 2\pi t).\]

Indeed, letting t=1 we recover $B_{t_1}$ and letting $t=1/5$ we recover $B_{t_3}.$[We can also recover $B_{t_2}$ by letting $t=3/5$.]  From this we derive an explicit parametrization for our deformation surface

\[\rho:[0,2\pi]\times[.2,1]\longrightarrow\bR^3.\]

Giving target cylindrical coordinates, $\rho$ can be expressed by the rule

\[\rho(\theta,t) = (2+t\cos(5\theta+2\pi t),\theta,t).\]

The static animation $\cB$ is precisely the image of $\rho$, and by restricting to $\cB$ the projection $\bR^3\to\bR$ onto the $z$-axis, we recover the projection $\pi:\cB\to I$.  For $t_0\in I$ we compute that

\[\pi^{-1}(t_0) = \{\rho(\theta,t_0):\theta\in[0,2\pi]\} = \{(2+t_0\cos(5\theta+2\pi t_0),\theta,t_0):\theta\in[0,2\pi]\}.\]

This is precisely the polar flower given by $r = 2+t_0\cos(t\theta+2\pi t_0)-$that is, the frame $B_{t_0}-$floating at the fixed height $t_0$.

§.§ Higher Dimensional Stacks of Frames

It is worth returning to the observation that there were two varying parameters in our family of polar flowers, one changing amplitude and another controlling the phase shift.  We controlled both using the single parameter $t$, but of course each can move independantly, say, with parameters $s$ and $t$.  This would give a 2-dimensional parameter space of polar flowers, say $\bR^2$, and for every $(s,t)\in\bR^2$ we could consider the curve

\begin{equation}\label{2DFrame}
C_{s,t}:\hspace{10pt} r = 2 + s\cos(5\theta + 2\pi t).
\end{equation}

This puts us in a situation similar to Definition <ref>, except that we have a 2-dimensional stack of frames.  Therefore, trying to build the static animation would put us in 4-space,

\[{\cC} = \{(b,s,t):b\in B_{s,t}\}\subseteq \bR^2\times\bR^2 = \bR^4.\]

${\cC}$ can also be explicitly parametrized, via the rule:

\[(\theta,s,t)\mapsto (r,\theta, z,w) = (2 + s\cos(5\theta + 2\pi t),\theta,s,t)\subseteq\bR^2\times\bR^2,\]

where the first $\bR^2$ is given polar coordinates.  Of course, this attempted static animation is now a 3-dimensional hypersurface in $\bR^4$, making it impossible to 3D print.  This is exactly the type of situation studied in Section <ref> where we look at deformations of Julia Sets which have a 2-dimensional parameter space (the Mandelbrot set), so rather than just give up here, we will briefly describe a way to cut down a dimension.

Let $I$ be an interval and consider a path $\gamma:I\to\bR^2$, which we express coordinatewise as $\gamma(t) = (x(t),y(t))$.   Then for every $t\in I$ we can consider the frame $D_t:=C_{x(t),y(t)}$ defined as in Equation (<ref>).  We are now back in the setup of Definition <ref> and can form the static animation

\[\cD_\gamma:=\{(d,t):t\in I, d\in D_t = C_{x(t),y(t)}\}\subseteq\bR^3.\]

For example, if we let $\gamma:[.2,1]\to\bR^2$ be the path $\gamma(t) = (1-t,1-t)$ (consisting of travelling from $(1,1)$ to $(.2,.2)$ in a straight line along the diagonal), then $\cD_\gamma$ recovers the static animation $\cB$ from Section <ref>.

The key takeaway is the following: if one has a continuously varying family of 2D objects which vary over a higher dimensional parameter space $P$, one can still create 3D static animations by considering one parameter variations associated to paths within $P$.

§ DEFORMING THE PYTHAGOREAN TREE FRACTAL

The Pythagorean tree fractal is a recursive plane fractal, constructed as follows.  First, one attaches the the hypotenuse of a right triangle to the side of a square.

(0,0) – (0,3) – (3,3) – (3,0) – cycle;
at (4,1.5) $\Longrightarrow$;
(5,0) – (5,3) – (8,3) – (8,0) – cycle;
(5,3) – (6.5,4.5) – (8,3) – cycle;

Then two squares are attached to the legs of this triangle, and the process repeats, with the hypotenuse of a right triangle affixed to the opposite side of each new square.

(5,0) – (5,3) – (8,3) – (8,0) – cycle;
(5,3) – (6.5,4.5) – (8,3) – cycle;
(5,3) – (6.5,4.5) – (5,6) – (3.5,4.5) – cycle;
(6.5,4.5) – (8,6) – (9.5,4.5) – (8,3) – cycle;

at (4,1.5) $\Longrightarrow$;
[shift = (2,0)]
(5,0) – (5,3) – (8,3) – (8,0) – cycle;
(5,3) – (6.5,4.5) – (8,3) – cycle;
(5,3) – (6.5,4.5) – (5,6) – (3.5,4.5) – cycle;
(6.5,4.5) – (8,6) – (9.5,4.5) – (8,3) – cycle;
(3.5,4.5) – (3.5,6) – (5,6);
(8,6) – (9.5,6) – (9.5,4.5);

Squares are added to the legs of these triangles, and the process continues recursively forever.

[scale = 1]

[scale = 1]

Notice all of the triangles in the construction above were congruent, in particular they were isosceles right triangles with acute angles of $45^\circ$.  One could ask what happens if one used, say, a 30-60-90 triangle at each step instead.  As before, one begins by aligning the hypotenuse of a 30-60-90 triangle with the top of the square.

(0,0) – (0,3) – (3,3) – (3,0) – cycle;
at (4,1.5) $\Longrightarrow$;
(5,0) – (5,3) – (8,3) – (8,0) – cycle;
(5,3) – (5.75,4.29) – (8,3) – cycle;
at (6.9,3.25) $30$;

Then, following a similar process, two squares are attached to the legs of this triangle, although in this case they will not be the same size, instead matching the two differently sized legs.  Continuing, the hypotenuse of an appropriately sized 30-60-90 triangle is attached to the opposite side of each square, making sure to maintain consistent orientations.

(0,0) – (0,3) – (3,3) – (3,0) – cycle;
(0,3) – (.75,4.29) – (3,3) – cycle;
at (1.9,3.25) $30$;
(0,0) – (0,3) – (3,3) – (3,0) – cycle;

(0,0) – (0,3) – (3,3) – (3,0) – cycle;

at (6,1.5) $\Longrightarrow$;
(0,0) – (0,3) – (3,3) – (3,0) – cycle;
(0,3) – (.75,4.29) – (3,3) – cycle;
at (1.9,3.25) $30$;
(0,0) – (0,3) – (3,3) – (3,0) – cycle;
(0,3) – (.75,4.29) – (3,3) – cycle;
at (1.9,3.3) $30$;

(0,0) – (0,3) – (3,3) – (3,0) – cycle;
(0,3) – (.75,4.29) – (3,3) – cycle;
at (1.9,3.25) $30$;

As before, we continue this process recursively, this time observing that the global effect is somewhat different than when we used isosceles triangles.

[scale = 1]

[scale = 1]

A right triangle is determined (up to congruence), by its smallest angle $\theta$.  In particular, for each angle $0<\theta\le45$, one can follow this construction using a right triangle whose acute angles are $\theta$ and $90-\theta$.  By varying the angle $\theta$, one obtains a varying family of planar fractals parametrized by $\theta$.  Below we display the resulting fractal for some values of $\theta$ between $45$ and $30$.  The images seem to indicate that these fractals interpolate between the two that we have seen so far, suggesting the creation of a (static) animation.

at (.55,-.7) $\theta=45$;
at (.55,-.7) $\theta=40$;

[shift = (0,-6)]
at (.55,-.7) $\theta=35$;

[shift = (6,-6)]
at (.55,-.7) $\theta=30$;

Indeed, this puts us precisely in the setup of Definition <ref>.  For angles $\theta$ in the interval $[30,45]$, we can let $P_\theta\subseteq\bR^2$ be the Pythagorean tree fractal constructed using the right triangle of minimal angle $\theta$ in the recursive step, considered as a filled subset of the plane.  This gives us an entire deforming family of Pythagorean tree fractals, which one can consider as an animation on a screen [3].[See Matthew Conroy's animation as $\theta$ varies from 0 to 90: <https://www.madandmoonly.com/doctormatt/webimages/animations/pythagoreanFractal2.gif>.]

3D printing gives us a new way to observe the transformation of this fractal as $\theta$ varies.  Indeed, our approach is to follow the framework of Section <ref>, and to consider the deformation space as a solid static animation in $\bR^3$:

\[\cP = \{(x,\theta):\theta\in[30,45],x\in P_\theta\}\subseteq\bR^3.\]

We can use a computer to model $\cP$ using CAD software (we did this using the OpenSCAD [16]), and create a 3D object which can be 3D printed, held in your hand, and studied from all points of view.

In this way, 3D printing gives us a new avenue to illustrate perturbations and deformations of this family of fractals.  For one, we get a new dimension where we can observe the movement$-$for example, by considering where the slope becomes more exaggerated, it becomes immediately visible that the change in angle has more of an effect toward the outer edges of the fractal, and more of an effect on the 30 degree side than the 60 degree side.  Even more, one obtains a physical and tactile way to experience a mathematical concept, opening new doors for broad mathematical interaction.

§ DEFORMATIONS OF JULIA SETS

§.§.§ Joint with Bernat Espigulé

Let $\bC$ denote the field of complex numbers, fix a complex number $c\in\bC$, and consider the function $f_c:\bC\to\bC$ given by the rule $f_c(z) = z^2+c$.  What happens when we iterate this function?  More explicitly, if we fix some $z_0\in\bC$, what can we say about the following sequence of complex numbers?

\begin{eqnarray*}
z_1 &=& f_c(z_0)\\
z_2 &=& f_c(z_1) = f_c(f_c(z_0))\\
z_3 &=& f_c(z_2) = f_c(f_c(f_c(z_0)))\\
\vdots&&
\end{eqnarray*}

A fundamental question asks how the behavior of this sequence changes as we vary $z_0$, or as we vary $c$.  Let's get a basic idea of what can happen by letting $c=0$.

Consider $f_{0}(z) = z^2$.  Then

\[f_0^n(z_0) = \underbrace{f_0\circ f_0\circ\cdots\circ f_0}_{n-\text{times}}(z_0) =

If $|z_0|<1$ then this sequence converges to 0, if $|z_0|>1$ then this sequence escapes to infinity, and if $|z_0|=1$ then this sequence stays on unit circle $|z|=1$.  In particular, the set

\[\{z,f_0(z),f_0^2(z),f_0^3(z),\cdots\}\]

is bounded if and only if $|z|\le1$.

It appears that when iterating $f_c$ from a fixed starting point we can expect the outputs to either escape to infinity, or else to all be contained in a bounded region of $\bC$.  Filled Julia sets consist of the starting points which produce bounded sequences.

Let $c\in\bC$ be a complex number and consider the function $f_c(z) = z^2+c$.  The filled Julia set associated to $c$ is the set:

\[J_c:=\{z\in\bC : \text{ the set }\{z,f_c(z),f_c^2(z),...\}\text{ is bounded}\}\subseteq\bC.\]

In Example <ref> we determined that the filled Julia set $J_0$ is precisely the complex unit disk $\{z:|z|\le1\}\subseteq\bC$.  For other complex numbers, the situation becomes much richer.  For example, below are plots of $J_c$ for various values of $c$[The 2D images of filled Julia sets in this section were rendered in Mathematica [26]].

$c = -1.25$

$c = .365 + .1*i$$c = -0.12+0.74i$

A central question asks how $J_c$ changes as $c$ does.  For example, observe the variation of the filled in Julia sets for polynomials $f_c$ with values $c=-.75$, $c=-.74+.05i$, and $c=-.73+.1i$ respectively.

It appears that if $c$ changes a little bit, then so does $J_c$,[In fact, the situation is quite a bit more subtle than this.  We discuss this a little further in Section <ref>. The question is studied systematically by Douady [7].] suggesting that we can use these filled Julia sets as the frames of a static animation.  But, as discussed in Section <ref>, there is an obstacle.   This is because the varying parameter $c$ is a complex number, and so we obtain a 2-dimensional stack of frames.  When we form a deformation space exactly as in Definition <ref>, we obtain a subset of $\bC\times\bC$:

\[\cJ = \{(c,z):c\in\bC,z\in J_c\}\subseteq\bC^2.\]

As a hypersurface in 4D space, this object is impossible to 3D print.  In order to get something we can print, we will choose a 1-dimensional path through the parameter space (c.f. Section <ref>).  But which paths should we choose?  To answer this it is important to understand a little more about the geometry of the parameter space of filled Julia sets.

Each filled Julia set we have looked at so far has been connected and has a clearly defined interior, but for many values of $c$ the set $J_c$ is a totally disconnected Cantor set <cit.>.  For the purposes of 3D printing it makes sense to focus on values $c$ where $J_c$ is connected, which is a fascinating and celebrated collection of complex numbers in its own right.

The Mandelbrot set is the set:

\[\cM:=\{c\in\bC: J_c\text{ is connected}\}\subseteq\bC\]

A more familiar definition of the Mandlebrot set is perhaps

\[\{c\in\bC: \{0,f_c(0),f_c^2(0),...\}\text{ is bounded}\} = \{c\in\bC: 0\in J_c\}.\]

The fact that these are equivalent is <cit.>.  The Mandlebrot set famously has the following shape:

To cut down our deformation space $\cJ$ by one dimension and obtain something 3D printable, we mirror Section <ref> and pick a 1-dimensional subset of the Mandlebrot set defined by a path.  As the behavior of the Julia set $J_c$ is closely related to the location of $c$ in the Mandlebrot set, we will use the geometry of the Mandlebrot set to inform our choice of path.

§.§ Turning Paths into Prints

Let $I\subseteq\bR$ be an interval, and consider a continuous map

\[\gamma:I\to\bC,\]

whose image is contained in $\cM$.  For each $t\in I$, we can consider the filled Julia set $J_{\gamma(t)}$.  This gives us a family of filled Julia sets parametrized by $t$, putting us squarely in the situation of Definition <ref>.  Therefore the deformation space associated to this data is a static animation:

\[\cJ_\gamma:=\{(z,t):t\in I, z\in J_{\gamma(t)}\}\subseteq\bC\times\bR\approx\bR^3.\]

This sets up a general framework for creating static animations illustrating the deformation of filled Julia sets associated to a path in the Mandlebrot set.  It only remains to choose a path.  The boundary of the Mandlebrot set partitions the collection of filled Julia sets into two types: the filled Julia sets that are connected, and those that are not.  In particular, it behaves like a critical point, and we begin by studying how filled Julia sets deform near the boundary, or along it.  Our initial exploration was to traverse the upper boundary of the main cardioid.

One can deduce from <cit.> that the boundary of the main cardioid can be parametrized as follows.
\begin{eqnarray*}
x(t) &=& \frac{1}{2}\cos(t)(1-\cos(t)) + \frac{1}{4},\\
y(t) &=& \frac{1}{2}\sin(t)(1-\cos(t)).
\end{eqnarray*}

Then the path $\gamma:[0,\pi]\to\bC$ given by $\gamma(t) = x(t) + iy(t)$ parametrizes the upper boundary of the main cardioid.  Associated to this path there is a static aimation $J_\gamma\subseteq\bC\times\bR$, a 3D print of which is pictured below.

§.§ Parabolic Implosions

We suggested in Section <ref> that we were mostly interested in static animations with continuously deforming frames, but is this true for Julia sets?  This question was studied by Douady [7] who showed that in general, it is not.  The fact is, there are some special points at which $c\mapsto J_c$ does not vary continuously in $c$, even along our path $\gamma$ traversing the main cardioid.  To illustrate an example of this we look at the filled Julia sets associated to $\gamma(0)$ and $\gamma(\varepsilon)$ for some very small $\varepsilon$ (in this case $\varepsilon = .075$).  These are the frames $J_c$ associated with the complex numbers $c_0 = .25$ and $c_{\varepsilon} =0.251402 + 0.000105321 i$ respectively.

Although the overall shape is similar, there is a sudden appearance of a number of gaps.  This is an example of a phenomenon called parabolic implosion <cit.>, which is of considerable interest [20], but this curious phenomenon isn't clearly visible in our model, as the filled top layer of the print obstructs the view of the suddenly appearing gaps.  One way to see the parabolic implosions is to instead consider a restriction of our path to a slightly smaller interval.  There is a parabolic implosion also at $\gamma(\pi)$, so we consider the path $\gamma:[\varepsilon,\pi-\varepsilon]$ for very small $\varepsilon$.   The static animation associated to this path allows a better view of formation and deformation of the parabolic implosions, and produces a striking visual effect.

Having a 2-dimensional parameter space initially seemed like an obstacle in creating static animations, but the philosophy of Section <ref> turns this into a richness instead.  The freedom of having a static animation associated to each path allows us to focus on interesting areas and study them from multiple points of view.  We exploited this by choosing paths that remain in the vicinity of these parabolic points, allowing us to do local explorations around the parabolic implosions of the cusp ($c\approx\frac{1}{4}$).

We do a similar exploration around the parabolic implosion near the so-called Basilica ($c\approx-\frac{3}{4}$).

§.§ Non-Connected Julia Sets

So far we have restricted ourselves to considering filled Julia sets $J_c$ for values of $c$ in the Mandlebrot set.  What happens if we try to also consider $J_c$ for $c\in\bC\setminus\cM$?  These Julia sets will not be connected, and in fact, the Julia-Fatou theorem says that $J_c$ is either connected, or else totally disconnected and homeomorphic to a Cantor set <cit.>.  We experimented with how our methods (cf. Section <ref>) adapted to a path which exited the Mandlebrot set, and obtained a curious result.

Where do these towering mountains come from?  In some sense, they are an artifact of computer graphics.  What happens is that our method for generating images of Julia sets uses density to color pixels.  That is, if a pixel contains many points in the Julia set, it will be colored black, and otherwise it will be left white.  As a result, for complex numbers $c$ which are very close to the boundary of $\cM$, the screen displays $J_c$ as a union of several connected components, since these are where there are high densities of points in the Julia set.  For example, here is our image for $c=.28+.01i$, a complex number outside of the Mandlebrot set, but just barely.

We emphasize that this is an artifact of computer graphics, the set $J_c$ is totally disconnected.  What appears to be a connected region is merely a part of the plane where there is a high concentration of points which are bounded under the dynamical system of iterating $f_c$.  As $c$ progresses further away from $\cM$, these high concentration areas spread out and we can observe a rapid progression toward a more sparse Julia set.  It is worth remarking that this print, although not mathematically rigorous in the purest sense, provides mathematical insights about concentrations of points in Cantor-like Julia sets and how they deform and spread out

The fact that an artifact of computer graphics can both create interesting aesthetics and contribute meaningfully to our understanding is interesting in its own right, although beyond the scope of this article.  Discussions about techniques, artifacts, and limitations of generating Julia sets using computer graphics are discussed by Milnor <cit.>, and in great depth by Chéritat [2].

§.§ Other Directions

The static animations of filled Julia sets exhibited so far have been limited to paths on and near the boundary of the main cardioid,  but there is an entire menagerie of static animations that can be constructed simply by choosing paths through the Mandlebrot set.  We invite the reader to explore these possibilities, and to this end we have included a Mathematica notebook called  in the public github repository accompanying this paper [5].   This notebook allows the user to specify a parametrized path through $\bC$, and will generate the frames of the associated static animation.  The reader can then follow the procedure described in Section <ref> to construct the resulting static animation. For now, we suggest a few of the many paths worth exploring, and highlight a related project led by Caroline Davis [1].

§.§.§ Traversing Other Boundaries

To get a more complete picture of how filled Julia sets deform along the boundary of the Mandlebrot set we should also consider static animations associated to the boundaries of other bulbs.  For example, the bulb directly left of the main cardioid (called the period 2 hyperbolic component) is a perfect circle of radius 1/4 centered at $-1$.  Below is the static animation associated to traversing the upper semicircle of this component.

The other hyperbolic components do not have boundaries that are so easy to parametrize.  Giarruso and Fisher explicitly parametrized the period 3 hyperbolic component [11], and the parametrization is much more complicated than the two we have seen so far.  Beyond that, little has been worked out explicitly, leaving the question of how to correctly traverse the boundary as an interesting and difficult one.  In fact, the question of whether the boundary of the Mandelbrot set can be parametrized as a curve at all is related to the famous Mandlebrot Locally Connected (MLC) conjecture of Douady and Hubbard [9],[8].  It probably isn't quite as difficult to find the boundary of a specific hyperbolic component, nevertheless, the problem of finding parametrizations of subsets of the boundary of the Mandlebrot set is connected to deep research questions in complex dynamics.

§.§.§ Traversing Radii

The boundary is not the only interesting place to look.  Another direction (suggested by Caroline Davis), would be to consider the radii of the hyperbolic components.  Indeed, each bulb of the Mandlebrot Set gives rise to Julia sets with different discrete invariants (such as the periods of the attracting cycles <cit.>).  It is also known that each component is canonically holomorphic to a circle [18], so that each bulb has a canonically defined center and radii.  Tracing the radii between the center of one bulb to the center of another gives rise to a canonical deformation between two Julia sets with distinct discrete invariants, something that we are interested in 3D printing.  An important step is to obtain explicit parametrizations of these radii.  As with the boundaries, these are easy to find for the main cardioid, and the hyperbolic component of period 2 (which is already a disk), and can be extracted from [11] in the period 3 case.  Below is an image made by Walter Hannah [12] depicting these radii in the main cardioid.

[12]

§.§.§ Related Objects in Holomorphic Dynamics

We have only considered Julia sets associated to quadratic polynomials, but one can also iterate other functions$-$such as polynomials of higher degree or rational functions$-$and study their associated Julia sets.  An undergraduate research group at Indiana University, led by Caroline Davis, studied a family of polynomials whose associated Julia sets interpolate between the Apolonian and Sierpinski gaskets [1].  Below is a print of the associated static animation, printed by Bethany Mussman.

§.§ Static vs. Computer Animations of Julia sets

We are certainly not the first to visualize how Julia sets deform according to carefully chosen paths in the Mandlebrot set.  There are many computer animated visualizations of deforming Julia sets, including witnessing degenerations from connected Julia sets to Cantor sets, and exploratations near parabolic implosions.[In fact, the suggestion to 3D print parabolic implosions was made to the author by John H. Hubbard as he was demonstrating an animation he made of this phenomenon.]  The shift from computer to static animations brings with it advantages and disadvantages.  For example, computer animations give a more direct view of the internal structure, which can be hidden by the walls of the static animation, although this can be mostly remedied by adjusting the start and endpoints of the chosen path like we did in Section <ref>.   Nevertheless, there is a completely novel perspective the static animation brings to the table, something that is missing from the computer animation: the side view.  From the horizontal perspective we can observe how the different bulbs on the edge of these Julia sets move at different rates as the Julia set deforms, and these bulbs braid around each other in a way that seems to interpolate between the mathematical and the organic.  The following image gives a particularly good illustration of this phenomenon.

§ CREATING A STATIC ANIMATION

When entering the world of 3D printing, it is easy to get overwhelmed by the large number of CAD programs and 3D modeling tools, many of which can be difficult to learn, especially as a beginner.  On the other hand, since static animations consist of stacks of 2-dimensional objects, one can hope to build a 3D static animation using purely 2-dimensional techniques.  In this section we describe a workflow that does exactly that,[This particular workflow was first shared with the author by Bernat Espigulé, and is how we created the static animations of filled Julia sets described in Section <ref>.] so that anyone who can make a 2-dimensional computer animation can begin designing 3D prints.  After describing the general process, we will then run through an explicit tutorial, building the static animation of deforming polar flowers from Section <ref>.  We intend this section to serve as an invitation for beginners to mathematical 3D printing, while also providing a new modeling technique to the experts.

The main idea is simple, and can be thought of as a sort of reverse tomography.  While tomography is an imaging technique which extracts a collection of 2D cross sections from a 3D object, we begin with a collection of 2-dimensional frames, which we stitch together into a 3D model whose horizontal cross sections will be precisely the frames we begin with.

A summary of the workflow is the following: create a collection of 2-dimensional frames however you like, and upload them to a piece of software that can bond them together into a 3D model.  We use the open source software, UCSF Chimera [21],[25].  It has many uses, but we are most interested in its use as medical imaging software.  Roughly, MRI tomography takes a series of images of an organ at various depths inside the body, and these images can be fed to Chimera to reverse engineer a 3D model of the organ.  We can exploit this functionality to make our static animations.

§.§ Overview

Fix a collection of frames $B_t$, parametrized by values $t$ in some interval $I$.  We will describe a general process that can be followed to model the static animation

\[\cB = \{(b,t):b\in B_t\}\subseteq\bR^3,\]

as a 3D object file.  These steps assume the reader can already produce a 2-dimensional image of each frame $B_t$.

§.§.§ Step 1: Make a collection of frames

Subdivide your interval into a reasonable number of frames.  Keep in mind that a common layer thickness for a 3D print is about .2mm, and the larger at home 3D printers allow a maximum height of no more than 300mm, so the resolution of your 3D printed object will probably not benefit from having more than 1500 frames (although industrial 3D printers can do better).  This gives a finite number of $t_i\in I$ which you will sample from.

For each $t_i$, produce a PNG image of $B_{t_i}$, coloring points $p\in B_{t_i}$ white, and points $q\notin B_{t_i}$ black.  This is the color scheme and file type that Chimera recognizes as cross sections of a 3D object.  Save the entire collection of images in a single folder and choose filenames so that the alphabetical ordering of the images corresponds to the correct ordering of the frames, this way they will be imported in the correct order during the next step.

You can create the frames of the animation any way you would like.  For the prints in Section <ref>, we fed a sampling of points on our parametrized path to Mathematica's JuliaSetPlot function [26], which we modified to color points in white and black as needed.

Warning: You may need to remove the alpha channel.  The PNG graphics file format can store color images with an RGB color palette (storing numerical values for the red, green, and blue of each pixel), or an RGBA color palette, which has a fourth alpha channel for transparency.  Chimera will reject PNG graphics with an alpha channel, so  you may need to manually remove the alpha channel from your frames.  To do this we wrote a short script called , included in the repository accompanying this paper [5].

§.§.§ Step 2: Convert the frames to a solid in Chimera

From the Chimera home screen, choose .  Navigate to and highlight the entire set of images from Step 1, and select .  A rough version of your model should become visible, as well as a volume viewer menu with various settings.  There are 3 important settings to be aware of:

* Step: A dropdown menu labelled step, with choices of values 1,2,4,8, and 16 (4 is the default).  This roughly controls the vertical resolution.  Choose 1 (the highest vertical resolution).  Your model should immediately look much sharper.

* Level: A slider to select the level, which is a number between 0 and 255.  This roughly controls the tolerance for edge detection.  We often have to play around with this a bit.

* Style: Radio buttons to select either surface, mesh, or solid.  We always choose surface.

Once you have found settings to your liking, you can export the model.  Select   You can choose from a number of file types.  The standards for 3D printing are OBJ or STL.  Depending on the size of the image and number of frames, this can take a long time.  For the particularly detailed animations of Section <ref> we needed to let it run for several hours.[For reference, the computation was done on a modestly powered Dell XPS 15 7590 with an Intel Core i7-9750H (12MB Cache, up to 4.5 GHz, 6 cores).]

§.§.§ Step 3: Slice and Print

Your model should now be ready to open up in a slicer (for example Cura [4]), and then sliced and sent to a 3D printer.  The full process depends on your 3D printer, and is beyond the scope of this article, but we'd like to include a few tips to help avoid failure at this step.  We point to a few resources for slicing and printing 3D models at the end of this subsection.

Avoid steep slopes and overhangs.  An extrusion based 3D printer will print your static animation by printing each frame on top of the last, using the previously printed frames as a support base for the next ones.  If the frames $B_t$ vary too wildly with $t$, the printer may end up trying to extrude part of a frame with no supporting structure underneath.  These unsupported layers of plastic building up will lead to unsatisfactory results and failed prints.  As a rule of thumb, most 3D printers can cleanly print overhanging structures with an angle up to somewhere between 45 and 60 degrees, so if you would like to print a wildly deforming family, you may need to do some experimentation to get something that works.[Alternatively, you can employ the use of support material, extra printed material which holds up excessive overhangs and can be removed after.]

Lay your model flat.  Chimera won't automatically export a model which lies flat, instead exporting the model aligned to the perspective in the viewing window.  Sending the tilted model to an extrusion based 3D printer will likely cause issues with adhesion and support, and could result in a failed print.  This can be avoided if you don't adjust the view of the model in Chimera's viewscreen between importing your frames and exporting the STL.  Otherwise you must rotate your exported model so that it lies flat on the buildplate.  This can be done in Cura using the Lay Flat command.

More resources for the printing process.  Although beyond the scope of this paper, there are many places to turn for a thorough introduction to the use of slicers and 3D printers to turn a digital model like your static animation into a physical model.  For example, the author's lecture during the 2022 Joint Mathematics Meetings mini-course 3D Printing: Challenges and Applications was a tutorial on slicing and printing, and the slides are publicly available [6].  There are also many detailed written tutorials (for example [10]) and youtube tutorials (for example [19]).

§.§ Tutorial

We now demonstrate this process step-by-step with an example, constructing the static animation of deforming polar flowers introduced in Section <ref>.  The reader is invited to follow along on their own machine.  To this end we have provided code snippets within this document, as well as the complete collection of code and files in a public Github repository [5].  For more detail on where to find what you need to follow along, see Section <ref> below.

In Section <ref> we computed a complete parametrization of this particular deformation space in $\bR^3$.  With this parametrization there are many more direct ways to produce a 3D model (for example, plotting in Mathematica or Sagemath and then exporting as an STL).  Nevertheless, one often encounters deformation spaces which do not admit such a parametrization (for example, those in Sections <ref> and <ref>).  One should think of this tutorial as an introduction to this particular process, rather than the most efficient way to make this particular model.

The frames of the static animation in this tutorial are polar curves whose equations we computed in Section <ref> to be $r = 2 + t\cos(5\theta + 2\pi t)$ for values $t$ in the interval $I = [0.2,1]$.  Since we are making a 3-dimensional solid, we will replace these curves with the regions they bound, so that our resulting object will be filled rather than hollow.  The frames we will use will therefore be the following sets

\[B_t = \{(r,\theta): r\le 2+t\cos(5\theta + 2\pi t)\}\subseteq\bR^2.\]

For example, $B_1$ is now the entire shaded region below, rather than just the boundary curve.

\[\includegraphics[width=.25\textwidth]{images/regionStar.pdf}\]

We now run through all the steps to model the static animation $\cB = \{(b,t):t\in I,b\in B_t\}$ as a 3D object file.

§.§.§ (Tutorial) Step 1

We subdivide our interval to create 161 frames, so that $t$ progresses between $0.2$ and $1$ at steps of $.005$ units.  Our list of $t$ values is therefore

\[T = \{.2, .205, .21, .215,\cdots,.990,.995,1\}.\]

We need to make a PNG image for each frame, drawing each region $B_t$ in white, against a black background.  Let's walk through a way to do this in python.  We need the  and  libraries.

import matplotlib.pyplot as plt
import numpy as np

Chimera expects white images against a black background, so we first set the color of the background we are drawing against.

For each frame, we'd like to draw the polar curve $r = 2t\cos(\theta+2\pi t)$.  To do this,  will sample a number of points and interpolate.  We use numpy's  command to create the list of $\theta$ values to sample from.

NUM_PTS = 1000
theta = np.linspace(0,2*np.pi,NUM_POINTS)

We then calculate the corresponding $r$ values to finish computing the list of points.  Notice we set $t=1$ in this snippet, but this will change from frame to frame.

t = 1 #This varies on each frame
r = 2 + t*np.cos(5*theta + 2*np.pi*t)

Now whats left is to let  draw and fill the curve.

plt.polar(theta,r,alpha=0) #Draws the curve
plt.fill_between(theta,0,r,color='white') #Fills the curve
plt.axis('off') #We just want the flower, no axes

Here is the output with $t$ values 1 and 1/2 respectively.

To create the entire stack of frames we loop through each value of $t\in T$ and export a PNG built using the the code above, choosing filenames in alphabetical order so that they are arranged correctly when imported by Chimera.  We can do this by looping through the last two snippets of the code above and saving the output for each frame.

NUM_FRAMES = 161 #Set the number of frames

for i in range(NUM_FRAMES):
t = 0.2 + .005*i
r = 2 + t*np.cos(5*theta + 2*np.pi*t) #The boundary curve
plt.polar(theta,r,alpha=0) #Draws the curve
plt.fill_between(theta,0,r,color='white') #Fills the curve
plt.axis('off') #We just want the flower, no axes

#Now export and clear the plot.
imageName = 'flowerStack/flowers_'+str(1000+i)+'.png'
plt.savefig(imageName) #Export the image
plt.cla() #clear the plot before looping back

The images are now saved as  through .  The reason we started counting at 1000 is because Chimera reads alphabetically in such a way that 10 comes before 9 because it starts with a 1.  In particular, if we started counting at 0 Chimera would arrange the imported frames in an incorrect order during Step 2.  We now have our frames images all saved as PNG images.

In our accompanying github repository [5] we include a python program  for creating the entire stack of frames, as well as an implementation in Mathematica, .

§.§.§ (Tutorial) Step 2

Open Chimera, select browse, and highlight all of the PNG images from step 1.

Click open.  Chimera stitches all of the polar flowers into a single solid.

It looks a bit rough.  By changing step to 1 and adjusting the level we arrive at a much finer model.  We also choose surface in the style setting (if it isn't already selected).

This looks much better.  Here are a few other points of view of the model Chimera produced.

We now export the model from Step 2 as an STL.  To do this we merely export the scene, choosing the STL file extension.  We should make sure to export the model before rotating it at all in the viewscreen.  This way it will be oriented correctly when opened in a slicer, lying flat against the buildplate (see the Step 3 outline above).

For more complicated animations this can take a long time$-$for example, if you need more frames, or a higher resolution to capture the complex behavior you are interested in.  This animation, though, is relatively simple, and on the author's current machine[HP ProBook 445 G8, AMD Ryzen 5 5600U with Radeon Graphics 2.30 GHz, RAM: 16 GB] it took about 3 minutes to export.

Warning. When importing your frames to Chimera at the beginning of this step, you may have encountered the following error message:

This is because the frames were saved as PNG files with an alpha layer.  To fix this, you must remove the alpha layer.  One can do this in Python using the Python Imaging Library (PIL).

from PIL import Image
image = Image.open(fileName).convert('RGB') #conversion step

In the snippet above,  should be a path to your image.  In the github repository accompanying this manuscript [5] we include a script, , which can loop through your entire image stack and remove the alpha layer from each frame.  The program  accompanying this tutorial already removes the alpha layers from your stack, so this step will only be necessary if you created your image stack elsewhere.

§.§.§ (Tutorial) Step 3

The model is now ready to be opened in slicing software (we use Cura [4]), sliced, and then sent to a printer.

§.§ Where do I get...

Our intent with this tutorial is to have the reader follow along on their own machine.  This way they learn the workflow and can integrate it into their own mathematical practice.  In order to follow along, the reader will need the following tools.

* Chimera: At the time of writing of this document, UCSF Chimera [25] is available for free download here: <https://www.cgl.ucsf.edu/chimera/>.

* Cura: At the time of writing of this document, Untimaker Cura [4] is available for free download here: <https://ultimaker.com/software/ultimaker-cura>.

* Tutorial Files: All the files generated while making this tutorial are available at our public github repository [5].  This includes:

* A python program generating the image stack.

* A Mathematica notebook generating the image stack.

* A python program for removing the alpha layer from a PNG stack.

* The PNG images that serve as the frames of the animation.

* The STL exported from Chimera.

If the reader would only like to follow along for one step, for example, stitching together the image stack in Chimera, they are welcome to download the frames from [5] and experiment in Chimera themselves.

§ FUTURE DIRECTIONS

This article only begins to scratch the surface of what the author believes to be a vast array of possibilities for 3D printed static animations in mathematical illustration and research. Indeed, with how ubiquitous the notion of deformations are in mathematics, and how useful animations on a screen have proved to be, we look forward to exciting new directions for this technique of illustration and experimentation.  We hope that with the techniques developed in Section <ref>, the reader will have the tools necessary both to incorporate 3D printing into their practice, and to find new applications for static animations.  We'd like to conclude with a couple of potential next steps.

§.§ Exploring the Mandlebrot set

As we saw in Section <ref>, the field of holomorphic dynamics offers a rich menagerie of static animations to investigate, including many stemming from paths within the Mandlebrot set. The process for creating static animations from paths in the Mandlebrot set is already established (and can be accessed in our public github repository [5]), so all one needs are parametrized paths and to begin to explore.

§.§ Chaos and the Double Pendulum

The double pendulum [23] is an example of a chaotic system.  In this system, attached to the mass at the end of a pendulum there is a second pendulum with another mass.

[15]

This is a dynamical system governed by a set of ordinary differential equations which is chaotic, that is, it is highly sensitive to its initial conditions.  Animations and experiments with double pendulums show that, even with two very similar starting points, the behavior of the moving pendulums will quickly diverge quite radically.  We imagine turning the computer animation into a static animation, building the subsequent positions of the double pendulum on top of the last.  We have begun such an exploration with the help of Daniel Piker using his physics engine, Kangaroo [22].

§.§ And more...

We invite the reader to include static animations in their own investigations.   Any place that a 2D animation has been useful, explore adapting it to a 3D static animation.  We hope that, with the gentle learning curve and connections to many mathematical disciplines, folks will find this to be an accessible avenue to the world of 3D printing in mathematics.  Indeed, as we showcased in Section <ref>, we have begun to see researchers lead undergraduate research investigations which use these techniques [1].

There are many fascinating applications of 3D printing to exploring and illustrating mathematics, and we hope that static animations will both become a part of that canon, and also serve as an invitation to the broader world of 3D printing in math.

[1]
Baltz, J., Davis, C., Mussman, B., and Ronnenberg, M.
3d printing julia sets.

[2]
Chéritat, A.
Techniques for computer generated pictures in complex dynamics.
<https://www.math.univ-toulouse.fr/ cheritat/wiki-draw/index.php/Mandelbrot_set>.

[3]
Conroy, M.
Animation of the Pythagorean tree fractal.

[4]
David Braam, Ultimaker.
Ultimaker Cura.

[5]
Dorfsman-Hopkins, G.
Deformation spaces and static animations: Accompanying code and
A github repository.

[6]
Dorfsman-Hopkins, G.
Tutorial and demonstration (slides).
Part of the 2022 JMM Mini Course, 3D Printing: Challenges and

[7]
Douady, A.
Does a Julia set depend continuously on the polynomial?
In Complex dynamical systems (Cincinnati, OH, 1994),
vol. 49 of Proc. Sympos. Appl. Math. Amer. Math. Soc., Providence, RI,
1994, pp. 91–138.

[8]
Douady, A., and Hubbard, J. H.
Étude dynamique des polynômes complexes. Partie I,
vol. 84 of Publications Mathématiques d'Orsay.
Université de Paris-Sud, Département de Mathématiques,
Orsay, 1984.

[9]
Douady, A., and Hubbard, J. H.
Étude dynamique des polynômes complexes. Partie II,
vol. 85 of Publications Mathématiques d'Orsay.
Université de Paris-Sud, Département de Mathématiques,
Orsay, 1985.
With the collaboration of P. Lavaurs, Tan Lei and P. Sentenac.

[10]
Dwamena, M.
How to use cura for beginners – step by step guide & more.
[Online; accessed November 12, 2022].

[11]
Giarrusso, D., and Fisher, Y.
A parameterization of the period $3$ hyperbolic components of the
Mandelbrot set.
Proc. Amer. Math. Soc. 123, 12 (1995), 3731–3737.

[12]
Hannah, W.
Internal rays of the mandlebrot set, 2006.
Cornell University, Undergraduate Honors Thesis.

[13]
Hartshorne, R.
Algebraic geometry.
Springer-Verlag, New York, 1977.
Graduate Texts in Mathematics, No. 52.

[14]
Hatcher, A.
Algebraic topology.
Cambridge Univ. Press, Cambridge, 2000.

[15]
Double pendulum svg.
[Online; accessed June 28, 2022].

[16]
Marius Kintel, Claire Wolf.

[17]
Milnor, J.
Dynamics in One Complex Variable. (AM-160): Third Edition.
Princeton University Press, 2006.

[18]
Milnor, J.
Hyperbolic components.
In Conformal dynamics and hyperbolic geometry, vol. 573 of Contemp. Math. Amer. Math. Soc., Providence, RI, 2012, pp. 183–232.
With an appendix by A. Poirier.

[19]
Now, D.
Cura 3d slicer for beginners.
[Online; accessed November 12, 2022].

[20]
Peters, H., and Vivas, L.
Parabolic implosion.
Notices Amer. Math. Soc. 67, 8 (2020), 1095–1103.

[21]
Pettersen, E., Goddard, T., Huang, C., Couch, G., Greenblatt, D., Meng,
E., and TE., F.
UCSF Chimera–a visualization system for exploratory research and
J Comput Chem 25, 13 (2004), 1605–12.

[22]
Piker, D.

[23]
Richter, P. H., and Scholz, H.-J.
Chaos in classical mechanics: the double pendulum.
In Stochastic phenomena and chaotic behaviour in complex systems
(Flattnitz, 1983), vol. 21 of Springer Ser. Synergetics. Springer,
Berlin, 1984, pp. 86–97.

[24]
Strogatz, S. H.
Nonlinear Dynamics and Chaos.
Addison Wesley, Reading, MA, 1994.

[25]
UCSF Resource for Biocomputing, Visualization, and Informatics.
Ucsf chimera.

[26]
Wolfram Research Incorporated.
Mathematica, Version 13.0.0.
Champaign, IL, 2021.

§ APPENDIX: DEFORMATIONS IN ALGEBRAIC GEOMETRY

The author's initial motivation to create 3D printed static animations of continuously deforming mathematical objects came from encountering these types of constructions in algebraic geometry, and more specifically within deformation theory.  Here, one often will think of parameters of variation as a spacial variables and produce a deformation space of higher dimension to which witnesses the interpolation.  Indeed, if $C$ and $C'$ are varieties which can deform to one another, exhibiting them both as living within some larger deformation space $D$ gives a geometric framework wherein information can be passed between them via the geometry of $D$.  As this technique was the authors inspiration to begin creating static animations, we'd like to use this appendix to describe this connection in some more detail, using examples of flat families of algebraic curves.

Suppose you are studying a family of curves $C_t$ parametrized by $t\in k$ where $k$ is a field (little is lost imagining $k$ to be $\bC$ or $\bR$).  Suppose further that each curve is the zero set in $k^2$ of a polynomial equation in 2 variables $x$ and $y$.  The key observation is that this $t$ can be introduced as a new spacial variable, defining an equation, now in $x,y,$ and $t$, whose zero set in $k^3$ we will call $D$.  Then $D$ comes equipped with a projection $\pi:D\to k$, which plucks out the $t$-coordinate, and each member of the original family of curves can be recovered as the fibers of this projection:

\[C_t = \pi^{-1}(t)\subseteq D.\]

The remarkable thing about this viewpoint is that relationships between the various curves $C_t$ in the family can be understood in terms of the global geometry of the deformation space $D$, which can be studied as an algebraic variety in its own right. For example, one can use the geometry of $D$ to show that the genus of the curves in such a family is constant <cit.>.[The genus of a curve can be defined over any field,  but if $k=\bC$, then the curve lives as a subspace of 1 complex dimenion in $\bC^2$ and can therefore be thought of as a Riemann surface, allowing one to determine genus in the usual way] To connect with our 3D printing philosophy, the 3D printed static animation plays the role of the total deformation space $D$, which can be viewed as a vertical stack of the $C_t$ at height $t$.  Much like how the global geometry of $D$ provides new insights into the connections between the parametrized curves $C_t$, this paper begins to ask how the global structure of a 3D printed static animation reflects the connections between the parametrized objects it contains.

We give one explicit example of the deformation space of a family of curves.  To do this let $k=\bR$ and begin by considering the hyperbola $xy = 1$.

[->] (-3, 0) – (3, 0) node[right] $x$;
[->] (0, -3) – (0, 3) node[above] $y$;
[domain=.33:3, smooth, variable=, blue] plot (, 1/);
[domain=-3:-.33, smooth, variable=, blue] plot (, 1/);

Now let us perturb the equation, considering $xy=\frac{1}{2}$, or $xy=\frac{1}{4}$.

[->] (-3, 0) – (3, 0) node[right] $x$;
[->] (0, -3) – (0, 3) node[above] $y$;
[domain=.17:3, smooth, variable=, blue] plot (, .5/);
[domain=-3:-.17, smooth, variable=, blue] plot (, .5/);
[domain=.08:3, smooth, dashed, variable=, red] plot (, .25/);
[domain=-3:-.08, smooth, dashed, variable=, red] plot (, .25/);

We can continue decreasing this constant term, considering the curve $C_t$ given by the equation $xy=t$ for $t=\frac{1}{8},\frac{1}{16},...$.

[samples = 50]
[->] (-3, 0) – (3, 0) node[right] $x$;
[->] (0, -3) – (0, 3) node[above] $y$;
[domain=.04:3, smooth, variable=, blue] plot (, .125/);
[domain=-3:-.04, smooth, variable=, blue] plot (, .125/);
[domain=.02:3, smooth, dashed, variable=, red] plot (, .0625/);
[domain=-3:-.02, smooth, dashed, variable=, red] plot (, .0625/);

As $t$ approaches zero the curves appear to be converging to the axes.  And indeed, $xy=0$ if and only if $x=0$ or $y=0$, so that $C_0$ is precisely the union of the $x$ and $y$ axes.  We have therefore constructed a family of curves $C_t$, parametrized by $t\in\bR$, which continuously interpolate between hyperbolas and a union of 2 lines.  Of course, nothing is stopping us from treating $t$ as a spacial variable.  That is, we can consider the the surface in $\bR^3$ defined by the equation $xy=z$.

The horizontal cross section at height $z = t$ is precisely the curve $C_t$ given by the equation $xy=t$, and if one considers the associated contour lines, one recovers the picture of the deformation of hyperbolas to the planar axes.

In this way, the surface $xy=z$ is a higher dimensional space which witnesses the continuous deformation of hyperbolas to a union of 2 lines.


# Adaptive Local Implicit Image Function for Arbitrary-scale Super-resolution

###### Abstract

Image representation is critical for many visual tasks. Instead of
representing images discretely with 2D arrays of pixels, a recent study,
namely local implicit image function (LIIF), denotes images as a continuous
function where pixel values are expansion by using the corresponding
coordinates as inputs. Due to its continuous nature, LIIF can be adopted for
arbitrary-scale image super-resolution tasks, resulting in a single effective
and efficient model for various up-scaling factors. However, LIIF often
suffers from structural distortions and ringing artifacts around edges, mostly
because all pixels share the same model, thus ignoring the local properties of
the image. In this paper, we propose a novel adaptive local image function
(A-LIIF) to alleviate this problem. Specifically, our A-LIIF consists of two
main components: an encoder and a expansion network. The former captures
cross-scale image features, while the latter models the continuous up-scaling
function by a weighted combination of multiple local implicit image functions.
Accordingly, our A-LIIF can reconstruct the high-frequency textures and
structures more accurately. Experiments on multiple benchmark datasets verify
the effectiveness of our method. Our codes are available at
https://github.com/LeeHW-THU/A-LIIF.

Index Terms—  image super-resolution, implicit image function, image
restoration, deep learning

##  1 Introduction

Deep neural networks (DNNs) have been widely and successfully adopted in many
applications, such as image super-resolution. Currently, most existing main-
stream DNN-based super-resolution methods treat images as discrete 2D arrays
and learn an encoder-decoder style model for each up-scaling ratio. As such,
developers need to retrain or adopt different networks for different up-
scaling ratios.

Fig. 1: An example of amplified images generated by LIIF [1] and our A-LIIF at
the resolution scale $\times 12$. As denoted in the red box areas, LIIF [1]
has ringing artifacts around edges, while our A-LIIF can relieve this problem.

Recently, the local implicit image function (LIIF) [1] explores this task from
another perspective, opening the new era of image super-resolution.
Specifically, LIIF [1] shows that high-resolution (HR) images can be learned
as continuous functions of low-resolution (LR) features and pixel coordinates,
where the function can be easily modeled by a simple multi-layer perception
(MLP). Due to its continuous manner, users can exploit LIIF [1] for arbitrary-
scale super-resolution, without creating and training different models.

However, we observe that LIIF [1] often generates ringing artifacts around
edges, especially when the resolution scale is relatively large (as shown in
Figure 1 (b)). This problem is mainly because all predicted pixels share the
same MLP whereas pixels from different regions (such as smooth and edge areas)
may have different characteristics. It motivates us to predict pixels by
considering the pixel characteristics.

In this paper, we propose an adaptive local image function (A-LIIF) by
automatically modeling pixel differences with multiple MLPs. As shown in
Figure 2, our A-LIIF mainly consists of encoder, expansion network, and
multiple MLPs. Specifically, we design an expansion network to generate $K$
different adaptive scores of each element in the feature map, where $K$ is the
expansion rate. After that, the predicted pixel is produced by a weighted
linear combination of $K$ MLPs with the adaptive scores. Unlike the LIIF [1]
that suffers from ringing artifacts, our A-LIIF can relieve such problem and
reconstruct the details more accurately (seen in Figure 1(c)).

Fig. 2: The main pipeline of our A-LIIF. The red dots are the positions of
targeted pixels that we need to render.

In conclusion, the main contributions of this paper can be summarized as
follows:

  * •

We reveal that LIIF [1] usually suffers from structural distortions and
ringing artifacts around edges, mostly because all pixels share the same MLP.

  * •

We propose an adaptive local image function (A-LIIF) by modeling pixel
differences with multiple MLPs automatically to alleviate artifacts.

  * •

We conduct extensive experiments on five benchmark datasets to verify the
effectiveness of our A-LIIF.

##  2 RELATED WORKs

This section provides a brief overview of recent work on implicit neural
representations, expansion network, basis, and various super-resolution
methods related to our work.

###  2.1 Single Image Super-Resolution

Super-resolution is one of the low-level vision tasks that has been explored
in our vision field for decades, especially single image super-resolution
(SISR). With the rapid advancement of deep learning techniques, super-
resolution models implemented based on deep learning are also in the process
of breakthrough. Different types of deep learning strategies are widely
applied to target super-resolution related problems. The earliest were adopted
CNN algorithms SRCNN [2], and later, networks with larger residual structures
[VDSR [3], IRCNN [4], SRResNet [5], EDSR [6], RDN [7]] , while the concept of
non-localization is also very prevalent [SENet [8], Non-local [9], SAN [10]]
taking the SR quality to a new level, all these works have actively explored
the field of SR.

###  2.2 Arbitrary-Scale Super-Resolution

Arbitrary-scale super-resolution methods hold great promise for research. The
emergence of traditional single-image super-resolution methods, whose models
can only be used at a single scale, a phenomenon that makes it difficult to
land super-resolution modeling methods commercially in real life, greatly
surpasses previous SISR work in terms of convenience. The earliest arbitrary-
scale super-resolution methods can be traced back to MetaSR [11], which was
the first CNN-based arbitrary-scale SR method. The recent SOTA work LIIF [1]
proposes a new framework for arbitrary-scale super-resolution using implicit
neural representations, where feature maps, coordinates, and scaling factors
are input to the MLP after Encoder for computation to obtain RGB values.

Table 1: Quantitative comparison on DIV2K [12] validation set (PSNR (dB)). $\phi$ indicates ours implementation. EDSR-baseline [6] trains different models for different scales. Method | In-distribution | Out-of-distribution
---|---|---
×2 | ×3 | ×4 | ×6 | ×12 | ×18 | ×24 | ×30
bicubic [6] | 31.01 | 28.22 | 26.66 | 24.82 | 22.27 | 21.00 | 20.19 | 19.59
EDSR-baseline [6] | 34.56 | 30.82 | 28.84 | - | - | - | - | -
EDSR-baseline-LIIF [1]ϕ | 33.57 | 29.93 | 28.14 | 26.02 | 23.17 | 21.74 | 20.83 | 20.19
A-LIIF | 34.65 | 30.95 | 28.99 | 26.75 | 23.71 | 22.18 | 21.18 | 20.48

##  3 The Proposed METHOD

Our work is inspired by NeRF [13], which demonstrates that 3D scenes with fine
textures can be rendered with a neural rendering approach. A recent work named
LIIF [1] further demonstrates that implicit image functions can recover images
at arbitrary scales of super-resolution with accurate details. However, the
LIIF [1] approach usually suffers from ringing artifacts by using a shared MLP
for every pixel point. Thus, a natural question arises: does it use the same
implicit neural function even for very different edge features? For this
purpose, we design a novel adaptive local image function (A-LIIF) to relieve
such problem, which mainly consists of expansion network and multiple MLPs
modules. Besides, our method demonstrates that for each pixel point a
different MLPs should be used to recover high frequency textures. The main
diagram of our method is shown in Figure 2, which is detailed in the following
sections.

###  3.1 Expansion Network

The edge features of all pixel points are different when recovering an image,
however LIIF [1] ignores this importance and directly for all pixel points is
recovered using a shared implicit image function. We empirically note that
these cases of structural distortion and apparent artifacts arising from
neural representations can be resolved if an adaptive MLP layer is used. So
our A-LIIF feeds the output of backbone into the expansion network to
adaptively expansion $K$ (1).

$\omega_{K}^{i}=N(P(\tau,\xi)),\quad\xi\propto x-x_{r}.$ (1)

where $\tau$ refers to the feature map output by backbone, $\xi$ refers to the
normalized distance between the target pixel location $x$ and the reference
feature location $x_{r}$, $P$ refers to our expansion network, and $N$ refers
to the fact that we normalize the output of the expansion network so that the
sum is guaranteed to be $1$. So the final output $\omega_{K}^{i}$ refers to
the $K$ expansion values of the current $i$-th pixel point. Here our expansion
network uses a 5-layer MLP for expansion.

###  3.2 Multiple MLPs

The expansion network generates $K$ different adaptive scores of each element
in the feature map. After that, the predicted pixel is produced by a weighted
linear combination of $K$ MLPs with the adaptive scores. This is combined with
the $K$ designed basis, and this fusion ensures that all pixel points are the
result of the adaptive combination of the $K$ basis (2).

$SR=ReLU(\sum_{K}\tau\cdot(\omega_{K}^{i}\times f(\theta_{K}))).$ (2)

##  4 EXPERIMENTS

In this section, we describe the training details and dataset used for our
experiments. We then demonstrate the performance of our A-LIIF by comparing it
with other models.

###  4.1 Dataset

We use the DIV2K [12] dataset from the NTIRE 2017 challenge to train and
evaluate the model. In this, we use 800 images from the DIV2K [12] training
set for training and 100 images from the validation set for testing. Following
much of the previous work, Our model also compares performance on the [Set5
[14], Set14 [15], B100 [16], Urban100 [17]] benchmark dataset. Finally, we use
the widely adopted PSNR as an evaluation metric for the model.

Table 2: Quantitative comparison on benchmark datasets (PSNR (dB)). $\phi$ indicates ours implementation. EDSR-baseline [6] trains different models for different scales. Dataset | Method | In-distribution | Out-of-distribution
---|---|---|---
×2 | ×3 | 4 | ×6 | ×8
Set5 [14] | EDSR-baseline [6] | 34.99 | 31.07 | 29.14 | - | -
EDSR-baseline-LIIF [1]ϕ | 36.99 | 33.08 | 30.91 | 27.83 | 26.1
A-LIIF (ours) | 37.96 | 34.39 | 32.19 | 28.89 | 26.97
Set14 [15] | EDSR-baseline [6] | 31.41 | 28.21 | 26.59 | - | -
EDSR-baseline-LIIF [1]ϕ | 32.73 | 29.47 | 27.78 | 25.65 | 24.32
A-LIIF (ours) | 33.64 | 30.36 | 28.65 | 26.47 | 24.92
B100 [16] | EDSR-baseline [6] | 30.47 | 27.64 | 26.31 | - | -
EDSR-baseline-LIIF [1]ϕ | 31.47 | 28.49 | 27.05 | 25.40 | 24.41
A-LIIF (ours) | 32.18 | 29.11 | 27.60 | 25.84 | 24.80
Urban100 [17] | EDSR-baseline [6] | 27.91 | 25.00 | 23.59 | - | -
EDSR-baseline-LIIF [1]ϕ | 29.98 | 26.46 | 24.74 | 22.82 | 21.71
A-LIIF (ours) | 32.09 | 28.19 | 26.14 | 23.80 | 22.47

###  4.2 Training Details

In training A-LIIF, we follow the setup of the LIIF [1] model. We uniformly
down-sampled the GT HR images from the training set of 800 images in DIV2K
[12] by bicubic [6] interpolation (x2$\sim$x4) as our LR training set. For the
input model, we randomly crop the LR images into 48x48 patches, collect 2304
random pixels on the HR images, and use L1 loss to invert the ground truth.
For the choice of optimizer, we use an ADAM optimizer like LIIF [1]. The
initial learning rate is set to 10-4 and decays by half every 200 calendars.
We will train for 10-3 calendar times and perform 10-3 iterations per calendar
time. For the Encoder selection, we use EDSR [6] as our Encoder. for the
design of the expansion network, we use a 5-layer 256-neuron MLP with $K$
outputs of 10 expansion values. For the multi-MLP design, we use a 5-layer
16-neuron MLP for each pixel point. we Since our model has 16 neurons for each
pixel, we set the number of neurons in the later MLP layers to 16 when
replicating the LIIF [1].

Fig. 3: Qualitative comparison between LIIF [1] and A-LIIF.

###  4.3 Main Results

Table 1 compares the performance of LIIF [1] and our A-LIIF at different
resolution scales on DIV2K [12]. We note that the PSNR results of our model
are consistently higher than those of LIIF [1] at all scales.

We use the red box areas in Figure 1 as an example to compare LIIF and our
methods via SSIM and PSNR. The SSIM of LIIF and our method are 0.8572 and
0.8802, respectively. The PSNR of LIIF and our method are 22.1157 and 22.8543,
respectively. These results quantitatively verify that LIIF usually suffers
from distortions or ringing artifacts around edges whereas our method is not.

In Table 2, we also compare A-LIIF with other methods on fore standard
benchmark datasets: Set5 [14], Set14 [15], B100 [16] and Urban100 [17]. We
again demonstrate the superiority of our A-LIIF, and our model outperforms
both EDSR [6] and LIIF [1] in all experiments, even though EDSR [6] trains a
dedicated model for each scale. Finally, we show a qualitative comparison
between LIIF [1] and A-LIIF. as shown in Figure 3. Our A-LIIF avoids the
effect of LIIF [1] structure loss with high frequency textures by using a
fusion of expansion network and basis for each pixel.

We also compare our method on RDN [7] backbone. We only train models 100
epochs instead of 1,000 epochs. The PSNR of RDN-LIIF, and RDN-A-LIIF (Ours)
are 27.94 and 28.73, respectively. These results show that our method is still
better than LIIF [1].

###  4.4 The Ablation Study

In our ablation study, we trained several models in which we gradually added
our new design and explored the choice of hyper-parameter $K$. For a fair
comparison, we used EDSR [6] as our encoder in all models. The training setup
was kept constant, where the same ADAM optimizer, learning rate, scheduler,
and training length were used. The overall results are shown in Fig 4, Table
1-2 demonstrates that using a shared MLP layer for all pixel points for
recovery does not work well, but using an adaptive MLP for each pixel point
leads to promising results.

Fig. 4: The choice of hyper-parameter $K$.

##  5 Conclusion

In this paper, we introduced a novel arbitrary-scale super-resolution model,
dubbed A-LIIF, that combined implicit image functions with an expansion
network and the depth of multiple MLPs. It revealed the importance of using
different MLPs for different pixel points. We also verified the importance of
using non-shared MLPs by experiments.

## References

  * [1] Yinbo Chen, Sifei Liu, and Xiaolong Wang,  “Learning continuous image representation with local implicit image function,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2021, pp. 8628–8638.

  * [2] Chao Dong, Chen Change Loy, Kaiming He, and Xiaoou Tang,  “Image super-resolution using deep convolutional networks,”  IEEE transactions on pattern analysis and machine intelligence, vol. 38, no. 2, pp. 295–307, 2015.
  * [3] Jiwon Kim, Jung Kwon Lee, and Kyoung Mu Lee,  “Accurate image super-resolution using very deep convolutional networks,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 1646–1654.
  * [4] Kai Zhang, Wangmeng Zuo, Shuhang Gu, and Lei Zhang,  “Learning deep cnn denoiser prior for image restoration,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2017, pp. 3929–3938.
  * [5] Christian Ledig, Lucas Theis, Ferenc Huszár, Jose Caballero, Andrew Cunningham, Alejandro Acosta, Andrew Aitken, Alykhan Tejani, Johannes Totz, Zehan Wang, et al.,  “Photo-realistic single image super-resolution using a generative adversarial network,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2017, pp. 4681–4690.
  * [6] Bee Lim, Sanghyun Son, Heewon Kim, Seungjun Nah, and Kyoung Mu Lee,  “Enhanced deep residual networks for single image super-resolution,”  in Proceedings of the IEEE conference on computer vision and pattern recognition workshops, 2017, pp. 136–144.
  * [7] Yulun Zhang, Yapeng Tian, Yu Kong, Bineng Zhong, and Yun Fu,  “Residual dense network for image super-resolution,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2018, pp. 2472–2481.
  * [8] Jie Hu, Li Shen, and Gang Sun,  “Squeeze-and-excitation networks,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2018, pp. 7132–7141.
  * [9] Xiaolong Wang, Ross Girshick, Abhinav Gupta, and Kaiming He,  “Non-local neural networks,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2018, pp. 7794–7803.
  * [10] Tao Dai, Jianrui Cai, Yongbing Zhang, Shu-Tao Xia, and Lei Zhang,  “Second-order attention network for single image super-resolution,”  in Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, 2019, pp. 11065–11074.
  * [11] Xuecai Hu, Haoyuan Mu, Xiangyu Zhang, Zilei Wang, Tieniu Tan, and Jian Sun,  “Meta-sr: A magnification-arbitrary network for super-resolution,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2019, pp. 1575–1584.
  * [12] Eirikur Agustsson and Radu Timofte,  “Ntire 2017 challenge on single image super-resolution: Dataset and study,”  in Proceedings of the IEEE conference on computer vision and pattern recognition workshops, 2017, pp. 126–135.
  * [13] Ben Mildenhall, Pratul P Srinivasan, Matthew Tancik, Jonathan T Barron, Ravi Ramamoorthi, and Ren Ng,  “Nerf: Representing scenes as neural radiance fields for view synthesis,”  in European conference on computer vision. Springer, 2020, pp. 405–421.
  * [14] Marco Bevilacqua, Aline Roumy, Christine Guillemot, and Marie Line Alberi-Morel,  “Low-complexity single-image super-resolution based on nonnegative neighbor embedding,”  2012\.
  * [15] Roman Zeyde, Michael Elad, and Matan Protter,  “On single image scale-up using sparse-representations,”  in International conference on curves and surfaces. Springer, 2010, pp. 711–730.
  * [16] David Martin, Charless Fowlkes, Doron Tal, and Jitendra Malik,  “A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics,”  in Proceedings Eighth IEEE International Conference on Computer Vision. ICCV 2001. IEEE, 2001, vol. 2, pp. 416–423.
  * [17] Jia-Bin Huang, Abhishek Singh, and Narendra Ahuja,  “Single image super-resolution from transformed self-exemplars,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2015, pp. 5197–5206.