function [U]=SolveVIC(Uini,iscale,nmod)
check=1;
ttic=cputime;
load(fullfile('TMP',sprintf('sample0_%d',1-1)),'sizeim');
sizeim0=sizeim;
load(fullfile('TMP','params'),'param');
param0=param;
roi=param0.roi;
reverse=0;
if isfield(param0,'reverse_image')
    reverse=param0.reverse_image;
end
if iscell(param0.deformed_image)
    nim=length(param0.deformed_image);
else
    nim=1;
end
if isfield(param0,'restart')
    restart=param0.restart;
else
    restart=1;
end
maxiter=param0.iter_max;
conv=param0.convergance_limit;
load(fullfile('TMP',sprintf('%d_params',nmod)),'param');
disp(sprintf('Starting resolution for scale %d...',iscale));

pscale=2^(1-1);
load(fullfile('TMP',sprintf('sample0_%d',iscale-1)),'im0','sizeim','mean0','std0','nband');
dynamic=max(im0(:))-min(im0(:));
ndim=length(sizeim);

if check
    load(fullfile('TMP',sprintf('sample0_%d',1-1)),'lso','ls1','nx','ny');
    load(fullfile('TMP',sprintf('sample0_%d',iscale-1)),'on');
load(fullfile('TMP',sprintf('%d_phi_%d',nmod,(iscale-1))),'phii','ls1i');
[xon,yon]=ind2sub(sizeim,nband(on));
nxon=-nx(nband(on));
nyon=-ny(nband(on));
end

load(fullfile('TMP',sprintf('%d_phix_%d',nmod,(iscale-1))),'phix','wdetJ','Xi','Yi');
load(fullfile('TMP',sprintf('%d_phiy_%d',nmod,(iscale-1))),'phiy');
load(fullfile('TMP',sprintf('%d_mask_%d',nmod,iscale-1)),'mask');
load(fullfile('TMP',sprintf('%d_phidf_%d',nmod,iscale-1)),'phidf');
load(fullfile('TMP',sprintf('%d_operator_%d',nmod,iscale-1)),'M','R');
%R=0*R;
M=M+R;
im0=im0-mean0;
im1=0;
U=Uini;
xi=roi(1)-1+(1:sizeim(1));
yi=roi(3)-1+(1:sizeim(2));
fid=fopen(fullfile('TMP',sprintf('%d_error_%d.mat',nmod,iscale-1)),'w');
fwrite(fid,min(255,dynamic));

if restart
    nmax=nim;
else
    if iscale==1
        nmax=nim;
    else
        if nim<2
            nmax=nim;
        else
            nmax=2;
        end
    end
end
%nmax=nim;

for ijm=1:nmax

display(sprintf('Image %d/%d',ijm,nmax));
res=1;
    ii=1;
    if restart
        Ui=Uini(:,ijm);
    else
                if ijm<3
                    Ui=Uini(:,ijm);
                else
                            DU=U(:,ijm-1)-U(:,ijm-2);
        
                Ui=U(:,ijm-1)+DU;
                end
%         if iscale==param.nscale
%             Ui=U(:,max(1,ijm-1));
%         else
%             Ui=Uini(:,ijm);
%         end
    end
    if check
    fig=figure;
    title(sprintf('Image %d/%d',ijm,nmax));
    hold on
    end
    while ( res>conv && ii< maxiter)
        Ux=phix*Ui;
        Uy=phiy*Ui;
        NestedResidual(ijm);
        disc=mask*disc;
        merror=std(disc(nband))/dynamic;

        F=phidf'*(wdetJ*disc)-R*Ui;
        if isfield(param0,'solver')
            solv=param0.solver;

            if strcmp(solv,'lu')
                if ii==1
                    [LT,UT]=lu(sparse(M));
                end
                dU=UT\(LT\F);
            elseif strcmp(solv,'chol')
                if ii==1
                    R=chol(sparse(M));
                end
                dU=R\(R'\F);
            end
        else
            dU=(M\F);
        end
        disp(sprintf('At iteration # %d',ii));
        disp(sprintf('Discrepancy wrt dyn. =%6.2f %%',merror*100));
%        if ii==1
            nU0=norm(Ui+dU);
%        else
            res=norm(dU)/nU0;
            disp(sprintf('|dU|=%f',res));

%        end
        Ui=Ui+dU;
        if check
if exist('h')&&ii>1
    delete(h)
else
    h=plot(xon+roi(1)-1,yon+roi(3)-1,'r.','LineWidth',2);
%   [C,h]=contour(xi,yi,lso',[0 0],'Color','red','LineWidth',2);  
       print ('-djpeg', fullfile('FIG',sprintf('%s-scale-%d-iter-%d.jpg',param0.result_file,iscale,ii-1)));
%       print ('-depsc2', fullfile('FIG',sprintf('%s-scale-%d-iter-%d.eps',param0.result_file,iscale,ii-1)));
%pause
 delete(h)
 pause(0.1)
end
% uv=(phii*Ui);
% uvi=interp1(ls1i,uv,ls1(:),'linear','extrap');
%  cont=lso+reshape(uvi,sizeim);
%  [C,h]=contour(xi,yi,cont',[0 0],'Color','blue','LineWidth',2);
    h=plot(xon+roi(1)-1+nxon.*(phii*Ui),yon+roi(3)-1+nyon.*(phii*Ui),'r.','LineWidth',2);
        print ('-djpeg', fullfile('FIG',sprintf('%s-scale-%d-iter-%d.jpg',param0.result_file,iscale,ii)));
 %       print ('-depsc2', fullfile('FIG',sprintf('%s-scale-%d-iter-%d.eps',param0.result_file,iscale,ii)));
 pause(0.1)
        end
 %         figure
%         imagesc(reshape(phi*Ui,sizeim))
%         colorbar
%         pause
        ii=ii+1;
    end
    if check 
        close(fig)
 pause(0.1)
    end
    U(:,ijm)=Ui;
    if dynamic>255
        disc=255*disc/dynamic;
    end
    fwrite(fid,abs(disc));
end
disp(sprintf('Enlapsed time for resolution = %6.2f s',cputime -ttic));
fclose(fid);
    function NestedResidual(kkk)
        if ii==1
            if nim==1
                fildef=param0.deformed_image;
                im1=double(readim(fildef));
            else
                fildef=param0.deformed_image{kkk};
                im1=double(readim(fildef));
             end
            if length(size(im1))==3
                im1=mean(im1,3);
%                im1=0.299 * im1(:,:,1)  + 0.587 *   im1(:,:,2)+ 0.114 *  im1(:,:,3);
            end
                if reverse
                    im1=im1';
                end
        end
%         figure
%         imagesc(im1)
        %         if iscale==2
%                  figure
%                  imagesc(reshape(Ux,sizeim*pscale))
%                  pause
        %         end
%        disc=mexInterpSpline((Xi+Ux+roi(1)-1),(Yi+Uy+roi(3)-1),im1);
                    sizeim1=size(im1);
if check
    if ii==1
                    h2=imagesc(im1');
                    colormap(gray)    
axis off;
axis xy;
axis image;
%delete(h2)
    end
end
        disc=(mexInterpSpline(min(sizeim1(1)-2,max(3,(Xi+Ux+roi(1)-1))),min(sizeim1(2)-2,max(3,(Yi+Uy+roi(3)-1))),im1));

        %                  figure
%                  imagesc(reshape(disc,sizeim*pscale))

%         if iscale>1
%             disc=reshape(disc,sizeim*pscale);
%             NestedCoarseImage(pscale);
%         end

        mean1=mean(disc(nband));
        std1=std(disc(nband));
%std1=max(disc(nband))-min(disc(nband));
%std1=std(disc(:));
%std1=sqrt(std(disc(nband))^2*numel(nband)/prod(sizeim));
        disc=disc-mean1;
%std1=max(im1(:))-min(im1(:));
        st=std0/std1;
        disc=(im0(:)-st*disc(:));


        function NestedCoarseImage(scale)

            imsiz0=size(disc);
            imsiz1=floor(imsiz0/2);
            nn=2*imsiz1;
            ndim=length(nn);
            if ndim==2
                disc=disc(1:nn(1),1:nn(2));
            elseif ndim==3
                disc=disc(1:nn(1),1:nn(2),1:nn(3));
            end

            disc=reshape(disc,scale,prod(nn)/scale);
            disc=mean(disc,1);
            nn(1)=nn(1)/scale;
            disc=reshape(disc,nn);

            if ndim==2
                disc=disc';
                disc=reshape(disc,scale,prod(nn)/scale);
                disc=mean(disc,1);
                nn(2)=nn(2)/scale;
                disc=reshape(disc,nn([2,1]));
                disc=disc';
                %toc();
            elseif ndim==3
                disc=permute(disc,[2,3,1]);
                disc=reshape(disc,scale,prod(nn)/scale);
                disc=mean(disc,1);
                nn(2)=nn(2)/scale;
                disc=reshape(disc,nn([2,3,1]));
                disc=permute(disc,[3,1,2]);

                disc=permute(disc,[3,1,2]);
                disc=reshape(disc,scale,prod(nn)/scale);
                disc=mean(disc,1);
                nn(3)=nn(3)/scale;
                disc=reshape(disc,nn([3,1,2]));
                disc=permute(disc,[2,3,1]);
            end

        end
    end
end