function D = makeSID2Dmatrix(sid, Nv, Nh, verbose)
% makeSID2Dmatrix Make sparse matrix D from 2D Shift Invariant Dictionary 
%
% use:   D = makeSID2Dmatrix(sid, Nv Nh, verbose);
% -------------------------------------------------------------------------
%   D     A dictionary matrix of size NxK  (K is calculated from sid)
%   sid   The 2D Shift Invariant Dictionary, a struct as in initSID2D
%   Nv,Nh Size for signal (or signal block) y. 
% -------------------------------------------------------------------------

%--------------------------------------------------------------------------
% 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  28.06.2016  KS: made file 
%--------------------------------------------------------------------------

if ~isstruct(sid) || ~isfield(sid,'Type') || ~strcmp(sid.Type, 'SID2D')
    disp([mfilename,': First argument not a 2D SI-FSD structure.']);
    D = [];
    return;
end

if (nargin < 4); verbose = 0; end;

M = sid.M;
Q = sid.Q;
N = Nv*Nh;

if verbose >= 1
    fprintf('\n%s: 2D SI-FSD has %i atoms kernels and %i free variables.\n', mfilename, M, Q);
    fprintf('In dictionary matrix D atom length is %ix%i = %i.\n', Nv, Nh, N);
    fprintf(['2D SI-FSD atom length: Qv = [%3i',char(kron(ones(1,M-1),',%3i')),']\n'], sid.Qv);
    fprintf(['2D SI-FSD atom length: Qh = [%3i',char(kron(ones(1,M-1),',%3i')),']\n'], sid.Qh);
    fprintf(['2D SI-FSD atom shift:  Pv = [%3i',char(kron(ones(1,M-1),',%3i')),']\n'], sid.Pv);
    fprintf(['2D SI-FSD atom shift:  Ph = [%3i',char(kron(ones(1,M-1),',%3i')),']\n'], sid.Ph);
end

km1 = zeros(1,M);  % index for first index in D for block Dm
km2 = zeros(1,M);  % index for last index in D for block Dm
supportDm = cell(M,1);
nnzD = 0;
for m=1:M
    Qv = sid.Qv(m);
    Qh = sid.Qh(m);
    Qm = Qv*Qh;
    Pv = sid.Pv(m);
    Ph = sid.Ph(m);
    Kv = ceil((Nv+1-Qv)/Pv);
    Kh = ceil((Nh+1-Qh)/Ph);
    Km = Kv*Kh;
    if (m==1); km1(1)=1; else km1(m)=km2(m-1)+1; end;
    km2(m) = km1(m) + Km - 1;
    fprintf('m=%2i: Qm=Qv*Qh = %4i=%2i*%2i, Pv=%2i, Ph=%2i, --> Kv=%3i, Kh=%3i, Km=%5i.  k1:k2 %5i:%5i\n',...
        m, Qm, Qv, Qh, Pv, Ph, Kv, Kh, Km, km1(m), km2(m) );
    supportDm{m} = reshape(repmat((1:Qv)',1,Qh) + repmat((0:(Qh-1))*Nv,Qv,1),Qm,1);
    nnzD = nnzD + Qm*Km;
end
K = km2(M);
fprintf('total size of D is NxK:  %i x %i,  nnz(D) = %6i.\n', N,K,nnzD);
D = sparse([],[],[],N,K,nnzD);

end

% 
% 
% 
% Km = floor( (N+sid.Pm-sid.Qm)./sid.Pm );
% K = sum(Km);
% k2 = cumsum(Km);
% k1 = [1; k2(1:(sid.M-1))+1];
% nnZeros = sum((k2-k1+1).*sid.Qm);
% % disp(['Number of non-zeros in D (',int2str(N),'x',int2str(K),') is ',int2str(nnZeros),'.']);
% 
% D = sparse([],[],[],N,K,nnZeros);
% for m = 1:sid.M
%     atom = sid.freeVar(sid.i1(m):sid.i2(m));
%     n = 1:sid.Qm(m);
%     for k = k1(m):k2(m)
%         D(n,k) = atom; %#ok<SPRIX>
%         n = n+sid.Pm(m);
%     end
% end
% 
% end
%     
