Fast Image Deconvolution using Hyper-Laplacian Priors

function [w] = solve_image(v, beta, alpha) % % solve the following component-wise separable problem % min |w|^\alpha + \frac{\beta}{2} (w - v).^2 % % A LUT is used to solve the problem; when the function is first called % for a new value of beta or alpha, a LUT is built for that beta/alpha % combination and for a range of values of v. The LUT stays persistent % between calls to solve_image. It will be recomputed the first time this % function is called. % range of input data and step size; increasing the range of decreasing % the step size will increase accuracy but also increase the size of the % LUT range = 10; step = 0.0001; persistent lookup_v known_beta xx known_alpha ind = find(known_beta==beta & known_alpha==alpha); if isempty(known_beta | known_alpha) xx = [-range:step:range]; end if any(ind) fprintf('Reusing lookup table for beta %.3g and alpha %.3g\n', beta, alpha); %%% already computed if (exist('pointOp') == 3) % Use Eero Simoncelli's function to extrapolate w = pointOp(double(v),lookup_v(ind,:), -range, step, 0); else w = interp1(xx', lookup_v(ind,:)', v(:), 'linear', 'extrap'); w = reshape(w, size(v,1), size(v,2)); end; else %%% now go and recompute xx for new value of beta and alpha tmp = compute_w(xx, beta, alpha); lookup_v = [lookup_v; tmp(:)']; known_beta = [known_beta, beta]; known_alpha = [known_alpha, alpha]; %%% and lookup current v's in the new lookup table row. if (exist('pointOp') == 3) % Use Eero Simoncelli's function to extrapolate w = pointOp(double(v),lookup_v(end,:), -range, step, 0); else w = interp1(xx', lookup_v(end,:)', v(:), 'linear', 'extrap'); w = reshape(w, size(v,1), size(v,2)); end; fprintf('Recomputing lookup table for new value of beta %.3g and alpha %.3g\n', beta, alpha); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % call different functions to solve the minimization problem % min |w|^\alpha + \frac{\beta}{2} (w - v).^2 for a fixed beta and alpha % function w = compute_w(v, beta, alpha) if (abs(alpha - 1) < 1e-9) % assume alpha = 1.0 w = compute_w1(v, beta); return; end; if (abs(alpha - 2/3) < 1e-9) % assume alpha = 2/3 w = compute_w23(v, beta); return; end; if (abs(alpha - 1/2) < 1e-9) % assume alpha = 1/2 w = compute_w12(v, beta); return; end; % for any other value of alpha, plug in some other generic root-finder % here, we use Newton-Raphson w = newton_w(v, beta, alpha); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function w = compute_w23(v, beta) % solve a quartic equation % for alpha = 2/3 epsilon = 1e-6; %% tolerance on imag part of real root k = 8/(27*beta^3); m = ones(size(v))*k; % Now use formula from % http://en.wikipedia.org/wiki/Quartic_equation (Ferrari's method) % running our coefficients through Mathmetica (quartic_solution.nb) % optimized to use as few operations as possible... %%% precompute certain terms v2 = v .* v; v3 = v2 .* v; v4 = v3 .* v; m2 = m .* m; m3 = m2 .* m; %% Compute alpha & beta alpha = -1.125*v2; beta2 = 0.25*v3; %%% Compute p,q,r and u directly. q = -0.125*(m.*v2); r1 = -q/2 + sqrt(-m3/27 + (m2.*v4)/256); u = exp(log(r1)/3); y = 2*(-5/18*alpha + u + (m./(3*u))); W = sqrt(alpha./3 + y); %%% now form all 4 roots root = zeros(size(v,1),size(v,2),4); root(:,:,1) = 0.75.*v + 0.5.*(W + sqrt(-(alpha + y + beta2./W ))); root(:,:,2) = 0.75.*v + 0.5.*(W - sqrt(-(alpha + y + beta2./W ))); root(:,:,3) = 0.75.*v + 0.5.*(-W + sqrt(-(alpha + y - beta2./W ))); root(:,:,4) = 0.75.*v + 0.5.*(-W - sqrt(-(alpha + y - beta2./W ))); %%%%%% Now pick the correct root, including zero option. %%% Clever fast approach that avoids lookups v2 = repmat(v,[1 1 4]); sv2 = sign(v2); rsv2 = real(root).*sv2; %%% condensed fast version %%% take out imaginary roots above v/2 but below v root_flag3 = sort(((abs(imag(root))<epsilon) & ((rsv2)>(abs(v2)/2)) & ((rsv2)<(abs(v2)))).*rsv2,3,'descend').*sv2; %%% take best w=root_flag3(:,:,1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function w = compute_w12(v, beta) % solve a cubic equation % for alpha = 1/2 epsilon = 1e-6; %% tolerance on imag part of real root k = -0.25/beta^2; m = ones(size(v))*k.*sign(v); %%%%%%%%%%%%%%%%%%%%%%%%%%% Compute the roots (all 3) t1 = (2/3)*v; v2 = v .* v; v3 = v2 .* v; %%% slow (50% of time), not clear how to speed up... t2 = exp(log(-27*m - 2*v3 + (3*sqrt(3))*sqrt(27*m.^2 + 4*m.*v3))/3); t3 = v2./t2; %%% find all 3 roots root = zeros(size(v,1),size(v,2),3); root(:,:,1) = t1 + (2^(1/3))/3*t3 + (t2/(3*2^(1/3))); root(:,:,2) = t1 - ((1+i*sqrt(3))/(3*2^(2/3)))*t3 - ((1-i*sqrt(3))/(6*2^(1/3)))*t2; root(:,:,3) = t1 - ((1-i*sqrt(3))/(3*2^(2/3)))*t3 - ((1+i*sqrt(3))/(6*2^(1/3)))*t2; root(find(isnan(root) | isinf(root))) = 0; %%% catch 0/0 case %%%%%%%%%%%%%%%%%%%%%%%%%%%% Pick the right root %%% Clever fast approach that avoids lookups v2 = repmat(v,[1 1 3]); sv2 = sign(v2); rsv2 = real(root).*sv2; root_flag3 = sort(((abs(imag(root))<epsilon) & ((rsv2)>(2*abs(v2)/3)) & ((rsv2)<(abs(v2)))).*rsv2,3,'descend').*sv2; %%% take best w=root_flag3(:,:,1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function w = compute_w1(v, beta) % solve a simple max problem for alpha = 1 w = max(abs(v) - 1/beta, 0).*sign(v); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function w = newton_w(v, beta, alpha) % for a general alpha, use Newton-Raphson; more accurate root-finders may % be substituted here; we are finding the roots of the equation: % \alpha*|w|^{\alpha - 1} + \beta*(v - w) = 0 iterations = 4; x = v; for a=1:iterations fd = (alpha)*sign(x).*abs(x).^(alpha-1)+beta*(x-v); fdd = alpha*(alpha-1)*abs(x).^(alpha-2)+beta; x = x - fd./fdd; end; q = find((x)); x(q) = 0;isnan % check whether the zero solution is the better one z = beta/2*v.^2; f = abs(x).^alpha + beta/2*(x-v).^2; w = (f<z).*x;