function y = multSID2D(sid, w, Nv, Nh)
% multSID2D       Multiplication for 2D Shift Invariant Dictionary: y = D*w
%
% use:   y = multSID2D(sid, w, size(y));
%        y = multSID2D(sid, w, Nv, Nh);
% -------------------------------------------------------------------------
%   y       a Nv x Nh image (matrix)
%   sid     The 2D Shift Invariant Dictionary, a struct as in initSID2D
%   w       a (usually large and sparse) coefficient vector Kx1 
%   Nv, Nh  size of image y, given as two integers or an array with two elements
% -------------------------------------------------------------------------
% ex:  clear all;
%      sid = initSID2D([2,3], [2,4], [2,1], [2,1], [ones(4,1);(13:24)']); % atoms: 2x2 and 3x4 .
%      Nv = 6; Nh = 8; N = Nv*Nh;   % image is 6x8, 3x4 2x2 blocks, 4x5 3x4 blocks.
%      w = zeros(3*4+4*5, 1); w(1)=4; w(8)=2; w(end) = 1;
%      y = multSID2D(sid, w, Nv, Nh);
%      wmp = saSID2Domp(y, sid, 3, 'mexOMP', [], 1);

%----------------------------------------------------------------------
% Made by Karl Skretting, University of Stavanger (UiS), 
% Department of Electrical Engineering and Computer Science (IDE)
% Mail:  karl.skretting@uis.no   Homepage:  http://www.ux.uis.no/~karlsk/
% 
% HISTORY:  dd.mm.yyyy
% Ver. 1.0  12.08.2016  KS: made file 
%----------------------------------------------------------------------

if ~isstruct(sid) || ~isfield(sid,'Type') || ~strcmp(sid.Type, 'SID2D')
    disp([mfilename,': First argument not a 2D SID structure.']);
    y = [];
    return
end
if nargin < 4
    if (nargin == 3) && (numel(Nv) == 2)
        Nh = Nv(2);
        Nv = Nv(1);
    else
        disp([mfilename,': image size must be given, see help.']);
        y = [];
        return
    end
end

im = thisKfun(sid, Nv, Nh);    % how coefficients in w are for image y

if numel(w) ~= im.K 
    disp([mfilename,': size of image (Nv and Nh) and numel(w) do not match.']);
    y = [];
    return
end

% Assume that size of expanded dictionary (Nv*Nh)x(numel(w)) is too large
% to make it possible to build the huge dictionary D and find y as
%    y = reshape(D*w, Nv, Nh);
% saSID2Domp has a subfunction, D = thisDblk(sid, blk), 
% that makes D for a small block.

coeff = nonzeros(w);
idx = find(w);
y = zeros(Nv,Nh);
m = 1;
atom = reshape(sid.freeVar(sid.i1(m):sid.i2(m)),sid.Qv(m),sid.Qh(m));

for i = 1:numel(idx)
    k = idx(i);
    c = coeff(i);   % coefficient
    while (k > im.k2(m)); 
        m = m + 1; 
        atom = reshape(sid.freeVar(sid.i1(m):sid.i2(m)),sid.Qv(m),sid.Qh(m));
    end;
    k = (k-im.k1(m));          % in  0,1,2,...,im.Km(m)-1
    kh = floor( k/im.Kv(m) );  % in  0,1,2,...,im.Kh(m)-1
    kv = k - kh*im.Kv(m);      % in  0,1,2,...,im.Kv(m)-1
    nv = sid.Pv(m)*kv;
    nh = sid.Ph(m)*kh;
    v = nv + (1:sid.Qv(m));   % vertical indexes
    h = nh + (1:sid.Qh(m));   % horizontal indexes
    y(v,h) = y(v,h) + c*atom;
end

end


%%

% Find how coefficients are organized. 
% This subfunction is also in saSID2Domp.m.
function out = thisKfun(sid, Nv, Nh)
Kv = ceil( (Nv+1-sid.Qv)./sid.Pv );
Kh = ceil( (Nh+1-sid.Qh)./sid.Ph );
Km = Kv.*Kh;
K = sum(Km);
k2 = cumsum(Km);
k1 = [1, reshape(k2(1:(sid.M-1)),1,sid.M-1)+1];
out = struct('N',Nv*Nh, 'K',K, 'M',sid.M, ...
             'Kv',Kv,     'Kh',Kh,    'Km',Km, ...
             'k1',k1,     'k2',k2 );             
%             'Qv',sid.Qv, 'Qh',sid.Qh, ...
%             'Pv',sid.Pv, 'Ph',sid.Ph, ...
%             'Nv',Nv,     'Nh',Nh, ...
end
