Deblur: Richardson-Lucy deconvolution

If we know how an image was blurred, we can deblur it.

• The Script
• Stars
• References
• Scripts

A point spread function (PSF) defines how a point of light is spread across an image. If we know the PSF that formed an image, we can derive the image as it was (more or less) before the light was spread by the PSF. For example, if a photograph suffers from camera-shake, and we know the exact nature of the motion, we can restore the photograph as if the camera hadn't moved.

The technique shown here is the Richardson-Lucy deconvolution.

For some aspects of photographic blur, including the concept of PSFs, see my Camera blurs.

The Script

The script needs the PSF (or "convolution kernel") in the form of an image. The script then converts that to text, because ImageMagick currently requires kernels to be specified as text. This conversion, from image to text, is performed by my process module img2knl (see blah).

My script closely follows the algorithm given in Wikipedia. See Python documentation to understand how ...

psf_hat = psf[::-1,::-1]

... translates to IM:

convert psf.png -flip -flop psfhat.png

The script needs IM to be compiled with HDRI (I think: check! blah). It does not need FFTW.

A sample image cb_samp.png ...
... with this kernel cb_poly5.png ...
... is blurred in linear colorspace making this result drl_conv.png: for /F "usebackq" %%L in (`%IMDEV%convert ^ cb_poly5.png -process img2knl NULL:`) do set KNL=%%L %IM%convert ^ cb_samp.png ^ -colorspace RGB ^ -define convolve:scale="^!" ^ -morphology convolve %KNL% ^ -colorspace sRGB ^ drl_conv.png
We do 3 iterations of a Richardson-Lucy deconvolution of the blurred image with the kernel: call %PICTBAT%deconvRL drl_conv.png cb_poly5.png 3 drl_d3.png
We do 10 iterations of a Richardson-Lucy deconvolution of the blurred image with the kernel: rem call %PICTBAT%deconvRL drl_conv.png cb_poly5.png 10 drl_d10.png
We do 100 iterations of a Richardson-Lucy deconvolution of the blurred image with the kernel: rem call %PICTBAT%deconvRL drl_conv.png cb_poly5.png 100 drl_d100.png
We do 1000 iterations of a Richardson-Lucy deconvolution of the blurred image with the kernel: rem call %PICTBAT%deconvRL drl_conv.png cb_poly5.png 1000 drl_d1000.png

Stars

I took a photograph of the night sky with a Nikon D800 camera, Noct-Nikkor 58mm f/1.2 lens, ISO 6400, at f/1.2, 1/10 second, hand-held. The sky was not entirely dark.

Show the image resized, and a crop. set STAR_SRC=\pictures\20120922\DSC_0389_gms.tiff set STAR_CROP=-gravity Center -crop 600x400-300-200 +repage %IM%convert ^ %STAR_SRC% ^ ( +clone -resize 600x600 -auto-level -write drl_stars1.png +delete ) ^ %STAR_CROP% ^ -auto-level ^ drl_stars_crop1.png

Pick a bright star (or maybe it is a planet), and process it to make a kernel.

%IM%convert ^ %STAR_SRC% ^ -crop 19x14+2976+1828 +repage ^ -colorspace gray ^ -auto-level ^ -write drl_star_raw.png ^ ( +clone ^ -scale 1000%% ^ -write drl_star_raw_sc.png ^ +delete ) ^ -level 70%%,100%% ^ -trim +repage ^ -write drl_star_knl.png ^ -scale 1000%% ^ drl_star_knl_sc.png rem -sigmoidal-contrast 20,80%% ^ rem -level 10%%,100%% ^ %IM%convert ^ -size 13x13 xc:#000 ^ -fill #fff -draw "line 0,0 6,6" ^ -write drl_star_knlX.png ^ -scale 1000%% ^ drl_star_knl_scX.png

Deconvolve the photo with this kernel.

call %PICTBAT%deconvRL drl_stars_crop1.png drl_star_knl.png 10 drl_stars_c1_d.png
call %PICTBAT%deconvRL %STAR_SRC% drl_star_knl.png 3 drl_star_d.png %IMDEV%convert ^ drl_star_d.png ^ -auto-level ^ ( +clone -resize 600x600 -auto-level -write drl_stars_dc1.png +delete ) ^ %STAR_CROP% ^ -auto-level ^ drl_stars_dc1_sc.png

References

• Bayesian-Based Iterative Method of Image Restoration, Richardson, 1972.
• An iterative technique for the rectification of observed distributions, Lucy, 1974
• Wikipedia: Richardson-Lucy deconvolution, retrieved 10-November-2014.

Scripts

deconvRL.bat

rem Richardson-Lucy deconvolution.
rem %1 is image.
rem %2 is kernel (PSF).
rem %3 is number of iterations (default 10).
rem %4 is optional name of output file.
@rem
@rem Ref http://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution
@rem
@rem We need option to NOT convert to linear.

@if "%2"=="" findstr /B "rem @rem" %~f0 & exit /B 1

@setlocal

rem @call echoOffSave

call %PICTBAT%setInOut %1 drl

set NITER=%3
if "%NITER%"=="." set NITER=
if "%NITER%"=="" set NITER=10

if not "%4"=="" set OUTFILE=%4

set TMPDIR=%CYGTEMP%/
set TMPEXT=.miff

set PSF=%2
set PSFHAT=%TMPDIR%drl_psfhat%TMPEXT%
set EST=%TMPDIR%drl_est%TMPEXT%
set TMPIN=%TMPDIR%drl_in%TMPEXT%

%IMDEV%convert ^
  %PSF% ^
  -flip -flop ^
  %PSFHAT%

%IMDEV%convert ^
  %INFILE% ^
  -colorspace RGB ^
  +depth ^
  -write %TMPIN% ^
  %EST%

rem  -fill gray50 -colorize 100


rem for /F "usebackq" %%L in (`%IMDEV%convert ^
rem   %PSF% -process img2knl NULL:`) do set PSF_KNL=%%L

rem for /F "usebackq" %%L in (`%IMDEV%convert ^
rem   %PSFHAT% -process img2knl NULL:`) do set PSFHAT_KNL=%%L

%IMDEV%convert ^
  %PSF% -process img2knl NULL: >psf_knl.txt

%IMDEV%convert ^
  %PSFHAT% -process img2knl NULL: >psfhat_knl.txt

for /L %%i in (1,1,%NITER%) do %IMDEV%convert ^
  %EST% ^
  ( +clone ^
    ( +clone ^
      -define convolve:scale="^!" ^
      -morphology convolve @psf_knl.txt ^
      %TMPIN% ^
      -compose DivideDst -composite ^
      -morphology convolve @psfhat_knl.txt ^
    ) ^
    -compose Multiply -composite ^
  ) ^
  ( -clone 0-1 ^
    -metric RMSE ^
    -format "%%[distortion]\n" -compare -write info: ^
    -delete 0-1 ^
  ) ^
  -delete 0 ^
  +depth ^
  %EST%

%IMDEV%convert ^
  %EST% ^
  -set colorspace RGB ^
  -colorspace sRGB ^
  +depth ^
  %OUTFILE%

%IM%identify -version

Version: ImageMagick 6.9.0-0 Q16 x64 2014-11-14 http://www.imagemagick.org
Copyright: Copyright (C) 1999-2014 ImageMagick Studio LLC
Features: DPC OpenMP
DelegatesIM%identify -version

Version: ImageMagick 6.9.0-0 Q16 x64 2014-11-14 http://www.imagemagick.org
Copyright: Copyright (C) 1999-2014 ImageMagick Studio LLC
Features: DPC OpenMP
Delegates (built-in): bzlib cairo freetype jbig jng jp2 jpeg lcms lqr pangocairo png ps rsvg tiff webp xml zlib

Page version v1.0 10-November-2014.

Page created 05-Mar-2015 09:26:37.

Copyright © 2015 Alan Gibson.