function [A, b] = VM_3lifting_3(Np, Nd, x1, y1, x2, y2, x3, y3, phik, xs, Vr)
% compute the approximated linear system for the desired numbers of
% vanishing moments for the optimization problem L3n
% Np primal and Nd dual vanishing moments
% xi, yi -- size of the ith lifting filter
% Copyright (c) 2006 Yi Chen

[Ny, Nx] = meshgrid(1:y1, 1:x1);
n10 = reshape(Nx+Ny-x1/2-1, [x1*y1,1]);
n11 = reshape(Nx-Ny-x1/2, [x1*y1,1]);

[Ny, Nx] = meshgrid(1:y2, 1:x2);
n20 = reshape(Nx+Ny-x2/2-1, [x2*y2,1]);
n21 = reshape(Nx-Ny-x2/2, [x2*y2,1]);

[Ny, Nx] = meshgrid(1:y3, 1:x3);
n30 = reshape(Nx+Ny-x3/2-1, [x3*y3,1]);
n31 = reshape(Nx-Ny-x3/2, [x3*y3,1]);

i1 = ones(x1*y1,1);
i2 = ones(x2*y2,1);
i3 = ones(x3*y3,1);

I1 = [eye(x1*y1) zeros(x1*y1, x2*y2+x3*y3)];
I2 = [zeros(x2*y2, x1*y1) eye(x2*y2) zeros(x2*y2, x3*y3)];
I3 = [zeros(x3*y3, x1*y1+x2*y2) eye(x3*y3)];

a1_hat = I1*(xs+Vr*phik);
a2_hat = I2*(xs+Vr*phik);
a3_hat = I3*(xs+Vr*phik);

V1_hat = I1*Vr;
V2_hat = I2*Vr;
V3_hat = I3*Vr;

if Nd == 4
bd(1) = a1_hat'*(n10.*n11.*2.*i1) + a3_hat'*(n30.*n31.*2.*i3)...
    + a2_hat'*((n20.*n21.*2.*i2)*2*i3'+2*i2*(n30.*n31.*2.*i3)')*a3_hat...
    + a1_hat'*(n10.*n11.*2.*i1)*a2_hat'*2*i2*a3_hat'*2*i3...
    + a1_hat'*2*i1*a2_hat'*(n20.*n21.*2.*i2)*a3_hat'*2*i3...
    + a1_hat'*2*i1*a2_hat'*2*i2*a3_hat'*(n30.*n31.*2.*i3);

Ad(1, 1:length(phik)) = (n10.*n11.*2.*i1)'*V1_hat + (n30.*n31.*2.*i3)'*V3_hat...
    + a3_hat'*((n20.*n21.*2.*i2)*2*i3'+2*i2*(n30.*n31.*2.*i3)')'*V2_hat...
    + a2_hat'*((n20.*n21.*2.*i2)*2*i3'+2*i2*(n30.*n31.*2.*i3)')*V3_hat...
    + a1_hat'*(n10.*n11.*2.*i1)*a2_hat'*2*i2*(2*i3)'*V3_hat...
    + a1_hat'*(n10.*n11.*2.*i1)*a3_hat'*2*i3*(2*i2)'*V2_hat...
    + a2_hat'*2*i2*a3_hat'*2*i3*(n10.*n11.*2.*i1)'*V1_hat...
    + a2_hat'*(n20.*n21.*2.*i2)*a3_hat'*2*i3*(2*i1)'*V1_hat...
    + a2_hat'*(n20.*n21.*2.*i2)*a1_hat'*2*i1*(2*i3)'*V3_hat...
    + a1_hat'*2*i1*a3_hat'*2*i3*(n20.*n21.*2.*i2)'*V2_hat...
    + a3_hat'*(n30.*n31.*2.*i3)*a2_hat'*2*i2*(2*i1)'*V1_hat...
    + a3_hat'*(n30.*n31.*2.*i3)*a1_hat'*2*i1*(2*i2)'*V2_hat...
    + a2_hat'*2*i2*a1_hat'*2*i1*(n30.*n31.*2.*i3)'*V3_hat; 

bd(2) = a1_hat'*(n10.*n10.*2.*i1) + a3_hat'*(n30.*n30.*2.*i3)...
    + a2_hat'*((n20.*n20.*2.*i2)*2*i3'+2*i2*(n30.*n30.*2.*i3)')*a3_hat...
    + a1_hat'*(n10.*n10.*2.*i1)*a2_hat'*2*i2*a3_hat'*2*i3...
    + a1_hat'*2*i1*a2_hat'*(n20.*n20.*2.*i2)*a3_hat'*2*i3...
    + a1_hat'*2*i1*a2_hat'*2*i2*a3_hat'*(n30.*n30.*2.*i3);

Ad(2, 1:length(phik)) = (n10.*n10.*2.*i1)'*V1_hat + (n30.*n30.*2.*i3)'*V3_hat...
    + a3_hat'*((n20.*n20.*2.*i2)*2*i3'+2*i2*(n30.*n30.*2.*i3)')'*V2_hat...
    + a2_hat'*((n20.*n20.*2.*i2)*2*i3'+2*i2*(n30.*n30.*2.*i3)')*V3_hat...
    + a1_hat'*(n10.*n10.*2.*i1)*a2_hat'*2*i2*(2*i3)'*V3_hat...
    + a1_hat'*(n10.*n10.*2.*i1)*a3_hat'*2*i3*(2*i2)'*V2_hat...
    + a2_hat'*2*i2*a3_hat'*2*i3*(n10.*n10.*2.*i1)'*V1_hat...
    + a2_hat'*(n20.*n20.*2.*i2)*a3_hat'*2*i3*(2*i1)'*V1_hat...
    + a2_hat'*(n20.*n20.*2.*i2)*a1_hat'*2*i1*(2*i3)'*V3_hat...
    + a1_hat'*2*i1*a3_hat'*2*i3*(n20.*n20.*2.*i2)'*V2_hat...
    + a3_hat'*(n30.*n30.*2.*i3)*a2_hat'*2*i2*(2*i1)'*V1_hat...
    + a3_hat'*(n30.*n30.*2.*i3)*a1_hat'*2*i1*(2*i2)'*V2_hat...
    + a2_hat'*2*i2*a1_hat'*2*i1*(n30.*n30.*2.*i3)'*V3_hat; 

bd(3) = a1_hat'*(n11.*n11.*2.*i1) + a3_hat'*(n31.*n31.*2.*i3)...
    + a2_hat'*((n21.*n21.*2.*i2)*2*i3'+2*i2*(n31.*n31.*2.*i3)')*a3_hat...
    + a1_hat'*(n11.*n11.*2.*i1)*a2_hat'*2*i2*a3_hat'*2*i3...
    + a1_hat'*2*i1*a2_hat'*(n21.*n21.*2.*i2)*a3_hat'*2*i3...
    + a1_hat'*2*i1*a2_hat'*2*i2*a3_hat'*(n31.*n31.*2.*i3);

Ad(3, 1:length(phik)) = (n11.*n11.*2.*i1)'*V1_hat + (n31.*n31.*2.*i3)'*V3_hat...
    + a3_hat'*((n21.*n21.*2.*i2)*2*i3'+2*i2*(n31.*n31.*2.*i3)')'*V2_hat...
    + a2_hat'*((n21.*n21.*2.*i2)*2*i3'+2*i2*(n31.*n31.*2.*i3)')*V3_hat...
    + a1_hat'*(n11.*n11.*2.*i1)*a2_hat'*2*i2*(2*i3)'*V3_hat...
    + a1_hat'*(n11.*n11.*2.*i1)*a3_hat'*2*i3*(2*i2)'*V2_hat...
    + a2_hat'*2*i2*a3_hat'*2*i3*(n11.*n11.*2.*i1)'*V1_hat...
    + a2_hat'*(n21.*n21.*2.*i2)*a3_hat'*2*i3*(2*i1)'*V1_hat...
    + a2_hat'*(n21.*n21.*2.*i2)*a1_hat'*2*i1*(2*i3)'*V3_hat...
    + a1_hat'*2*i1*a3_hat'*2*i3*(n21.*n21.*2.*i2)'*V2_hat...
    + a3_hat'*(n31.*n31.*2.*i3)*a2_hat'*2*i2*(2*i1)'*V1_hat...
    + a3_hat'*(n31.*n31.*2.*i3)*a1_hat'*2*i1*(2*i2)'*V2_hat...
    + a2_hat'*2*i2*a1_hat'*2*i1*(n31.*n31.*2.*i3)'*V3_hat; 
bd = bd(:);
end

if Np == 4
bp(1) = 2*a1_hat'*(n10.*n11.*i1*i2'+i1*(n20.*n21.*i2)')*a2_hat - (n20.*n21.*i2)'*a2_hat; 
Ap(1, 1:length(phik)) = 2*a1_hat'*(n10.*n11.*i1*i2'+i1*(n20.*n21.*i2)')*V2_hat...
    + 2*a2_hat'*((n20.*n21.*i2)*i1'+i2*(n10.*n11.*i1)')*V1_hat - (n20.*n21.*i2)'*V2_hat;  

bp(2) = 2*a1_hat'*(n10.*n10.*i1*i2'+i1*(n20.*n20.*i2)')*a2_hat - (n20.*n20.*i2)'*a2_hat; 
Ap(2, 1:length(phik)) = 2*a1_hat'*(n10.*n10.*i1*i2'+i1*(n20.*n20.*i2)')*V2_hat...
    + 2*a2_hat'*((n20.*n20.*i2)*i1'+i2*(n10.*n10.*i1)')*V1_hat - (n20.*n20.*i2)'*V2_hat; 

bp(3) = 2*a1_hat'*(n11.*n11.*i1*i2'+i1*(n21.*n21.*i2)')*a2_hat - (n21.*n21.*i2)'*a2_hat; 
Ap(3, 1:length(phik)) = 2*a1_hat'*(n11.*n11.*i1*i2'+i1*(n21.*n21.*i2)')*V2_hat...
    + 2*a2_hat'*((n21.*n21.*i2)*i1'+i2*(n11.*n11.*i1)')*V1_hat - (n21.*n21.*i2)'*V2_hat; 

bp = bp(:);
end

if Np == 2 && Nd == 2
    A = [];
    b = [];
elseif Np == 4 && Nd == 2
    A = Ap; 
    b = bp;
elseif Np == 2 && Nd == 4
    A = Ad;
    b = bd;
elseif Np == 4 && Nd == 4
    A = [Ap; Ad];
    b = [bp; bd];
end