# obj_fun.m

```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
else
TH = 0;
end
if length(X) > 3
S = X(4);
else
S = 1;
end

% 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;
else
% 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);
end

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

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