function run_vic_job(paramo,modelo)
nmod=0;
param=paramo;
filres=paramo.result_file;
filreso=strrep(filres,'.res','');
nscale=modelo.nscale;
filo=param.reference_image;
[~,~,ext]=fileparts(filo);
if iscell(param.deformed_image)
fildefs=param.deformed_image;
fildefs=[filo,fildefs];
images=1:size(fildefs,2);
else
if strcmp(filo,param.deformed_image)
fildefs=filo;
images=1;
else
fildefs={filo,param.deformed_image};
images=1:size(fildefs,2);
end
end
param.deformed_image=fildefs;
filim=strrep(param.result_file,'.res','-vi.mat');
param.reference_image=filim;
[Yi,Xi]=meshgrid(1:param.sizeim(2),1:param.sizeim(1));
xo=[];
yo=[];
conn=[];
Ux=[];
Uy=[];
dU=[];
ddU=[];
Up=[];
nmax=0;
curv=cell(size(modelo.zone,2),1);
fleche=cell(size(modelo.zone,2),1);
s=cell(size(modelo.zone,2),1);
for iz=1:size(modelo.zone,2)
zone=modelo.zone(:,iz);
switch zone{4}
case 1
xyp=zone{2};
xyp=xyp*[1;1i];
roi=round([min(real(xyp)),max(real(xyp)),min(imag(xyp)),max(imag(xyp))]);
try
model.starting_point=zone{10};
catch
model.starting_point=[real(xyp(1)),imag(xyp(1))];
end
mx=1;
my=1;
model.continuity='periodic';
case 3
xyp=zone{5};
Zi=Xi+1i*Yi-(xyp(1)+1i*xyp(2));
lso=abs(Zi);
ls1=angle(Zi)*xyp(3);
nx=-real(Zi)./(max(lso,1));
ny=-imag(Zi)./(max(lso,1));
lso=-(lso-xyp(3));
roi=[1,param.sizeim(1),1,param.sizeim(2)];
mx=1;my=1;
model.continuity='periodic';
case 4
model.continuity='open';
xyp=zone{2};
xyp=xyp*[1;1i];
Zi=Xi+1i*Yi-xyp(1);
seg=diff(xyp([1;numel(xyp)]));
n=(-1i*seg)/abs(seg);
t=(seg)/abs(seg);
lso=(real((Zi)'*n))';
ls1=(real((Zi)'*t))';
nx=ones(param.sizeim(1:2))*real(n);
ny=ones(param.sizeim(1:2))*imag(n);
roi=round([min(real(xyp)),max(real(xyp)),min(imag(xyp)),max(imag(xyp))]);
mx=(abs(angle(n))<pi/4)||(abs(angle(n))>3*pi/4);
my=(abs(angle(n))>pi/4)&&(abs(angle(n))<3*pi/4);
end
if zone{4}>1
if zone{1}<0
im=readim(filo);
mim=mean(im((lso>0)&(lso<0.5*max(lso(:)))));
mam=mean(im((lso<0)&(lso>-0.5*max(lso(:)))));
if mim<mam
lso=-lso;
nx=-nx;ny=-ny;
end
model.contour_type='edge';
else
model.contour_type='line';
model.line_thickness=zone{9};
end
save(filim,'lso','ls1','nx','ny','-v7.3');
else
model.contour_type='line';
model.line_thickness=zone{9};
param.reference_image=filo;
end
thickness=zone{9};
tau=zone{8};
roi([1])=max(roi([1])-mx*tau*2^(nscale-1),1);
roi([3])=max(roi([3])-my*tau*2^(nscale-1),1);
roi([2])=min(roi([2])+mx*tau*2^(nscale-1),param.sizeim(1));
roi([4])=min(roi([4])+my*tau*2^(nscale-1),param.sizeim(2));
param.roi=roi;
param.transition_length=tau;
param.onflight=0;
LoadParameters(param);
model.nscale=nscale;
model.degree=zone{7};
model.basis='vic-nurbs';
model.Nelems=zone{6};
model.closed=zone{4}==3;
LoadParameters(model,nmod);
if zone{4}==1
InitializeVICContour(nmod,1,50,0);
else
if zone{4}==4&&size(xyp,1)>2
%pour Rana
si=cumsum(abs(diff(xyp)));
si=[0;si(:)];
ls1i=linspace(0,max(si(:)),1000);
if 0
xcoefs=polyfit(si,real(xyp),3);
ycoefs=polyfit(si,imag(xyp),3);
xi=polyval(xcoefs,ls1i);
yi=polyval(ycoefs,ls1i);
dxcoefs=(numel(xcoefs):-1:1).*xcoefs;
dxcoefs=dxcoefs(1:end-1);
dycoefs=(numel(ycoefs):-1:1).*ycoefs;
dycoefs=dycoefs(1:end-1);
dxi=polyval(dxcoefs,ls1i);
dyi=polyval(dycoefs,ls1i);
else
xi=interp1(si,real(xyp),ls1i,'spline');
yi=interp1(si,imag(xyp),ls1i,'spline');
dxi=gradient(xi);
dyi=gradient(yi);
end
ni=(-1i*(dxi+1i*dyi));
ni=ni(:)./abs(ni(:));
xyon=xi(:)+1i*yi(:);
figure
plot(xyon)
hold on
plot(xyp)
ini=inpolygon(real(xyon),imag(xyon),roi([1,2,2,1,1]),roi([3,3,4,4,3]));
xyon=xyon(ini);
ni=-ni(ini);
[lson,ls1n,nmesh]=ComputeLevelSetFromPoints(nmod,1,real(xyon(:))-roi(1)+1,imag(xyon(:))-roi(3)+1,real(ni(:)),imag(ni(:)),0);
lson(~nmesh)=100000;
figure
imagesc(lson)
lso(roi(1):roi(2),roi(3):roi(4))=lson;
ls1(roi(1):roi(2),roi(3):roi(4))=ls1n;
nxi=FDgradient(lson,1);
nyi=FDgradient(lson,2);
nx(roi(1):roi(2),roi(3):roi(4))=nxi;
ny(roi(1):roi(2),roi(3):roi(4))=nyi;
save(filim,'lso','ls1','nx','ny','-v7.3');
end
VirtualImage(nmod);
end
U1=[];
VICMeshes(nmod);
LoadMask(nmod);
for iscale=(nscale):-1:1
CreateVICBasisFunction(iscale,nmod);
ComputeGradFPhi(iscale,nmod);
AssembleCorrelationOperator(iscale,nmod);
Uini=InitializeVICSolution(U1,iscale,nmod);
[U1]=SolveVIC(Uini,iscale,nmod);
end
Up=[Up;U1];
load(fullfile('TMP',sprintf('sample0_%d',0)),'lso','ls1','nx','ny','nband','on');
load(fullfile('TMP',[num2str(nmod),'_phi_0']),'phii','dphii','ddphii','sizeim');
[xon,yon]=ind2sub(sizeim,nband(on));
ls1on=ls1(nband(on));
nxon=-nx(nband(on));
nyon=-ny(nband(on));
[si,ind]=sort(ls1(nband(on)));
ls1on=ls1on(ind);
s{iz}=si;
phii=phii(ind,:);
dphii=dphii(ind,:);
ddphii=ddphii(ind,:);
fleche{iz}=phii*U1;
curv{iz}=ddphii*U1;
xon=xon(ind)+param.roi(1)-1;
yon=yon(ind)+param.roi(3)-1;
nxon=nxon(ind);
nyon=nyon(ind);
uxi=diag(sparse(nxon))*(phii*U1);
uyi=diag(sparse(nyon))*(phii*U1);
nn=length(xon);
%xon=xon+uxi(:,1);
%yon=yon+uyi(:,1);
%uxi(:,1)=0;
%uyi(:,1)=0;
ddui=ddphii*U1;
fid=fopen([filreso,sprintf('-zone-%02d.csv',iz)],'w');
fprintf(fid,'Result file;%s\n',param.result_file);
fprintf(fid,'Zone gage;%d\n',iz);
fprintf(fid,'Position s [pixel];X [pixel];Y [pixel];Curvature [1/pixel]\n');
fprintf(fid,'Image file;');
if size(U1,2)>1
for iim=1:size(U1,2)
fprintf(fid,'"%s";;;',param.deformed_image{iim});
end
else
fprintf(fid,'"%s";;;',param.deformed_image);
end
fprintf(fid,'\n');
fprintf(fid,'Image number;');
for iim=1:size(U1,2)
fprintf(fid,'%d;;;',iim);
end
fprintf(fid,'\n');
for ix=1:length(xon)
fprintf(fid,'%f;',ls1on(ix));
for iim=1:size(U1,2)
fprintf(fid,'%f;%f;%f;',xon(ix)+uxi(ix,iim),yon(ix)+uyi(ix,iim),ddui(ix,iim));
end
fprintf(fid,'\n');
end
switch model.contour_type
case 'line'
xon=[xon-0.5*(tau+thickness)*nxon;xon+0.5*(tau+thickness)*nxon];
yon=[yon-0.5*(tau+thickness)*nyon;yon+0.5*(tau+thickness)*nyon];
case 'edge'
xon=[xon-(tau+thickness)*nxon;xon+0*(tau+thickness)*nxon];
yon=[yon-(tau+thickness)*nyon;yon+0*(tau+thickness)*nyon];
end
uxi=[uxi;uxi];
uyi=[uyi;uyi];
Ux=[Ux;uxi];
Uy=[Uy;uyi];
ddU=[ddU;ddui;ddui];
xo=[xo;xon];
yo=[yo;yon];
conni=repmat([1,2,2+nn,1+nn],nn-1,1)+repmat((0:nn-2)',1,4);
conn=[conn;conni+nmax];
nmax=max(conn(:));
end
paramo.deformed_image=param.deformed_image;
param=paramo;
model=modelo;
Nnodes=[length(xo),1,1];
Nelems=[size(conn,1),1,1];
elt=repmat(4,prod(Nelems),1);
U=[Ux;Uy];
delete([fullfile('VTK',filreso),'-0*.vtk']);
ExportVTKField(filreso,images,xo,yo,elt,conn,{U,ddU},{'Displacement','Curvature'});
rint=false;
ng=0;
rflag=0;
save(filres,'U','Up','s','curv','fleche','ddU','Nnodes','Nelems','xo','yo','param','model','nmod','conn','elt','rint','ng','rflag','-v7.3');
end
