import numpy as np
%matplotlib notebook
import matplotlib.pyplot as plt
import traceback

def circle(x, y, r, rot, steps=100):
    a = rot + np.linspace(0., 2. * np.pi, steps)
    return np.append((np.cos(a) * r) + x, x), np.append((np.sin(a) * r) + y, y)

class Element:
    
    def __init__(self, rate, edge_radius, path_radius, path_offset, tracers):
        self.centre_x     = 0
        self.centre_y     = 0
        self.path_x       = 0
        self.path_y       = 0
        self.edge_radius  = edge_radius
        self.edge_handle, = plt.plot([],[], '-')
            
        self.path_radius  = path_radius
        self.path_offset  = path_offset
        self.path_handle, = plt.plot([],[], ':')
        
        self.tracers = tracers
        self.tracer_handles = []
        self.tracer_paths   = []
        for colour, fmt, _ in self.tracers:
            hnd, = plt.plot([],[], fmt, color=colour)
            self.tracer_handles += [ hnd ]
            self.tracer_paths   += [ [[], []] ]
        
        self.parent    = None
        self.children  = []
        self.rate      = rate
        self.position  = 0
        self.rotation  = 0
        self.acc_rotation = 0
    
    def step(self, show_elements):
        
        self.rotation += self.rate
        
        if (self.parent is not None): 
            self.acc_rotation = self.rotation + self.parent.acc_rotation
            
            if (self.parent.path_radius <= self.parent.edge_radius):
                offset = self.parent.path_radius - self.edge_radius
                delta  = (self.rate * self.edge_radius) / self.parent.path_radius
            else:
                offset = self.parent.path_radius + self.edge_radius
                delta  = -((self.rate * self.edge_radius) / self.parent.path_radius)
                
            self.position -= delta*2
        
            self.centre_x = self.parent.path_x + np.cos(self.position) * offset
            self.centre_y = self.parent.path_y + np.sin(self.position) * offset
        
        if show_elements:
            self.edge_handle.set_data(*circle(self.centre_x, 
                                              self.centre_y, 
                                              self.edge_radius, 
                                              self.acc_rotation))
        else:
            self.edge_handle.set_data([],[])
        
        self.path_x = self.centre_x + np.cos(self.acc_rotation) * self.path_offset
        self.path_y = self.centre_y + np.sin(self.acc_rotation) * self.path_offset
        
        if show_elements:
            self.path_handle.set_data(*circle(self.path_x, 
                                              self.path_y, 
                                              self.path_radius, 
                                              self.acc_rotation))
        else:
            self.path_handle.set_data([],[])
        
        for hnd, path, (_, _, offset) in zip(self.tracer_handles, 
                                                    self.tracer_paths, 
                                                    self.tracers):
            
            path[0] += [ self.path_x + np.cos(self.acc_rotation) * offset ]
            path[1] += [ self.path_y + np.sin(self.acc_rotation) * offset ]
            hnd.set_data(*path)
            
        for child in self.children:
            child.step(show_elements)
            
    def add(self, element):
        self.children += [ element ]
        self.children[-1].parent = self
    
class UI:
    def __init__(self):
        
        self.show_elements = True
        figw, figh, figdpi = 950, 950, 50
        self.figure = plt.figure(facecolor='w', figsize=(figw/figdpi, figh/figdpi), dpi=figdpi)
        self.axis   = plt.gca()
        plt.axis('off')
        plt.xlim([-1,1])
        plt.ylim([-1,1])
        
        radius = 0.45
        num_pens = 10
        cmap = [plt.get_cmap('RdBu')(int(idx)) 
                for idx in np.linspace(64, 256-64, num_pens)]
        
        pens = [(colour, '-', radius-delta) 
                for colour, delta in zip(cmap, np.linspace(0.1,0.3,num_pens))]
                
        elem_0            = Element(0.05, radius, radius, 0.0, pens)
        self.root_element = Element(0.0, 1.0, 1.0, 0.0, [])
        self.root_element.add(elem_0)
        
        plt.show()
    
    def on_click(self, event):
        if not (event.inaxes == self.axis): return
        plt.sca(self.axis)
        
        try:
            self.show_elements = not self.show_elements
            self.root_element.step(self.show_elements)      
        except Exception:
            plt.title(traceback.format_exc())
        
    def on_move(self, event):
        if not (event.inaxes == self.axis): return
        plt.sca(self.axis)
        
        try:
            self.root_element.step(self.show_elements)      
        except Exception:
            plt.title(traceback.format_exc())
        
    def attach(self, key, func):
        self.figure.canvas.mpl_connect(key, func)

ui = UI()

def on_click(event):
    global ui
    ui.on_click(event)
    
def on_move(event):
    global ui
    ui.on_move(event)
        
ui.attach('button_press_event',  on_click)
ui.attach('motion_notify_event', on_move)

A video of my talk taken by my co-author Prof. Mark Jones on our paper "A Deep Learning Approach to No-Reference Image Quality Assessment For Monte Carlo Rendered Images" published at CGVC 2018.

In Full-Reference Image Quality Assessment (FR-IQA) images are compared with ground truth images that are known to be of high visual quality. These metrics are utilized in order to rank algorithms under test on their image quality performance. Throughout the progress of Monte Carlo rendering processes we often wish to determine whether images being rendered are of sufficient visual quality, without the availability of a ground truth image. In such cases FR-IQA metrics are not applicable and we instead must utilise No-Reference Image Quality Assessment (NR-IQA) measures to make predictions about the perceived quality of unconverged images. In this work we propose a deep learning approach to NR-IQA, trained specifically on noise from Monte Carlo rendering processes, which significantly outperforms existing NR-IQA methods and can produce quality predictions consistent with FR-IQA measures that have access to ground truth images.

Cite as:

@inproceedings{Whittle2018, 
    author    = "Whittle, Joss and Jones, Mark W.", 
    booktitle = "Computer Graphics and Visual Computing (CGVC) 2018", 
    title     = "A Deep Learning Approach to No-Reference Image Quality Assessment For Monte Carlo Rendered Images", 
    year      = "2018", 
    month     = "Sep",
    pages     = {23--31}, 
    doi       = {10.2312/cgvc.20181204} 
}

Initial work towards the creation of a CNC modification of a classic Etch-A-Sketch toy.

The stepper motors used for this project were super cheap, just 12 quid on amazon for 5 small 5v stepper motors which each came with their own control board allowing me to control them over PWM from a Raspberry Pi Model 3.

I created 3D printed parts for mounting the stepper motors to the Etch-A-Sketch and used a laser cutter to create gears out of 6mm birch plywood.

Being cheap and low powered the device suffers from being quite slow which is an issue for the Etch-A-Sketch. Slow movement tends to cause material build up on the drawing cursor making the drawn line become increasingly wide as drawing progresses.

With some software tweaking I think some of these issues can be alleviated which should allow for the project to be demonstrated at department open days for prospective students considering studying at Swansea.

This project was inspired by the fantastic pen plotter visualizations created by Michael Fogleman.

Poisson Disk Sampling is a technique for drawing batches of blue noise distributed samples from an n-dimensional domain. The method works by selecting an initial, seed, point and proposing k (the branching factor) random points within 1 to 2 radii r of the initial point. For each of these proposal points we test whether they are closer than the threshold radius r from any of the accepted points (initially just the seed point). If the point is far enough away from any of the accepted points, it becomes accepted, and we sample k new points around it for later processing. If the point is too close to another accepted point it is immediately discarded.

By continuing the process for a given number of accepted points, until the n-dimensional space can no longer be filled without points being closer than the threshold r, or some other criteria is met we end up with a set of well distributed points that have nice mathematical properties when used in stochastic approximation methods.

An interesting observation is that the set of sample points is "grown" outwards from the seed point, and that each accepted point can trace its origin to a single parent point which spawned it. If we connect the sampled points as a tree hierarchy we can visualize the growth pattern of sample set as a tree.

The implementation I used to generate the above image used the sampling algorithm described in Fast Poisson Disk Sampling in Arbitrary Dimension, Robert Bridson 2007 which can produce batches of well distributed samples in O(n) computation time.

I have released the code for this project open source as a jupyter notebook. The main bottleneck in this code is actually the line plotting of the sample tree due to limitations of Matplotlib. With a better method of drawing the generated trees larger and deeper growth patterns could easily be visualized.

In this project I implemented the paper Conditional Image Synthesis With Auxiliary Classifier GANs, Odena et. al. 2016 using the Keras machine learning framework.

Generative Adversarial Networks, Goodfellow et. al. 2014 represents a training regime for teaching neural networks how to synthesize data that could plausibly have come from a distribution of real data - commonly images with a shared theme or aesthetic style such as images of celebrity faces (CelebA), of handwritten digits (MNIST), or of bedrooms (LSUN-Bedroom).

In GANs two models are trained - a generative model that progressively learns to synthesize realistic and plausible images from a random noise input (the latent vector) - and a discriminative model that learns to tell these generated (fake) images from real images sampled from the target dataset. The two models are trained in lock-step such that the generative model learns to fool the discriminator model, and the discriminator adapts to become better at not being fooled by the generator.

This forms a minimax game between the two models which converges to a Nash equilibrium. At this point the generator should be able to consistently produce convincing images that appear to be from the original dataset, but are in-fact parameterized by the latent vector fed to the generative model.

Auxiliary Classifier GANs extend the standard GAN architecture by jointly minimizing the generators ability to fool the discriminative model, with the ability of the discriminator to correctly identify which digit it was shown. This allows the generative model to be parameterized not only by a random latent vector, but also a representative encoding of which digit we would like it to synthesize.

The above image shows the result of my AC-GAN implementation trained on the MNIST dataset. On the left we see real images sampled randomly from MNIST for each of the 10 digit classes, and on the right we see images synthesized by the generative model for each class. The generated images are not sampled completely randomly, in this image I was selecting a random value of the latent vector and sweeping it from a value of 0 to 1. We can see that for each digit class the had the subtle effect of adjusting rotation and "flair" or perhaps "serif-ness", showing that the generative model has mapped the space of possible values that exist in the latent vector to different stylistic traits of the produced digits.

The results of this experiment are satisfying but not great overall. I believe the model suffers from, at least partial, "mode collapse" where the generator learns to produce a subset of possible stylistic variations convincingly and so never attempts to learn how to produce other stylistic variants.

Since the publication of Goodfellow's seminal work on GANs many variations have been proposed that attempt to solve common issues such as mode collapse and training stability.

In the future I plan to revisit this project and implement some of the newer and more advanced methods. While the code for this project is written as a jupyter notebook I do not plan to release the code as it is not very clean or well documented. I will however release well documented code when I revisit this project.

This image was generated using a renderer developed during my PhD using the bidirectional path tracing algorithm.

Here I simulate the interactions of light coming from a diffuse area light source, scattering through a glass model of the Stanford Bunny which has been modified to have both inner and outer walls, and filled with a simulated wine like liquid. In the foreground, within the shadow of the bunny, we can see caustic illumination patterns where light has been tinted and focused as it undergoes refraction as it passes through the glass and liquid mediums. In the background of the image the grid pattern on the floor becomes blurred and out of focus due to the physical simulation of light interactions with a camera aperture and lens elements. Similarly, the aperture simulation can be seen in the form of small hexagonal specular highlights on the bunnies ears which occur due to the aperture being modeled as a six sided polygon.

This image was generated using a renderer developed during my PhD using the bidirectional path tracing algorithm.

In ray tracing based rendering algorithms the main performance bottleneck comes from the time it takes to perform intersection tests between rays and the scene geometry. As the number of elements used to represent surfaces within the scene increases so to does the computation complexity of performing these intersection tests.

Rather than testing a given ray against each surface element present in the scene, which would scale poorly with O(n) complexity, acceleration structures such as the KD-Tree and BHV-Tree can be used to reduce the number of intersection tests to roughly O(log n) on average. This is accomplished by traversing a binary-tree like hierarchy performing a comparatively cheap intersection test between a ray and a cutting plane or axis-aligned bounding-box. When the leaf nodes of the tree have been reached only the mesh elements contained within the leaf node need to be tested against the given ray. After an initial intersected element has been found all nodes of the tree which can only contain elements that would have to be further away than the current closest intersection can be ignored as the tree traversal continues, yielding significant performance gains.

The Stanford Dragon model shown in this image contains 871,414 triangular mesh elements. In the above image, 99 copies of the dragon are present. By a naive scene construction method this would total 86,269,986 triangles. This has the effect of greatly increasing the memory requirements and size of the resulting BVH-Tree. In an ideal world we would only store the dragon mesh and a BVH-Tree over its elements once, and have a system for replaying this data efficiently in arbitrary positions, rotations, and scales.

To this end I have implemented geometric instancing as a feature within the rendering software developed for my PhD research. This is done through a coordinate-space transformation and cross-linking injected into the acceleration structure traversal. The dragon mesh is loaded into the renderer, and a BVH-Tree is constructed for the single dragon model. The root node of the BVH-Tree's bounding box is then duplicated through independent affine-transformations into the size, rotation, and position of each of the desired dragons. Another BVH-Tree is then constructed over the these bounding boxes, along with the triangles making up the floor, walls, and light source of the scene. During traversal of the BVH-Tree, when one of the transformed bounding boxes is queried, the given ray is transformed by the inverse of the affine-transformation matrix for the current bounding box. This inverse transformed ray is then intersected with the BVH-Tree containing only the dragon mesh. The resulting intersection location can be found by scaling the original ray by the resulting distance to the closest intersection found, regardless of whether the intersection was within one of the nested bounding boxes.

With some additional work, each instanced mesh can have different materials applied to it through a stack based shading hierarchy which is also tracked during BVH-Tree traversal.

This image was generated using a renderer developed during my PhD using the bidirectional path tracing algorithm in submission of the 2015, SURF Research as Art Competition, hosted for reasearcher at Swansea University.

The image is rendered at 4k resolution and was rendered in approximately 8 hours to a high degree of convergence by utilizing 50 of the machines in the Swansea Computer Science Department's linux teaching lab as a computational cluster using the Message Passing Interface (MPI) framework to implement multi-node parallism into my rendering software.

This image was generated using a renderer developed during my PhD using the bidirectional path tracing algorithm.

The pig model is shaded with my interpretation of a glazed ceramic material. This combines a diffuse ceramic substrate with a glossy overlayer which are Fresnel blended during shading computation. The glaze overlay has a normal map applied to it for adding detailed surface variation.

This image was generated using a renderer developed during my PhD using the bidirectional path tracing algorithm.

Just a quick render after implementing an Anisotropic Metal material in my renderer. The test scene was inspired by one I saw on Kevin Beason's worklog blog.

Yesterday I stumbled upon a lesser known and far superior method for mapping points from a square to a disk. The common approach which is presented to you after a quick google for Draw Random Samples on a Disk is the clean and simple mapping from Cartesian to Polar coordinates; i.e.

Given a disk centered origin (0,0) with radius r

// Draw two uniform random numbers in the range [0,1)
R1 = RAND(0,1);
R2 = RAND(0,1);

// Map these values to polar space (phi,radius)
phi = R1 * 2PI;
radius = R2 * r;

// Map (phi,radius) in polar space to (x,y) in Cartesian space
x = cos(phi) * radius;
y = sin(phi) * radius;

The result of this sampling on a regular grid of samples is shown in the image below. The left plot shows the input points as simple ordered pairs in the range [0,1)^2, while the right plot shows these same points (colour for colour) mapped onto a unit disk using Polar mapping as described above.

As you can see the mapping is not ideal with many points being over-sampled at the poles (I wonder why they call is Polar coordinates), and with areas towards the radius left under-sampled. What we would actually like is a disk sampling strategy that keeps the uniformity seen in the square distribution while mapping the points onto the disk.

Enter, Concentric Disk Sampling. A Low Distortion Map Between Disk and Square, Shirley & Chiu 1997 presents the idea for warping the unit square into that of a unit circle. Their method is nice but it contains a lot of nested branching for determining which quadrant the current point lays within. Shirley mentions an improved variant of this mapping on his blog, accredited to Dave Cline. Cline's method only uses one if-else branch and is simpler to implement.

Again, given a disk centered origin (0,0) with radius r

// Draw two uniform random numbers in the range [0,1)
R1 = RAND(0,1);
R2 = RAND(0,1);

// Avoid a potential divide by zero
if (R1 == 0 && R2 == 0) {
    x = 0; y = 0;
    return;
}

// Initial mapping
phi = 0; radius = r;
a = (2 * R1) - 1;
b = (2 * R2) - 1;

// Uses squares instead of absolute values
if ((a*a) > (b*b)) { 
    // Top half
    radius  *= a;
    phi = (pi/4) * (b/a);
}
else {
    // Bottom half
    radius *= b;
    phi = (pi/2) - ((pi/4) * (a/b)); 
}

// Map the distorted Polar coordinates (phi,radius) 
// into the Cartesian (x,y) space
x = cos(phi) * radius;
y = sin(phi) * radius;

This gives a uniform distribution of samples over the disk in Cartesian space. The result of the mapping applied to the same set of uniform square samples is shown above. Notice how we now get full coverage of the disk using just as many samples, and that each point has (relatively) equal distance to all of it's neighbors, meaning no bunching at the poles, and no under-sampling at the fringe.

I've applied this sampling technique to my path tracer as a means of sampling the aperture of the virtual camera when computing depth of field. Convergence to the true out-of-focus light distribution is now much faster and more accurate than it was with Polar sampling which, due to oversampling at the poles, cause a disproportionate number of rays to be fired along trajectories very close to the true ray.