octave_packages/geometry-1.5.0/shape2d/shapecentroid.m

   1 %% Copyright (c) 2011 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
   2 %%
   3 %% This program is free software; you can redistribute it and/or modify
   4 %% it under the terms of the GNU General Public License as published by
   5 %% the Free Software Foundation; either version 3 of the License, or
   6 %% (at your option) any later version.
   7 %%
   8 %% This program is distributed in the hope that it will be useful,
   9 %% but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 %% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  11 %% GNU General Public License for more details.
  12 %%
  13 %% You should have received a copy of the GNU General Public License
  14 %% along with this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 %% -*- texinfo -*-
  17 %% @deftypefn {Function File} { @var{cm} =} shapecentroid (@var{pp})
  18 %%  Centroid of a simple plane shape defined with piecewise smooth polynomials.
  19 %%
  20 %% The shape is defined with piecewise smooth polynomials. @var{pp} is a
  21 %% cell where each elements is a 2-by-(poly_degree+1) matrix containing a pair
  22 %% of polynomials.
  23 %% @code{px(i,:) = pp@{i@}(1,:)} and @code{py(i,:) = pp@{i@}(2,:)}.
  24 %%
  25 %% The edges of the shape should not self-intersect. This function does not check for the
  26 %% sanity of the shape.
  27 %%
  28 %% @seealso{shapearea, shape2polygon}
  29 %% @end deftypefn
  30
  31 function cm = shapecentroid (shape)
  32
  33   cm = sum( cell2mat ( cellfun (@CMint, shape, 'UniformOutput', false)), 1);
  34   A = shapearea(shape);
  35   cm = cm / A;
  36
  37   [~,id] = lastwarn ('','');
  38   if strcmp (id ,'geom2d:shapearea:InvalidResult')
  39     lastwarn('Inverting centroid','geom2d:shapecentroid:InvalidResult');
  40     cm = -cm;
  41   end
  42
  43 endfunction
  44
  45 function dcm = CMint (x)
  46
  47     px = x(1,:);
  48     py = x(2,:);
  49     Px = polyint (conv(conv (-px , py) , polyder (px)));
  50     Py = polyint (conv(conv (px , py) , polyder (py)));
  51
  52     dcm = zeros (1,2);
  53     dcm(1) = diff(polyval(Px,[0 1]));
  54     dcm(2) = diff(polyval(Py,[0 1]));
  55
  56 endfunction
  57
  58 %!demo % non-convex bezier shape
  59 %! boomerang = {[ 0 -2 1; ...
  60 %!               -4  4 0]; ...
  61 %!              [0.25 -1; ...
  62 %!               0     0]; ...
  63 %!              [ 0 1.5 -0.75; ...
  64 %!               -3 3    0];
  65 %!              [0.25 0.75; ...
  66 %!               0 0]};
  67 %! CoM = shapecentroid (boomerang)
  68 %! Gcentroid = centroid(shape2polygon(boomerang))
  69 %! figure(1); clf;
  70 %! shapeplot(boomerang,'-o');
  71 %! hold on
  72 %! drawPoint(CoM,'xk;shape centroid;');
  73 %! drawPoint(Gcentroid,'xr;point centroid;');
  74 %! hold off
  75 %! axis equal
  76
  77 %!demo
  78 %! Lshape = {[0.00000   0.76635; -0.67579  -0.24067]; ...
  79 %!             [0.77976   0.76635; 0.00000  -0.91646]; ...
  80 %!             [0.00000   1.54611; 0.38614  -0.91646]; ...
  81 %!             [-0.43813   1.54611; 0.00000  -0.53032]; ...
  82 %!             [0.00000   1.10798; 0.28965  -0.53032]; ...
  83 %!             [-0.34163   1.10798; 0.00000  -0.24067]};...
  84 %! CoM = shapecentroid (Lshape)
  85 %! Gcentroid = centroid (shape2polygon (Lshape))
  86 %!
  87 %! shapeplot(Lshape,'-o');
  88 %! hold on
  89 %! drawPoint(CoM,'xk;shape centroid;');
  90 %! drawPoint(Gcentroid,'xr;point centroid;');
  91 %! hold off
  92 %! axis equal
  93
  94 %!test
  95 %! square = {[1 -0.5; 0 -0.5]; [0 0.5; 1 -0.5]; [-1 0.5; 0 0.5]; [0 -0.5; -1 0.5]};
  96 %! CoM = shapecentroid (square);
  97 %! assert (CoM, [0 0], sqrt(eps));
  98
  99 %!test
 100 %! square = {[1 -0.5; 0 -0.5]; [0 0.5; 1 -0.5]; [-1 0.5; 0 0.5]; [0 -0.5; -1 0.5]};
 101 %! square_t = shapetransform (square,[1;1]);
 102 %! CoM = shapecentroid (square_t);
 103 %! assert (CoM, [1 1], sqrt(eps));
 104
 105 %!test
 106 %! circle = {[1.715729  -6.715729    0   5; ...
 107 %!            -1.715729  -1.568542   8.284271    0]; ...
 108 %!            [1.715729   1.568542  -8.284271    0; ...
 109 %!             1.715729  -6.715729    0   5]; ...
 110 %!            [-1.715729   6.715729    0  -5; ...
 111 %!             1.715729   1.568542  -8.284271    0]; ...
 112 %!            [-1.715729  -1.568542   8.284271    0; ...
 113 %!            -1.715729   6.715729    0  -5]};
 114 %! CoM = shapecentroid (circle);
 115 %! assert (CoM , [0 0], 5e-3);
 116
 117 %!shared shape
 118 %! shape = {[-93.172   606.368  -476.054   291.429; ...
 119 %!          -431.196   637.253    11.085   163.791]; ...
 120 %!         [-75.3626  -253.2337   457.1678   328.5714; ...
 121 %!           438.7659  -653.6278    -7.9953   380.9336]; ...
 122 %!         [-89.5841   344.9716  -275.3876   457.1429; ...
 123 %!          -170.3613   237.8858     1.0469   158.0765];...
 124 %!         [32.900  -298.704   145.804   437.143; ...
 125 %!         -243.903   369.597   -34.265   226.648]; ...
 126 %!         [-99.081   409.127  -352.903   317.143; ...
 127 %!           55.289  -114.223   -26.781   318.076]; ...
 128 %!         [-342.231   191.266   168.108   274.286; ...
 129 %!           58.870   -38.083   -89.358   232.362]};
 130
 131 %!test % x-Reflection
 132 %! v = shapecentroid (shape)(:);
 133 %! T = createLineReflection([0 0 1 0]);
 134 %! nshape = shapetransform (shape, T);
 135 %! vn = shapecentroid (nshape)(:);
 136 %! assert(vn,T(1:2,1:2)*v);
 137
 138 %!test % Rotation
 139 %! v = shapecentroid (shape)(:);
 140 %! T = createRotation(v.',pi/2);
 141 %! nshape = shapetransform (shape, T);
 142 %! vn = shapecentroid (nshape)(:);
 143 %! assert(vn,v,1e-2);
 144
 145 %!test % Translation
 146 %! v = shapecentroid (shape)(:);
 147 %! nshape = shapetransform (shape, -v);
 148 %! vn = shapecentroid (nshape)(:);
 149 %! assert(vn,[0; 0],1e-2);