annotate toolboxes/AudioInpaintingToolbox/Solvers/inpaintFrame_consOMP_Gabor.m @ 244:5c8bcdadb380 unlocbox

added UNLocBox interface, example and a 128x128 Lena
author bmailhe
date Tue, 04 Sep 2012 11:00:36 +0100
parents 56d719a5fd31
children
rev   line source
ivan@138 1 function y = inpaintFrame_consOMP_Gabor(problemData,param)
ivan@138 2 % Inpainting method based on OMP with a constraint
ivan@138 3 % on the amplitude of the reconstructed samples an optional constraint
ivan@138 4 % on the maximum value of the clipped samples, and using the Gabor dictionary
ivan@138 5 % generated by Gabor_Dictionary.m. The method jointly selects
ivan@138 6 % cosine and sine atoms at the same frequency.
ivan@138 7 %
ivan@138 8 % Usage: y = inpaintFrame_consOMP_Gabor(problemData,param)
ivan@138 9 %
ivan@138 10 %
ivan@138 11 % Inputs:
ivan@138 12 % - problemData.x: observed signal to be inpainted
ivan@138 13 % - problemData.Imiss: Indices of clean samples
ivan@138 14 % - param.D - the dictionary matrix (optional if param.D_fun is set)
ivan@138 15 % - param.D_fun - a function handle that generates the dictionary
ivan@138 16 % matrix param.D if param.D is not given. See Gabor_Dictionary.m
ivan@138 17 % - param.wa - Analysis window
ivan@138 18 % - param.Upper_Limit - if present and non-empty this fiels
ivan@138 19 % indicates that an upper limit constraint is active and its
ivan@138 20 % integer value is such that
ivan@138 21 %
ivan@138 22 % Outputs:
ivan@138 23 % - y: estimated frame
ivan@138 24 %
ivan@138 25 % Note that the CVX library is needed.
ivan@138 26 %
ivan@138 27 % -------------------
ivan@138 28 %
ivan@138 29 % Audio Inpainting toolbox
ivan@138 30 % Date: June 28, 2011
ivan@138 31 % By Valentin Emiya, Amir Adler, Michael Elad, Maria Jafari
ivan@138 32 % This code is distributed under the terms of the GNU Public License version 3 (http://www.gnu.org/licenses/gpl.txt).
ivan@138 33
ivan@138 34
ivan@138 35
ivan@138 36
ivan@138 37 x = problemData.x;
ivan@138 38 IObs = find(~problemData.IMiss);
ivan@138 39 p.N = length(x);
ivan@138 40 E2 = param.OMPerr^2;
ivan@138 41 E2M=E2*length(IObs);
ivan@138 42 wa = param.wa(param.N);
ivan@138 43
ivan@138 44 % build the dictionary matrix if only the dictionary generation function is given
ivan@138 45 if ~isfield(param,'D')
ivan@138 46 param.D = param.D_fun(param);
ivan@138 47 end
ivan@138 48
ivan@138 49
ivan@138 50 % clipping level detection
ivan@138 51 clippingLevelEst = max(abs(x(:)./wa(:)));
ivan@138 52 IMiss = true(length(x),1);
ivan@138 53 IMiss(IObs) = false;
ivan@138 54 IMissPos = find(x>=0 & IMiss);
ivan@138 55 IMissNeg = find(x<0 & IMiss);
ivan@138 56
ivan@138 57 DictPos=param.D(IMissPos,:);
ivan@138 58 DictNeg=param.D(IMissNeg,:);
ivan@138 59
ivan@138 60 % Clipping level: take the analysis window into account
ivan@138 61 wa_pos = wa(IMissPos);
ivan@138 62 wa_neg = wa(IMissNeg);
ivan@138 63 b_ineq_pos = wa_pos(:)*clippingLevelEst;
ivan@138 64 b_ineq_neg = -wa_neg(:)*clippingLevelEst;
ivan@138 65 if isfield(param,'Upper_Limit') && ~isempty(param.Upper_Limit)
ivan@138 66 b_ineq_pos_upper_limit = wa_pos(:)*param.Upper_Limit*clippingLevelEst;
ivan@138 67 b_ineq_neg_upper_limit = -wa_neg(:)*param.Upper_Limit*clippingLevelEst;
ivan@138 68 else
ivan@138 69 b_ineq_pos_upper_limit = Inf;
ivan@138 70 b_ineq_neg_upper_limit = -Inf;
ivan@138 71 end
ivan@138 72
ivan@138 73 %%
ivan@138 74 Dict=param.D(IObs,:);
ivan@138 75 W=1./sqrt(diag(Dict'*Dict));
ivan@138 76 Dict=Dict*diag(W);
ivan@138 77 Dict1 = Dict(:,1:end/2);
ivan@138 78 Dict2 = Dict(:,end/2+1:end);
ivan@138 79 Dict1Dict2 = sum(Dict1.*Dict2);
ivan@138 80 n12 = 1./(1-Dict1Dict2.^2);
ivan@138 81 xObs=x(IObs);
ivan@138 82 %K = size(param.D,2);
ivan@138 83
ivan@138 84 residual=xObs;
ivan@138 85 maxNumCoef = param.sparsityDegree;
ivan@138 86 indx = [];
ivan@138 87 % currResNorm2 = sum(residual.^2);
ivan@138 88 currResNorm2 = E2M*2; % set a value above the threshold in order to have/force at least one loop executed
ivan@138 89 j = 0;
ivan@138 90 while currResNorm2>E2M && j < maxNumCoef,
ivan@138 91 j = j+1;
ivan@138 92 proj=residual'*Dict;
ivan@138 93 proj1 = proj(1:end/2);
ivan@138 94 proj2 = proj(end/2+1:end);
ivan@138 95
ivan@138 96 alpha_j = (proj1-Dict1Dict2.*proj2).*n12;
ivan@138 97 beta_j = (proj2-Dict1Dict2.*proj1).*n12;
ivan@138 98
ivan@138 99 err_j = sum(abs(repmat(residual,1,size(Dict1,2))-Dict1*sparse(diag(alpha_j))-Dict2*sparse(diag(beta_j))).^2);
ivan@138 100 [dum pos] = min(err_j);
ivan@138 101
ivan@138 102 indx(end+1)=pos;
ivan@138 103 indx(end+1)=pos+size(Dict1,2);
ivan@138 104 a=pinv(Dict(:,indx(1:2*j)))*xObs;
ivan@138 105 residual=xObs-Dict(:,indx(1:2*j))*a;
ivan@138 106 currResNorm2=sum(residual.^2);
ivan@138 107 end;
ivan@138 108
ivan@138 109 %% Constrained reestimation of the non-zero coefficients
ivan@138 110 j = length(indx);
ivan@138 111 if isinf(b_ineq_pos_upper_limit)
ivan@138 112 %% CVX code
ivan@138 113 cvx_begin
ivan@138 114 cvx_quiet(true)
ivan@138 115 variable a(j)
ivan@138 116 minimize(norm(Dict(:,indx)*a-xObs))
ivan@138 117 subject to
ivan@138 118 DictPos(:,indx)*(W(indx).*a) >= b_ineq_pos
ivan@138 119 DictNeg(:,indx)*(W(indx).*a) <= b_ineq_neg
ivan@138 120 cvx_end
ivan@138 121 if cvx_optval>1e3
ivan@138 122 cvx_begin
ivan@138 123 cvx_quiet(true)
ivan@138 124 variable a(j)
ivan@138 125 minimize(norm(Dict(:,indx)*a-xObs))
ivan@138 126 cvx_end
ivan@138 127 end
ivan@138 128 else
ivan@138 129 %% CVX code
ivan@138 130 cvx_begin
ivan@138 131 cvx_quiet(true)
ivan@138 132 variable a(j)
ivan@138 133 minimize(norm(Dict(:,indx)*a-xObs))
ivan@138 134 subject to
ivan@138 135 DictPos(:,indx)*(W(indx).*a) >= b_ineq_pos
ivan@138 136 DictNeg(:,indx)*(W(indx).*a) <= b_ineq_neg
ivan@138 137 DictPos(:,indx)*(W(indx).*a) <= b_ineq_pos_upper_limit
ivan@138 138 DictNeg(:,indx)*(W(indx).*a) >= b_ineq_neg_upper_limit
ivan@138 139 cvx_end
ivan@138 140 if cvx_optval>1e3
ivan@138 141 cvx_begin
ivan@138 142 cvx_quiet(true)
ivan@138 143 variable a(j)
ivan@138 144 minimize(norm(Dict(:,indx)*a-xObs))
ivan@138 145 cvx_end
ivan@138 146 end
ivan@138 147 end
ivan@138 148
ivan@138 149 %% Frame Reconstruction
ivan@138 150 indx(length(a)+1:end) = [];
ivan@138 151
ivan@138 152 Coeff = sparse(size(param.D,2),1);
ivan@138 153 if (~isempty(indx))
ivan@138 154 Coeff(indx) = a;
ivan@138 155 Coeff = W.*Coeff;
ivan@138 156 end
ivan@138 157 y = param.D*Coeff;
ivan@138 158
ivan@138 159 return
ivan@138 160