function xy = solver(kd,bx)
npts = size(kd,1);
dx=bx(2);
dy=bx(4);
D=dx+1i*dy;
lD=abs(D);
mxmx=max(kd(:));
kd(kd>lD)=lD;
xy = zeros(npts,1);
mx=max(kd-eye(npts));
ismart=find(mx>=0);
if ~isempty(ismart)&dx+dy>1e-3
npts1=length(ismart);
if npts1==2
p=[0;D/lD*kd(ismart(1),ismart(2))];
xy=zeros(npts,1);
xy(ismart)=p;
xy=[real(xy) imag(xy)];
return
end
ipts=1:npts1;
kd=kd(ismart,ismart);
S=(kd-eye(npts1))>=0;
p=zeros(npts1,1);
if npts1>5
p = solvecomplex(kd,S,dx,dy,p,0,npts1,inf);
xy(ismart)=p;
else
pm=D/2;
zerolen=lD/1e6;
%k=kd;
%k(k<0)=0;
%mnd=mean(sum(k))/sqrt(npts1);
%dx=min(mnd,dx);
%dy=min(mnd,dy);
if npts1==2
xy(ismart(2))=D*kd(1,2)/lD;
else
k=kd-eye(npts1);
S=k>=0;
[mx,imx]=max(k);
[mx1,jmx]=max(mx);
i1=imx(jmx);
j=ipts(jmx);
i=ipts(i1);
n=3;
if mx1<=dx
p(i)=pm+mx1/2;
p(j)=pm-mx1/2;
elseif mx1<=dy
p(i)=pm+1i*mx1/2;
p(j)=pm-1i*mx1/2;
else
p1=D*min(1,mx1/lD)/2;
p(i)=pm+p1;
p(j)=pm-p1;
end
if jmx==1
if i1>2
ipts(i1)=ipts(2);
ipts(2)=i;
end
else
if i1~=1
ipts(i1)=ipts(1);
ipts(1)=i;
end
if jmx~=2
ipts(jmx)=ipts(2);
ipts(2)=j;
end
end
contloop=1;
while contloop>0&n<=npts1
sS=sum(S(ipts(1:n-1),ipts(n:end)),1);
[n1,i1]=max(sS);
i2=find(sS==n1);
[mx1,jmx]=max(sum(kd(ipts(1:n-1),ipts(i2+n-1)),1));
i1=i2(jmx);
i=ipts(i1+n-1);
if i>n
ipts(i1+n-1)=ipts(n);
ipts(n)=i;
end
if n1
l=find(S(ipts(1:n-1),i));
j=ipts(l);
if n1==2
dp=diff(p(j));
d1=abs(dp);
if d1<=zerolen
n1=1;
l=l(1);
j=j(1);
else
d2=kd(j(1),i);
d3=kd(j(2),i);
if d2<=zerolen
p(i)=p(j(1));
n1=0;
elseif d3<=zerolen
p(i)=p(j(2));
n1=0;
elseif d1>=d2+d3
% contloop=0;
p(i)=p(j(1))+dp/d1*d2; %verbeteren?
n1=0;
elseif d2>=d1+d3
% contloop=0;
p(i)=p(j(2))+dp/d1*d3; %(?)
n1=0;
elseif d3>=d1+d2
% contloop=0;
n1=0;
p(i)=p(j(1))-dp/d1*d2; %?
end
end
end
switch n1
case 0
% do nothing (already processed)
case 1
% contloop=0;
% break;
% Now the direction is taken to have a point closest
% to the center. This is not always(/often?) best.
if abs(p(j)-pm)<zerolen
p(i)=pm+D/lD*k(j,i);
else
p1=pm-p(j);
p(i)=p(1)+p1/abs(p1)*k(j,i);
end
p(i)=min(dx,max(0,real(p(i))))+1i*min(dy,max(0,imag(p(i))));
case 2
p(i)=calctriang(d1,d2,d3,p(j(1)),dp,pm);
otherwise
P=meshgrid(p(j));
dP=abs(P-P.');
[I,J]=find(triu(dP));
IJ=(I-1)*n1+J;
[d,id]=sort(-dP(IJ));
I=I(id);
J=J(id);
n2=min(6,n1);
pp=zeros(n2,2);
ok=zeros(n2,1);
for k=1:n2
d1=-d(k);
d2=kd(j(I(k)),i);
d3=kd(j(J(k)),i);
p1=p(j(I(k)));
p2=p(j(J(k)));
if abs((d1*d1+d2*d2-d3*d3)/2/d1/d2)<1
ok(k)=1;
[pp(k,1),pp(k,2)]=calctriang(d1,d2,d3,p1,p2-p1,pm);
end
end
iok=find(ok);
if isempty(iok)
contloop=0;
elseif length(iok)==1
if abs(pp(iok,1)-pm)<=abs(pp(iok,2)-pm)
p(i)=pp(iok);
else
p(i)=pp(iok);
end
contloop=-1;
else
pp=pp(iok,:);
pp=pp(:);
[P1,P2]=meshgrid(p(j),pp);
dP=sum(abs(abs(P1-P2)-kd(i*ones(length(pp),1),j)),2);
[P,n1]=min(dP);
p(i)=pp(n1);
if P>0.1
contloop=-1;
end
end
end
else % n1==0
% fprintf('!!!!4\n');
% p(i)=pm; % terug zoeken naar twee elementen?
contloop=0;
% vermits er al getest is of er gescheiden blokken
% zijn, en dat de "dummy's" er al uit zijn,
% zou dit niet kunnen voorkomen.
end
if contloop
p(ipts(1:n))=limitpoint(p(ipts(1:n)),dx,dy);
n=n+1;
end
end % place points
% [n npts1]
if n<=npts1|1 % (doing something else if all points a set ?)
n=n-1;
i=ipts(1:n);
pi=p(i);
mnx=min(real(pi));
mxx=max(real(pi));
mny=min(imag(pi));
mxy=max(imag(pi));
p(i)=pi+(pm-((mxx+mnx)+1i*(mxy+mny))/2);
kd1=kd(ipts,ipts);
S1=S(ipts,ipts);
p=p(ipts);
s1=findstrain(S1,p,kd1);
s2=findstrain(S1,zeros(npts1,1),kd1);
if s1>s2
n=0;
s1=s2;
end
p(ipts) = solvecomplex(kd1,S,dx,dy,p(ipts),n,npts1,s1);
end
xy(ismart) = p;
end
end
end
xy=[real(xy) imag(xy)];
function [p2,p22]=calctriang(d1,d2,d3,p0,dp,pm)
g=acos((d1*d1+d2*d2-d3*d3)/2/d2/d1);
p21=p0+exp( 1i*g)*dp/d1*d2;
p22=p0+exp(-1i*g)*dp/d1*d2;
if nargout>1
p2=p21;
elseif abs(p21-pm)<abs(p22-pm)
p2=p21;
else
p2=p22;
end
function p=limitpoint(p0,dx,dy)
p=p0;
rp2=real(p(end));
ip2=imag(p(end));
if rp2<0|rp2>dx
dtot=max(real(p))-min(real(p));
% (ook konstant bij te houden)
if dtot<dx
if rp2<0
p=p-max(rp2,dtot-dx);
p(end)=max(0,real(p(end)))+1i*imag(p(end));
else
p=p-min(rp2-dx,dx-dtot);
p(end)=min(dx,real(p(end)))+1i*imag(p(end));
end
elseif rp2<0
p(end)=1i*imag(p(end));
elseif rp2>dx
p(end)=dx+1i*imag(p(end));
end
end
if ip2<0|ip2>dy
dtot=max(imag(p))-min(imag(p));
% (ook konstant bij te houden)
if dtot<dy
if ip2<0
p=p-1i*max(ip2,dtot-dy);
p(end)=real(p(end))+1i*max(0,real(p(end)));
else
p=p-1i*min(ip2-dy,dy-dtot);
p(end)=real(p(end))+1i*min(dy,imag(p(end)));
end
elseif ip2<0
p(end)=real(p(end));
elseif ip2>dy
p(end)=real(p(end))+1i*dy;
end
end
function xy = solvecomplex(kd,S,dx,dy,xy0,npt0,npts,olds)
D=dx+1i*dy;
pm=D/2;
xy = xy0(:).';
xy(npt0+1:end)=pm;
d=0.1*max(dx,dy);
x=0:d:dx;
y=0:d:dy;
[xx,yy] = ndgrid(x,1i*y);
[nx,ny]=size(xx);
XX=[0;xx(:)+yy(:)];
nX=length(XX);
xOnes=ones(nX,1);
xZeros=zeros(nX,1);
start=1;
iones=ones(1,npts);
sCut = 1;
[dum,p]=sort(-max(kd));
S=S(p,p);
kd=kd(p,p);
[v,w]=sort(p);
for iter=1:7
for index = npt0+1:npts
I = find(S(:,index));
if ~isempty(I)
XX(1)=xy(index);
st = sum(abs(abs(XX(:,ones(size(I)))-xy(xOnes,I))-kd(xZeros+index,I)),2);
[null,minloc] = min(st);
xy(index) = XX(minloc);
end
end
[s,M]=findstrain(S,xy,kd);
if((abs(s-olds)/(s+1)) < 0.002) | s < sCut
break;
end
olds = s;
npt1=0;
end
xy=xy.';
xy=xy(w,:);
kd=kd(w,w);
S=S(w,w);
olds = 1.0e50;
oldxy = xy;
alpha = 1;
for iter=1:40
[s,M]=findstrain(S,xy,kd);
if( iter > 5 & abs(s-olds) / (s+1) < 0.0001 ) | s < sCut
xy = oldxy;
break;
end
if( s > olds )
xy = oldxy;
break;
end
oldxy = xy;
olds = s;
for indx = 1:npts
I = find(S(:,indx)>0);
dX=xy(I)-xy(indx);
Dx=mean(dX.*(1-kd(I,indx)./(abs(dX)+eps)));
% Dx=mean(dX./(abs(dX)+eps) .* M(I,indx));
xy(indx) = xy(indx) + alpha*Dx;
xy(indx) = min(dx,max(0,real(xy(indx))))+1i*min(dy,max(0,imag(xy(indx))));
end
end
alpha = 0.1;
ngrad = 30;
obj = 1.0e20*ones(ngrad,1);
xybig = zeros(npts,ngrad);
[objnew,strainMatrix]=findstrain(S,xy,kd);
obj(1) = objnew;
xybig(:,1) = xy;
gradient = zeros(npts,1);
for ij=2:ngrad
[s,sM,dM]=findstrain(S,xy,kd);
gradient(:) = sum(dM./(abs(dM)+eps).*sign(sM));
xy = xy - alpha * gradient;
xy = min(max(real(xy),0),dx) + 1j * min(max(imag(xy),0),dy);
[objnew,strainMatrix]=findstrain(S,xy,kd);
obj(ij) = objnew;
xybig(:,ij) = xy;
if( obj(ij) > obj(ij-1) )
xy = xybig(:,ij-1);
alpha = alpha / 10;
else
alpha = alpha * 2;
end
[bestobj,index] = min(obj);
if( ij > 3 & index < 2 ) | bestobj < sCut
break
end
end
[bestobj,index] = min(obj);
xy = xybig(:,index);
function [s,strainMatrix,dX]=findstrain(S,xy,kd)
sI=find(S);
strainMatrix=S;
dX=S;
[X1,X2]=meshgrid(xy);
dX(sI)=X1(sI)-X2(sI);
strainMatrix(sI) = abs(dX(sI))-kd(sI);
s = sum(abs(strainMatrix(sI)));
|
