function f = obj_fun(X)
 % Although matlab provides better ways to overload the objective function
 % with additional data, I'll keep things simple here and use the following
 % trick to bring the two shapes into the function workspace
target = evalin('base', 'target');
c2vsim = evalin('base', 'c2vsim');

 % Because the shapefiles contain islands I split them into polygons and use
 % only the first which is the outline.
[X1, Y1] = polysplit(c2vsim.X, c2vsim.Y);
[X2, Y2] = polysplit(target.X, target.Y);

DX = X(1); % This is the x translation
DY = X(2); % This is the y translation
 % first I run the optimzation without rotation and scaling
 % So I give here some default values
if length(X) > 2
    TH = deg2rad(X(3)); % This is the rotation
    TH = 0;
if length(X) > 3
    S = X(4);
    S = 1;

 % translate shape to origin
xc = mean(X1{1,1});
yc = mean(Y1{1,1});
X1{1,1} = X1{1,1} - xc;
Y1{1,1} = Y1{1,1} - yc;

 % first scale
X1{1,1} = S*X1{1,1};
Y1{1,1} = S*Y1{1,1};

 % then rotate around the origin (0,0)
 % in case the shape is away from the origin it may be necessary to
 % translate first so that the centroid of the shape is on the origin,
 % then rotate and translate back again
XR = X1{1,1}*cos(TH) -  Y1{1,1}*sin(TH);
YR = X1{1,1}*sin(TH) +  Y1{1,1}*cos(TH);

XR = XR + xc;
YR = YR + yc;

 % and then translate according to the desicion variables
XR = XR + DX;
YR = YR + DY;

 % find their intersection
[x, y] = polybool('intersection', XR, YR, X2{1,1}, Y2{1,1});

if isempty(x)
    % if the two shapes do not intersect they are still far away therefore
    % return the distance between their barycenters plus a constant that
    % will make the objective function value greater than any possible
    % value if they are intersected.
    f = sqrt((mean(XR) - mean(X2{1,1}))^2 + (mean(YR) - mean(Y2{1,1}))^2) + 1000;
    % If they intersect compute the area that remains after we subtract the
    % two shapes. We want this area minimum. Note that we have divided the
    % area by a 10000 sq km. to make this value small.
    % We perform two subtraction operators, to avoid the case where one is
    % completely inside the other and the subtracted area is zero.
    [x1, y1] = polybool('subtraction', X2{1,1}, Y2{1,1}, XR, YR);
    [x2, y2] = polybool('subtraction', XR, YR, X2{1,1}, Y2{1,1});

    [xs, ys] = polysplit(x1, y1);
    f1 = 0;
    for i = 1:size(xs,1)
        f1 = f1 + polyarea(xs{i,1}, ys{i,1})/(10000*10000);

    [xs, ys] = polysplit(x2, y2);
    f2 = 0;
    for i = 1:size(xs,1)
        f2 = f2 + polyarea(xs{i,1}, ys{i,1})/(10000*10000);
    % we return the maximum subtracted area.
    f = max(f1, f2);


Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.