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

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

Subsurface

George Kourakos

Leave a Reply Cancel reply