Diffing " Grid More4" and "try1"

function xy = solver(kd,bx)

npts = length(kd); 

if npts == 1
    xy = bx([1 3]);
    return
end

nones=ones(npts,1); 
xy = zeros(npts,1); 
mx=max(kd(:));
dx=min(bx(2),mx); 
dy=min(bx(4),mx);

if any(kd(:)>0)
    D=dx+dy;
    d=D/30;
    while max(2,ceil(dx/d))*max(2,ceil(dy/d))<150
        d=d/2;
    end
    x=0:d:dx;
    if x(end)<dx
        x(end+1)=dx;
    end
    y=0:d:dy;
    if y(end)<dy
        y(end+1)=dy;
    end
    [x,y]=meshgrid(x,y);
    xy0=x(:).'+i*y(:).';
    k = length(xy0);
    [dum,p]=sort(-max(kd));
    kd=kd(p,p);
    S=((kd-eye(npts))>=0); % Remember positive position
    sI=find(S); % Positive index
    strainMatrix = S; % Create strainMatrix;
    M=fix(8+npts/16);
    N0=min(M,npts);
    xy(1:N0) = solverl(S(1:N0,1:N0),kd(1:N0,1:N0),dx,dy,N0);
    for N=(M+1):npts
        xy(N) = solver2(S(1:N,1:N),kd(1:N,1:N),N,xy0,...
            xy(1:N,:),k,d,dx,dy);
    end
    for N=1:N0
        xy(N) = solver2(S,kd,N,xy0,xy,k,d,dx,dy);
    end
    i0=0;
    for N=1:npts
        [XY1,XY2] = meshgrid(xy,xy.');
        strainMatrix(sI) = abs(abs(XY1(sI)-XY2(sI))-kd(sI));
        [m,i]=max(sum(strainMatrix));
        if i==i0, break; else, i0=i; end
        xy(i) = solver2(S,kd,i,xy0,xy,k,d,dx,dy);
    end
    [v,w]=sort(p);
    xy=xy(w,:);
end 
xy=xy+(bx(1)+1i*bx(3)); 
xy=[real(xy) imag(xy)];


function xy = solverl(S,kd,dx,dy,npts) 
[I,J] = find(triu(S)); 
y = kd(find(triu(S)));
k = 450+npts*15; % change k according to memory capacity.
x = rand(npts,k)*dx+i*rand(npts,k)*dy;
[tv,id] = min(sum(abs(abs(x(I,:)-x(J,:))-y(:,ones(1,k)))));
xy=x(:,id);
for iC=1:5,
   [X1,X2]=meshgrid(xy);
   % Direction of force between every point and every other point
   % * magnitude of strain
   % * whether or not S(i,j) contains a don't-care for xy(i) and xy(j)
   forces = sum((X1-X2)./(abs(X1-X2)+eps).*(kd - abs(X1-X2)).*S,2); 
   xy = xy - .05*forces;
   xy=max(zeros(npts,1),real(xy))+j*max(zeros(npts,1),imag(xy));
   xy=min(dx*ones(npts,1),real(xy))+j*min(dy*ones(npts,1),imag(xy));
end;
            

function xy = solver2(S,kd,N,xy,xy0,k,d,dx,dy) 
I = find(S(:,N)); 
if isempty(I),
  xy = xy0(N);
  return;
end 
z = kd(:,N); 
kones=ones(1,k); 
[XY1,XY2]=meshgrid(xy,xy0(I));
[tv,id]=min(sum(abs(abs(XY1-XY2)-z(I,kones)))); 
x1=real(xy(id)); 
y1=imag(xy(id));
x=max(0,x1-d):d/10:min(dx,x1+d); 
y=max(0,y1-d):d/10:min(dy,y1+d); 
[x,y]=meshgrid(x,y);
xy=(x+i*y); 
xy=xy(:).'; 
k = length(xy); 
kones=ones(1,k); 
[XY1,XY2]=meshgrid(xy,xy0(I));
[tv,id]=min(sum(abs(abs(XY1-XY2)-z(I,kones)))); 
xy=xy(id);