PengStarobinets

Back to Psych 221 Projects 2013

Background

The performance of long-distance imaging systems can often be strongly affected by atmospheric turbulence caused by variation of refractive index along the optical transmission path. Such turbulence can produce geometric distortion, space and time-variant defocus blur, and motion blur. An example is shown in the following video of the moon.

Below is the exmaple of atmospheric turbulence that was used for this project:

Recording of the moon with effects of turbulence

There have been many approaches to solving this problem that attempt to restore a single high-quality image from an observed frame sequence distorted by air turbulence. As in the video, these approaches, and the approach addressed in this paper, work under the assumption that the scene and the image sensor are both static and that observed motions are due to the air turbulence alone.

The imaging process can be modeled as:

$g_{k}=H_{k}F_{k}f+n_{k}$ where $f$ denotes the ideal image, $F_{k}$ and $H_{k}$ represent the geometric deformation and blurring matrices respectively, $n_{k}$ denotes additive noise, and $g_{k}$ is the k-th observed frame.

The key then becomes to basically reverse this process so that we can find the desired corrected image.

Methods

Existing restoration algorithms for this problem can generally be categorized in two ways.

Multi-Frame Reconstruction Framework

First, a nonrigid image registration technique is used to register each observed frame with respect to a fixed reference grid (usually generated from an average of the image sequence to eliminate spatial variation). Next, the registration parameters are used to estimate the corresponding motion field for each frame and one frame is transformed back to a non-deformed position. Finally, a sharp image is formed through a Bayesian reconstruction filter.

Lucky Imaging

Image selection and fusion methods are used to reduce the blurring effects caused by turbulence. The image selection algorithm attempts to find frames of the best quality (lucky frames) from the image stream. The output image is then produced by fusing these best-quality images together. The principle behind this method is that for short-exposuer images, turbulence creates 'mutations' in image quality and randomly makes some images sharp. Lucky-image-based weighting schemes can be used to select the sharpest image or images with certain sharp regions, and combine them to create an image that is sharp everywhere.

Our Method

To correct geometric distortion and reduce space and time-varying blur, a new approach is proposed in this paper capable of restoring a single high-quality image from a given image sequence distorted by atmospheric turbulence. This approach uses two major steps. For the first step, we use a B-spline based non-rigid image registration algorithm to register each observed frame with respect to a reference image. To improve registration accuracy, a symmetry constraint is introduced to the cost function to penalize inconsistency between the forward and backward deformation parameters during the estimation process. A fast Gauss-Newton method is used to reduce the computational cost of the registration algorithm. We improved the computational efficiency further by applying our method to multiple small blocks and stitching them together. In addition to reducing computational load and time, this improves on the original method by reducing the effect of systematic errors on the system, because an error will be isolated to the particular block rather than affecting the whole image. In the second step, we restore a high quality image from the registered observed frames using a high-pass sharpening filter.

Proposed Approach

Non-rigid Image Registration

First, we obtain a reference image without image deformation, generally by averaging the frame sequence. The math is specified below in the equations section. However, in general, the geometric distortion between an observed distorted image and the reference image is obtained by the movement of control points placed in a grid on the images. B-spline interpolation is then used to map all the points from the reference image to the deformed image. A cost function is then calculated to minimize both the error between the mapped points and maximize symmetry between forward and backward deformations. After this cost function is minimized, a motion vector from the deformed image to the reference image is generated, and we can then register the deformed image and apply a sharpening filter to create a high-quality image.

Equations

The initial test points are taken with equal spacing and can be represented as: ${\hat {\mathbf {x} }}_{0i}=({\hat {x}}_{0i},{\hat {y}}_{0i})^{T}$

Running correlation algorithm, we found the deformed locations of the test points. This difference between the original positions and the deformed position are then stored in the deformation vector:

${\vec {\mathbf {p} }}=[{\hat {x}}_{0i}-{\hat {x}}_{i},{\hat {y}}_{0i}-{\hat {y}}_{i}]^{T}$ = $[\Delta {\hat {x}}_{1},...,\Delta {\hat {x}}_{i},\Delta {\hat {y}}_{1},...,\Delta {\hat {y}}_{i}]^{T}$

Using B-spline interpolation, we defined the spline basis $c_{i}$ and the basis function matrix $\mathbf {A{(x)}}$ for each pixel where $\mathbf {x} =(x,y)^{T}$.

By definition, $c_{i}$ is:

$c_{i}=\mathrm {B} ({\frac {x-{\hat {x}}_{0i}}{h_{x}}})$ $\mathrm {B} ({\frac {y-{\hat {y}}_{0i}}{h_{y}}})$, where $\mathrm {B} (z)={\begin{cases}2/3-(1-\left|z\right|/2)z^{2}{\mbox{, if }}0\leq \left|z\right|\leq 1\\(2-\left|z\right|)^{3}/6{\mbox{, if }}1\leq \left|z\right|\leq 2\\0{\mbox{, otherwise}}\\\end{cases}}$

The basis function matrix is then:

$\mathbf {A(x)} ={\begin{bmatrix}c_{1}&\cdots &c_{n}&0&\cdots &0\\0&\cdots &0&c_{1}&\cdots &c_{n}\end{bmatrix}}$

With the basis function matrix, we can then define the deformed coordinates for every pixel:

$\mathbf {W(x;} {\vec {\mathbf {p} }}\mathbf {)} =\mathbf {x} +\mathbf {A(x)} {\vec {\mathbf {p} }}$

As suggested in Professor Milanfar's paper, instead of using a classic B-spline registration approach to the estimation of the deformation vector through minimization of the following cost function,

$C({\overrightarrow {\mathbf {p} }})=\sum _{x}\left|R(\mathbf {W(x;} {\overrightarrow {\mathbf {p} }}\mathbf {)} )-G(\mathbf {x} )\right|^{2}$

which is prone to local minima traps, one can use the following cost function:

$C({\overrightarrow {\mathbf {p} }})=\sum _{x}\left|R(\mathbf {W(x;} {\overrightarrow {\mathbf {p} }}\mathbf {)} )-G(\mathbf {x} )\right|^{2}+\sum _{x}\left|G(\mathbf {W(x;} {\overleftarrow {\mathbf {p} }}\mathbf {)} )-R(\mathbf {x} )\right|^{2}+\gamma ({\overrightarrow {\mathbf {p} }}+{\overleftarrow {\mathbf {p} }})^{T}({\overrightarrow {\mathbf {p} }}+{\overleftarrow {\mathbf {p} }})$

by initial approximation that:

$\mathbf {p^{T}} =[{\overrightarrow {\mathbf {p} }}^{T},{\overleftarrow {\mathbf {p} }}^{T}]$

where we initially assume that

${\overleftarrow {\mathbf {p} }}=-{\overrightarrow {\mathbf {p} }}$

Then, after each iteration, the cost deformation vector can be updated as:

$\mathbf {p} ^{l+1}=\mathbf {p} ^{l}-\mathbf {H} ^{-1}\mathbf {b}$

The following equations describe the system:

$\mathbf {H} ={\begin{bmatrix}{\overrightarrow {\mathbf {H} }}+\mathbf {\gamma I} &\mathbf {\gamma I} \\\mathbf {\gamma I} &{\overleftarrow {\mathbf {H} }}+\mathbf {\gamma I} \end{bmatrix}}$, where $\gamma =5000$. ${\overrightarrow {\mathbf {H} }}=\sum _{x}{\overrightarrow {\mathbf {d} }}\mathbf {(x)} {\overrightarrow {\mathbf {d} }}\mathbf {(x)} ^{T},{\overleftarrow {\mathbf {H} }}=\sum _{x}{\overleftarrow {\mathbf {d} }}\mathbf {(x)} {\overleftarrow {\mathbf {d} }}\mathbf {(x)} ^{T}$

The faster implementation, as suggested in the paper is calculating ${\overrightarrow {\mathbf {H} }}$ and ${\overleftarrow {\mathbf {H} }}$ the following way:

${\overrightarrow {\mathbf {H} }}=\sum _{x}{\overrightarrow {\mathbf {H_{x}} }}$

where:

${\overrightarrow {\mathbf {H_{x}} }}={\begin{bmatrix}r_{x}(\mathbf {W} )^{2}\mathbf {C_{x}} &r_{x}(\mathbf {W} )r_{y}(\mathbf {W} )\mathbf {C_{x}} \\r_{x}(\mathbf {W} )r_{y}(\mathbf {W} )\mathbf {C_{x}} &r_{y}(\mathbf {W} )^{2}\mathbf {C_{x}} \end{bmatrix}}$ ${\overrightarrow {\mathbf {d} }}\mathbf {(x)^{T}} ={\frac {\partial R(\mathbf {W(x;} {\overrightarrow {\mathbf {p} }}^{l}\mathbf {)} )}{\partial \mathbf {W} }}\mathbf {A(x)} ,{\overleftarrow {\mathbf {d} }}\mathbf {(x)^{T}} ={\frac {\partial G(\mathbf {W(x;} {\overleftarrow {\mathbf {p} }}^{l}\mathbf {)} )}{\partial \mathbf {W} }}\mathbf {A(x)}$ ${\overrightarrow {\mathbf {b} }}=\sum _{x}{\overrightarrow {\mathbf {d} }}\mathbf {(x)} [R(\mathbf {W(x;} {\overrightarrow {\mathbf {p} }}^{l}))-G(\mathbf {x} )],{\overleftarrow {\mathbf {b} }}=\sum _{x}{\overleftarrow {\mathbf {d} }}\mathbf {(x)} [G(\mathbf {W(x;} {\overleftarrow {\mathbf {p} }}^{l}))-R(\mathbf {x} )]$

where ${\frac {\partial R(\mathbf {W(x;} {\overrightarrow {\mathbf {p} }}^{l}\mathbf {)} )}{\partial \mathbf {W} }}=[r_{x}(\mathbf {W} ),r_{y}(\mathbf {W} )]$ and is the gradient of the averaged image R at the deformed pixel location $\mathbf {W(x;} {\overrightarrow {\mathbf {p} }}^{l}\mathbf {)} )$.

MR Analysis

The MR data was analyzed using mrVista software tools.

Pre-processing

All data were slice-time corrected, motion corrected, and repeated scans were averaged together to create a single average scan for each subject. Et cetera.

PRF model fits

PRF models were fit with a 2-gaussian model.

MNI space

After a pRF model was solved for each subject, the model was trasnformed into MNI template space. This was done by first aligning the high resolution t1-weighted anatomical scan from each subject to an MNI template. Since the pRF model was coregistered to the t1-anatomical scan, the same alignment matrix could then be applied to the pRF model.
Once each pRF model was aligned to MNI space, 4 model parameters - x, y, sigma, and r^2 - were averaged across each of the 6 subjects in each voxel.

Et cetera.

Results - What you found

Retinotopic models in native space

Some text. Some analysis. Some figures.

Retinotopic models in individual subjects transformed into MNI space

Some text. Some analysis. Some figures.

Retinotopic models in group-averaged data on the MNI template brain

Some text. Some analysis. Some figures. Maybe some equations.

Equations

The initial test points are taken with equal spacing and can be represented as: ${\hat {\mathbf {x} }}_{0i}=({\hat {x}}_{0i},{\hat {y}}_{0i})^{T}$

Running correlation algorithm, we found the deformed locations of the test points. This difference between the original positions and the deformed position are then stored in the deformation vector:

${\vec {\mathbf {p} }}=[{\hat {x}}_{0i}-{\hat {x}}_{i},{\hat {y}}_{0i}-{\hat {y}}_{i}]^{T}$ = $[\Delta {\hat {x}}_{1},...,\Delta {\hat {x}}_{i},\Delta {\hat {y}}_{1},...,\Delta {\hat {y}}_{i}]^{T}$

Using B-spline interpolation, we defined the spline basis $c_{i}$ and the basis function matrix $\mathbf {A{(x)}}$ for each pixel where $\mathbf {x} =(x,y)^{T}$.

By definition, $c_{i}$ is:

$c_{i}=\mathrm {B} ({\frac {x-{\hat {x}}_{0i}}{h_{x}}})$ $\mathrm {B} ({\frac {y-{\hat {y}}_{0i}}{h_{y}}})$, where $\mathrm {B} (z)={\begin{cases}2/3-(1-\left|z\right|/2)z^{2}{\mbox{, if }}0\leq \left|z\right|\leq 1\\(2-\left|z\right|)^{3}/6{\mbox{, if }}1\leq \left|z\right|\leq 2\\0{\mbox{, otherwise}}\\\end{cases}}$

The basis function matrix is then:

$\mathbf {A(x)} ={\begin{bmatrix}c_{1}&\cdots &c_{n}&0&\cdots &0\\0&\cdots &0&c_{1}&\cdots &c_{n}\end{bmatrix}}$

With the basis function matrix, we can then define the deformed coordinates for every pixel:

$\mathbf {W(x;} {\vec {\mathbf {p} }}\mathbf {)} =\mathbf {x} +\mathbf {A(x)} {\vec {\mathbf {p} }}$

As suggested in Professor Milanfar's paper, instead of using a classic B-spline registration approach to the estimation of the deformation vector through minimization of the following cost function,

$C({\overrightarrow {\mathbf {p} }})=\sum _{x}\left|R(\mathbf {W(x;} {\overrightarrow {\mathbf {p} }}\mathbf {)} )-G(\mathbf {x} )\right|^{2}$

which is prone to local minima traps, one can use the following cost function:

$C({\overrightarrow {\mathbf {p} }})=\sum _{x}\left|R(\mathbf {W(x;} {\overrightarrow {\mathbf {p} }}\mathbf {)} )-G(\mathbf {x} )\right|^{2}+\sum _{x}\left|G(\mathbf {W(x;} {\overleftarrow {\mathbf {p} }}\mathbf {)} )-R(\mathbf {x} )\right|^{2}+\gamma ({\overrightarrow {\mathbf {p} }}+{\overleftarrow {\mathbf {p} }})^{T}({\overrightarrow {\mathbf {p} }}+{\overleftarrow {\mathbf {p} }})$

by initial approximation that:

$\mathbf {p^{T}} =[{\overrightarrow {\mathbf {p} }}^{T},{\overleftarrow {\mathbf {p} }}^{T}]$

where we initially assume that

${\overleftarrow {\mathbf {p} }}=-{\overrightarrow {\mathbf {p} }}$

Then, after each iteration, the cost deformation vector can be updated as:

$\mathbf {p} ^{l+1}=\mathbf {p} ^{l}-\mathbf {H} ^{-1}\mathbf {b}$

The following equations describe the system:

$\mathbf {H} ={\begin{bmatrix}{\overrightarrow {\mathbf {H} }}+\mathbf {\gamma I} &\mathbf {\gamma I} \\\mathbf {\gamma I} &{\overleftarrow {\mathbf {H} }}+\mathbf {\gamma I} \end{bmatrix}}$, where $\gamma =5000$. ${\overrightarrow {\mathbf {H} }}=\sum _{x}{\overrightarrow {\mathbf {d} }}\mathbf {(x)} {\overrightarrow {\mathbf {d} }}\mathbf {(x)} ^{T},{\overleftarrow {\mathbf {H} }}=\sum _{x}{\overleftarrow {\mathbf {d} }}\mathbf {(x)} {\overleftarrow {\mathbf {d} }}\mathbf {(x)} ^{T}$

The faster implementation, as suggested in the paper is calculating ${\overrightarrow {\mathbf {H} }}$ and ${\overleftarrow {\mathbf {H} }}$ the following way:

${\overrightarrow {\mathbf {H} }}=\sum _{x}{\overrightarrow {\mathbf {H_{x}} }}$

where:

${\overrightarrow {\mathbf {H_{x}} }}={\begin{bmatrix}r_{x}(\mathbf {W} )^{2}\mathbf {C_{x}} &r_{x}(\mathbf {W} )r_{y}(\mathbf {W} )\mathbf {C_{x}} \\r_{x}(\mathbf {W} )r_{y}(\mathbf {W} )\mathbf {C_{x}} &r_{y}(\mathbf {W} )^{2}\mathbf {C_{x}} \end{bmatrix}}$ ${\overrightarrow {\mathbf {d} }}\mathbf {(x)^{T}} ={\frac {\partial R(\mathbf {W(x;} {\overrightarrow {\mathbf {p} }}^{l}\mathbf {)} )}{\partial \mathbf {W} }}\mathbf {A(x)} ,{\overleftarrow {\mathbf {d} }}\mathbf {(x)^{T}} ={\frac {\partial G(\mathbf {W(x;} {\overleftarrow {\mathbf {p} }}^{l}\mathbf {)} )}{\partial \mathbf {W} }}\mathbf {A(x)}$ ${\overrightarrow {\mathbf {b} }}=\sum _{x}{\overrightarrow {\mathbf {d} }}\mathbf {(x)} [R(\mathbf {W(x;} {\overrightarrow {\mathbf {p} }}^{l}))-G(\mathbf {x} )],{\overleftarrow {\mathbf {b} }}=\sum _{x}{\overleftarrow {\mathbf {d} }}\mathbf {(x)} [G(\mathbf {W(x;} {\overleftarrow {\mathbf {p} }}^{l}))-R(\mathbf {x} )]$

where ${\frac {\partial R(\mathbf {W(x;} {\overrightarrow {\mathbf {p} }}^{l}\mathbf {)} )}{\partial \mathbf {W} }}=[r_{x}(\mathbf {W} ),r_{y}(\mathbf {W} )]$ and is the gradient of the averaged image R at the deformed pixel location $\mathbf {W(x;} {\overrightarrow {\mathbf {p} }}^{l}\mathbf {)} )$.

Retinotopic models in group-averaged data projected back into native space

Some text. Some analysis. Some figures.

Conclusions

Here is where you say what your results mean.

References - Resources and related work

References

Software

Appendix I - Code and Data

Code

File:CodeFile.zip

Data

zip file with my data

Appendix II - Work partition (if a group project)

Brian and Bob gave the lectures. Jon mucked around on the wiki.