1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 |
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 |