Continuing with experiments based on ADMM I implemented and experimented with a few different deconvolution priors. I chose several $\ell_2$ and $\ell_1$ priors based on sparse gradients, sparse curvature and simply the norm of the solution vector as well as an $\ell_{0.5}$ HyperLaplacian prior.
For all results the ground-truth signal is shown in black, the blurred and corrupted signal in red and the deblurred signal in blue. As a ground-truth signal I'm using a simple piecewise linear signal with a slanted region, some constant regions and discontinuous jumps. These features help to show the issues with different regularizers. The input signal is 128 samples long, blurred by a Gaussian filter with $\sigma=4$ and then corrupted by Gaussian noise with $\sigma=2.5%$ and $\mu=0$.
For the $\ell_2$ results I'm solving for the optimal solution directly in the Fourier domain, while for the $\ell_1$ and $\ell_0.5$ I'm using ADMM, solving the data-subproblem in the Fourier domain and the splitting variable subproblem with shrinkage/proximal operators. For all problems I've played with the prior weights a bit to try to get good results, but don't claim they are optimal.
$\ell_2$ norm of solution vector
First up is the $\ell_2$ norm of the solution. This does sharpen up the signal somewhat, but has bad ringing artifacts. These are expected for the $\ell_2$ solution, particularly given the high noise levels in the corrupted signal. The relatively large blur also leaves lots of freedom for low-frequency ringing.
$\ell_2$ norm of solution first derivative
Results for the $\ell_2$-norm of the signal derivatives do a bit better, particularly for the slanted segment, but still show low-frequency ringing. The same is true for the $\ell_2$-norm of the second derivative, although the ringing is a bit more pronounced in the flat regions:
$\ell_2$ norm of solution second derivative
Moving on to the $\ell_1$ priors, the total-variation prior is simply the $1$-norm of the signal first derivatives. It favours piecewise constant solutions which do well on the large flat regions but introduce staircase artifacts for the slanted region:
$\ell_1$ norm of solution first derivative
The total variation prior produces a much sharper signal than the $\ell_2$ priors, however for the slanted region it also introduces spurious features in the slanted region that are hard to distinguish from (correct) features elsewhere in the signal. This occurs because the assumption of the total variation prior is that the input is piecewise constant, which is not true for this signal. A better assumption is that the solution is piecewise linear, leading to the sparse-curvature prior. The sparse-curvature assumes that the second derivative of the signal is non-zero in only a few locations:
$\ell_1$ norm of solution second derivative
The sparse-curvature prior does much better in the slanted region than the total variation prior and still improves the signal, but the edges are not as sharp at discontinuities. Often people combine these priors to try to achieve a compromise between the two behaviours.
The final prior is the Hyper-Laplacian prior, which is simply the $\ell_{0.5}$-norm of the signal first derivatives:
$\ell_{0.5}$ norm of solution first derivative
The Hyper-Laplacian prior is non-convex, making it difficult to optimize for. Here it appears that the optimization got stuck: some edges are extremely sharp, however the trailing edge of the high-amplitude constant region is lost entirely, possible because of the slanted region adjacent to it. I played with the parameters quite a bit but did not get a result better than this.
Matlab example code generating these results is available at my GitHub repository:
The title says it all, I've written some basic code for loading VTK rectilinear grid files into Matlab. The code supports uniformly spaced meshes in up to four dimensions for both point and cell data.
I've posted some sample code on Github for performing image deblurring in Matlab using Mex. I wrote it as a way to play around with the ADMM algorithm for sparse signal reconstruction, as described in Stephen Boyd's ADMM paper, as well as to get some experience using C++ code from Matlab.
The code uses matrix-free linear operators to evaluate the forward and reverse measurement (blurring) process, solving the data-term subproblem with gradient descent. While you can do this much faster using FFTs for deconvolution, the code should generalize quite well to tomography with non-orthographic projections. I've tried to comment it well so that I can figure out what I did if I ever need to pick it up again.
In the image above, the left plot is the ground-truth image and the right is blurred by a point-spread function and corrupted by Gaussian noise. After running the demo, the following results are obtained:
The reconstructed image is shown on the right, while a slice through the center scanline shows the reconstruction (blue dots), ground-truth (solid black line) and input (red line). From the results it's clear that the method does well at preserving discontinuities even for fairly large blurs and noise levels.
The demo is pretty quick to run and is helpful for choosing parameters since you can see how the behavior of the method changes with the regularizer weight and the penalty term weight.
A recent project focused on reducing the cost of Transient Imaging. Transient images are effectively videos that depict light-transport at the timescales of light-transport; they allow you to see light 'in-flight' as it bounces around the scene.
This has been done in the past with extremely expensive optical setups costing hundreds of thousands of dollars that involve very short laser pulses and so-called streak-cameras. However, it is possible to accomplish much the same effect with modifications to off-the-shelf hardware for on the order of only a $1000. The paper describing this has been accepted to SIGGRAPH 2013 and there will be an accompanying technical demo in the Emerging Technologies area.
You can get an idea what the project does from the following video. If you happen to be attending SIGGRAPH this year, please stop by the booth to find out more!
My most recent submission on image deblurring has been accepted as a poster at CVPR2013, to be held in Portland Oregon. The approach is an adaptation of the stochastic random walk used in my earlier Stochastic Tomography paper to image deblurring and handles saturated pixels and boundary conditions naturally. The submission video gives an overview of the method:
More info is available at the Stochastic Deconvolution website. The site also includes sample C++ code illustrating the basic method.
Implementing Bresenham's line drawing algorithm is a pain and has some drawbacks for motion control applications. In particular, it relies on swapping endpoints of the line-segments to achieve specific preconditions and has eight configurations (in 2D alone!) that must be implemented to draw arbitrarily oriented lines. This has big disadvantages for implementation in high dimension (6D in my case) and even bigger disadvantages for motion control, where you want to move from point A to point B, not arbitrarily (and incorrectly) assume you're already at B and then move backward to A because, hey, order doesn't matter!
I recently ran across Cris Luengo's blog post on a Bresenham-equivalent algorithm in arbitrary dimension based on floating point rounding. This is likely to be faster on modern machines than the integer only version, simply because i) floating point ops and casting are cheap and ii) branches are expensive, which is effectively the opposite of the case when the original Bresenham algorithm was developed. He claims the method reproduces the original Bresenham algorithm's plotted pixels, although in my tests I seem to get occasional errant pixels even with double precision. It may well be an error on my part.
Regardless, this post gave me the idea of just implementing Bresenham's algorithm in a parametric form, accumulating error and advancing pixels in both (all) coordinates in the same manner. This is considerably more costly in 2D, where you double the operations required, but the relative cost decreases as the dimension increases. Since I'm targeting 6D, the relative cost is only 1/5th, and the simplicity of implementation and lack of endpoint switching more than makes up for it.
I've put up simple C++ code for 2D and 3D integer-only line drawing using a test implementation that is freely available for any use. It draws 50 random lines with the integer only algorithm described above into the green channel and the same 50 lines with the floating point rounding method (in double precision) into the red channel. The results are below, if the algorithms matched exactly, all the lines would be yellow, with no red or green pixels anywhere.
It looks okay at first, but if you zoom in and look carefully you can see a few red and green pixels around. Regardless, they both do a good job of plotting lines, and if my code is off by a step or two occasionally it's not the end of the world.
I frequently use CImg to perform basic image operations like loading, saving and resizing. However trying to debug code using the CImg library under XCode has been problematic, since i) CImg uses ImageMagick executables to perform format conversions and ii) XCode doesn't check the system $PATH variable when debugging.
To get around this, you can explicitly specify the path to the ImageMagick executables in your code. This avoids the command-not-found error that otherwise occurs when using CImg through the debugger. You can dink around with trying to bend XCode to your will, but I've found this is just way easier.
Here's the related code, note you will need to change the path to suit your system:
As noted in a previous post, my Stochastic Tomography paper was accepted to SIGGRAPH 2012. Last Tuesday I was in Los Angeles for the conference to present the paper, including present the synopsis in the conference 'fast-forward' to an appallingly large audience. The photo below shows the seating but during the actual event it was standing room only.
A bit nerve-wracking to say the least. However it went well and after presenting my main talk to a MUCH smaller crowd, I'd like to post some photos of the setup that we used for the paper. I should point out that this project actually did not contribute much to the capture setup, this was previously in place from work done by Borislav Trifonov, Michael Krimerman, Brad Atcheson, Derek Bradley and a slew of others. My work on this paper focused on the algorithms primarily, but I thought people might be interested in a quick overview of the tomography capture apparatus.
The goal of the paper was to build 3D animated models of mixing liquids from multiview video. To accomplish this, we used an array of roughly 16 Hi-Def Sony camcorders arranged in a semi-circle around our capture volume to record video streams of the two fluids mixing.
You can see the cameras in the photo above, all focused on the capture volume which is inside the glass cylinder. Each of these records a video of the mixing process, producing 16 streams of video that look more or less like the photo shown below:
You can see one of the cameras peeking out at the right side of the frame. The cameras are controlled by an Arduino based box that talks to each camcorder using the SONY LANC protocol. This is an unpublished protocol used by SONY editing consoles, however it has been reverse-engineered by others to allow people to control SONY equipment. We implemented this protocol on an Arduino, which allows us to start all the cameras recording, turn them on and off as arrays, switch them between photo and video mode and so on. Unfortunately we can't easily set the exposure levels, transfer files to and from the device, instead we have to painstakingly do this by hand through the on-camera menus, which is error-prone and time-consuming.
The two fluids we use are water, for the clear liquid, and Fluoroscein-Sodium fluorescing dye for the mixing liquid. This fluorescent dye is available in a powder that is soluble in water, which allows us to perform several types of capture. The image above is dye powder dropped onto the surface of the water, this mixes with the water and is slightly denser, forming the Rayleigh-Taylor mixing process you see in that shot. We can also pre-mix the dye powder and simple pour or inject it into the domain, this was the process used for the following two captures that were used in the paper.
This shows an unstable vortex propaging downwards, leaving a complex wake. I recommend watching in high-def (720p or 1080p). The next is alcohol mixed with the dye powder, injected into the cylinder from the bottom. Since alcohol is less dense than water, it rises under buoyancy, mixing as it goes.
In the video above you can see a laminar to turbulent transition as well as lots of complex eddies that form as part of the mixing process.
The captures are illuminated with a set of white LED concert strobe panels. These panels serve two purposes. First they let us get LOTS of light into the scene in a controlled fashion. Second we actually use a strobed illumination at about 30Hz to optically synchronize the cameras and remove rolling-shutted shear effects.
All captures start in darkness so we can tell the time offset from the start of the video to the first frame where there is significant illumination. In fact we can do better than alignment to a single frame, since with the rolling shutter used by these cameras, we can actually determine the first scanline that is exposed. Using a 30Hz illumination pattern, we can also determine the exposure setting of the camera by looking for the last scanline before the light goes off again.
We then have a rolling shutter compensation program that scans through each video and reassembles a new video from the exposed and dark scanlines. The result is a set of videos that are optically synchronized and that have minimal shearing introduced.
This gives us a set of input data, however we also need to perform some geometric calibration of the scene in order to know from what angle each video was recorded and to be able to obtain the ray inside of the capture volume that corresponds to every observed pixel.
To do this, we use an optical calibration library called CalTag that detects self-identifying marker patterns similar to QR codes in images. We print a calibration target on overheard transparencies and mount this pattern to a 3D printed calibration jig that is placed in the glass cylinder.
This jig fits tightly in the cylinder and is registered with a set of detents that fit into recesses in a registration plate that is glued to the inside of the capture cylinder. The marker pattern that you see in the photo above is also registered to a set of registration tabs. We have a calibrated pattern on the front of the target as shown above, but also on the back.
When a camera takes an image of this jig after placing it into the capture domain (filled with water), an image similar to the following is obtained, although generally with far less blur due to condensation.
CalTag can then give us the corners of the marker patterns, which can be interpolated to associate with every image pixel, the corresponding 3D point on the calibration target that it 'sees'. We then rotate the target 180 degrees to face away from the camera and take a picture of an identical and carefully aligned target on the back side of the jig, giving another 3D point for each pixel. Connecting the points gives a ray in 3D space inside the cylinder, without having to account for any optical interactions between the interior liquid and cylinder.
We do this for every camera, by mounting the capture domain on a rotation stage, which is again controlled by an Arduino. An automated calibration procedure rotates the stage and triggers each camera to image the front calibration plane, then rotates an additional 180 degrees to repeat the process. The whole mess is controlled by a python script using pySerial, including the strobes, the rotation stage and the embedded camera controller.
This gives us the needed calibration data to express our scene as a tomographic inverse problem. Here we look for the scene content that would reproduce the measurements (videos) we obtained, given a physical model for the scene. In this capture case, the scene is simply emissivity adding up along a ray-path, so we get a linear inverse problem, that we solve using our new Stochastic Tomography algorithm. The result is volumetric 3D fields that you can animate, inspect and slice through and re-render however you like, as seen below in the submission video.
My most recent project has just been accepted to SIGGRAPH 2012! The paper is called Stochastic Tomography and its Applications in 3D Imaging of Mixing Fluids, more information (including the paper, submission video and datasets) can be found on the project site: http://www.cs.ubc.ca/labs/imager/tr/2012/StochasticTomography/
The paper presents a new method for tomographic reconstruction based on a point-sampling process similar to the Metropolis-Hastings algorithm. In the paper this is applied to building reconstructions of a fluorescent dye mixing with alcohol from multi-view video. The algorithm is well suited to this task since it allows arbitrary projections and is naturally adaptive. It is also very straightforward to incorporate effectively any regularizer into the algorithm, allowing application specific prior-knowledge to be exploited during the reconstruction process.
I hope to post a bit more about the hardware setup that is used for this project in the new week or so.
Recently I had need for code to interpolate data stored on regular grids in 1D, 2D and 3D. This code is straightforward but time-consuming to write and I often find myself re-implementing it for each project that I work on.
The linked header file includes templated interpolation routines for 1D, 2D and 3D interpolation on uniformly gridded data with linear, cosine and cubic (Catmull-Rom) interpolants. It is fully documented with DOXYGEN and provides a simple fairly consistent API. The code works by interpolating values that are passed in either as individual values or as arrays allowing it to work with arbitrarily defined grids and also optionally returns the weights associated with each value. I have used the code to implement
The following code snippet taken from one of the 2D test functions shows the use of the API. In all API functions, points are ordered lexographically by axis, i.e. p[0] = [x, y], p[1] = [x+1, y], p[2] = [x, y+1], p[3] = [x+1, y+1].
/**
@brief test function for the 2D linear interpolation routines, currently just
checks that the methods are consistent, i.e. that they all produce the same
output for the same input
@return true if test passes false otherwise
*/
static bool interpolation_test_lerp_2D(){
const double t0=0.325f, t1=0.675, v[] = { 1.0, 2.0, 2.0, 5.0 };
double o0, o1, o2, o3, w[4];
bool result = true;
// check that API functions are consistent
o0 = interpolation_lerp_2D( v[0], v[1], v[2], v[3], t0, t1 );
o1 = interpolation_lerp_2D( v[0], v[1], v[2], v[3], t0, t1, w[0], w[1], w[2], w[3] );
o2 = interpolation_lerp_2D( v, t0, t1 );
o3 = interpolation_lerp_2D( v, t0, t1, w );
result &= fabs(o0-o1) < 1e-9 && fabs(o0-o2) < 1e-9 && fabs(o0-o3) < 1e-9;
// check that input values are reproduced
result &= fabs( interpolation_lerp_2D( v[0], v[1], v[2], v[3], 0.0, 0.0 )-v[0] ) < 1e-9;
result &= fabs( interpolation_lerp_2D( v[0], v[1], v[2], v[3], 1.0, 0.0 )-v[1] ) < 1e-9;
result &= fabs( interpolation_lerp_2D( v[0], v[1], v[2], v[3], 0.0, 1.0 )-v[2] ) < 1e-9;
result &= fabs( interpolation_lerp_2D( v[0], v[1], v[2], v[3], 1.0, 1.0 )-v[3] ) < 1e-9;
return result;
}
The code is freely available for all purposes and can be downloaded below:
When working with calibrated cameras it is often useful to be able to display things on screen for debugging purposes. However the camera model used by OpenGL is quite different from the calibration parameters from, for example, OpenCV. The linear parameters that OpenCV provides are the following:
where (from http://en.wikipedia.org/wiki/Camera_resectioning) is the skew between the x and y axes, are the image principle point , with f being the focal length and being scale factors relating pixels to distance. Multiplying a point by this matrix and dividing by resulting z-coordinate then gives the point projected into the image.
The OpenGL parameters are quite different. Generally the projection is set using the glFrustum command, which takes the left, right, top, bottom, near and far clip plane locations as parameters and maps these into "normalized device coordinates" which range from [-1, 1]. The normalized device coordinates are then transformed by the current viewport, which maps them onto the final image plane. Because of the differences, obtaining an OpenGL projection matrix which matches a given set of intrinsic parameters is somewhat complicated.
Roughly following this post, (update: a much-improved update from Kyle, the post's author is available here) the following code will produce an OpenGL projection matrix and viewport. I have tested this code against the OS-X OpenGL implementation (using gluProject) to verify that for randomly generated intrinsic parameters, the corresponding OpenGL frustum and viewport reproduce the x and y coordinates of the projected point. The code works by multiplying a perspective projection matrix by an orthographic projection to map into normalized device coordinates, and setting the appropriate box for the glViewport command.
/**
@brief basic function to produce an OpenGL projection matrix and associated viewport parameters
which match a given set of camera intrinsics. This is currently written for the Eigen linear
algebra library, however it should be straightforward to port to any 4x4 matrix class.
@param[out] frustum Eigen::Matrix4d projection matrix. Eigen stores these matrices in column-major (i.e. OpenGL) order.
@param[out] viewport 4-component OpenGL viewport values, as might be retrieved by glGetIntegerv( GL_VIEWPORT, &viewport[0] )
@param[in] alpha x-axis focal length, from camera intrinsic matrix
@param[in] alpha y-axis focal length, from camera intrinsic matrix
@param[in] skew x and y axis skew, from camera intrinsic matrix
@param[in] u0 image origin x-coordinate, from camera intrinsic matrix
@param[in] v0 image origin y-coordinate, from camera intrinsic matrix
@param[in] img_width image width, in pixels
@param[in] img_height image height, in pixels
@param[in] near_clip near clipping plane z-location, can be set arbitrarily > 0, controls the mapping of z-coordinates for OpenGL
@param[in] far_clip far clipping plane z-location, can be set arbitrarily > near_clip, controls the mapping of z-coordinate for OpenGL
*/
void build_opengl_projection_for_intrinsics( Eigen::Matrix4d &frustum, int *viewport, double alpha, double beta, double skew, double u0, double v0, int img_width, int img_height, double near_clip, double far_clip ){
// These parameters define the final viewport that is rendered into by
// the camera.
double L = 0;
double R = img_width;
double B = 0;
double T = img_height;
// near and far clipping planes, these only matter for the mapping from
// world-space z-coordinate into the depth coordinate for OpenGL
double N = near_clip;
double F = far_clip;
// set the viewport parameters
viewport[0] = L;
viewport[1] = B;
viewport[2] = R-L;
viewport[3] = T-B;
// construct an orthographic matrix which maps from projected
// coordinates to normalized device coordinates in the range
// [-1, 1]. OpenGL then maps coordinates in NDC to the current
// viewport
Eigen::Matrix4d ortho = Eigen::Matrix4d::Zero();
ortho(0,0) = 2.0/(R-L); ortho(0,3) = -(R+L)/(R-L);
ortho(1,1) = 2.0/(T-B); ortho(1,3) = -(T+B)/(T-B);
ortho(2,2) = -2.0/(F-N); ortho(2,3) = -(F+N)/(F-N);
ortho(3,3) = 1.0;
// construct a projection matrix, this is identical to the
// projection matrix computed for the intrinsicx, except an
// additional row is inserted to map the z-coordinate to
// OpenGL.
Eigen::Matrix4d tproj = Eigen::Matrix4d::Zero();
tproj(0,0) = alpha; tproj(0,1) = skew; tproj(0,2) = u0;
tproj(1,1) = beta; tproj(1,2) = v0;
tproj(2,2) = -(N+F); tproj(2,3) = -N*F;
tproj(3,2) = 1.0;
// resulting OpenGL frustum is the product of the orthographic
// mapping to normalized device coordinates and the augmented
// camera intrinsic matrix
frustum = ortho*tproj;
}
The code uses the Eigen linear algebra library, which conveniently stored matrices in column-major order, so applying the resulting frustum matrix is as simple as:
I recently got to visit a local studio and while I was there they asked if I wanted my face 3D scanned. Of course I did. Who wouldn't? Well at least one person. But I did. And here's the result, after runing it through the quadric-collapse decimator of Meshlab.
It's pretty damn good quality wise, even on challenging areas like my beard. Now I just have to merge the scans and create a watertight mesh, ready for 3D printing!
is the skew between the x and y axes, 
are the image principle point 
, 
with f being the focal length and 
being scale factors relating pixels to distance. Multiplying a point 
by this matrix and dividing by resulting z-coordinate then gives the point projected into the image.
