reconstruct.m

I = I(:);
    reconstruct = zeros(size(I));
    for j = 0:order-1
        reconstruct = reconstruct + ((I-reconstruct)\vp(:,size(vp,2)-j))*vp(:,size(vp,2)-j);
    end