X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Fgeometry-1.5.0%2Fshape2d%2Fshapecentroid.m;fp=octave_packages%2Fgeometry-1.5.0%2Fshape2d%2Fshapecentroid.m;h=8094b0b1967488878959cdebd27df9482a0c8bd1;hp=0000000000000000000000000000000000000000;hb=c880e8788dfc484bf23ce13fa2787f2c6bca4863;hpb=1705066eceaaea976f010f669ce8e972f3734b05 diff --git a/octave_packages/geometry-1.5.0/shape2d/shapecentroid.m b/octave_packages/geometry-1.5.0/shape2d/shapecentroid.m new file mode 100644 index 0000000..8094b0b --- /dev/null +++ b/octave_packages/geometry-1.5.0/shape2d/shapecentroid.m @@ -0,0 +1,149 @@ +%% Copyright (c) 2011 Juan Pablo Carbajal +%% +%% This program is free software; you can redistribute it and/or modify +%% it under the terms of the GNU General Public License as published by +%% the Free Software Foundation; either version 3 of the License, or +%% (at your option) any later version. +%% +%% This program is distributed in the hope that it will be useful, +%% but WITHOUT ANY WARRANTY; without even the implied warranty of +%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +%% GNU General Public License for more details. +%% +%% You should have received a copy of the GNU General Public License +%% along with this program; if not, see . + +%% -*- texinfo -*- +%% @deftypefn {Function File} { @var{cm} =} shapecentroid (@var{pp}) +%% Centroid of a simple plane shape defined with piecewise smooth polynomials. +%% +%% The shape is defined with piecewise smooth polynomials. @var{pp} is a +%% cell where each elements is a 2-by-(poly_degree+1) matrix containing a pair +%% of polynomials. +%% @code{px(i,:) = pp@{i@}(1,:)} and @code{py(i,:) = pp@{i@}(2,:)}. +%% +%% The edges of the shape should not self-intersect. This function does not check for the +%% sanity of the shape. +%% +%% @seealso{shapearea, shape2polygon} +%% @end deftypefn + +function cm = shapecentroid (shape) + + cm = sum( cell2mat ( cellfun (@CMint, shape, 'UniformOutput', false)), 1); + A = shapearea(shape); + cm = cm / A; + + [~,id] = lastwarn ('',''); + if strcmp (id ,'geom2d:shapearea:InvalidResult') + lastwarn('Inverting centroid','geom2d:shapecentroid:InvalidResult'); + cm = -cm; + end + +endfunction + +function dcm = CMint (x) + + px = x(1,:); + py = x(2,:); + Px = polyint (conv(conv (-px , py) , polyder (px))); + Py = polyint (conv(conv (px , py) , polyder (py))); + + dcm = zeros (1,2); + dcm(1) = diff(polyval(Px,[0 1])); + dcm(2) = diff(polyval(Py,[0 1])); + +endfunction + +%!demo % non-convex bezier shape +%! boomerang = {[ 0 -2 1; ... +%! -4 4 0]; ... +%! [0.25 -1; ... +%! 0 0]; ... +%! [ 0 1.5 -0.75; ... +%! -3 3 0]; +%! [0.25 0.75; ... +%! 0 0]}; +%! CoM = shapecentroid (boomerang) +%! Gcentroid = centroid(shape2polygon(boomerang)) +%! figure(1); clf; +%! shapeplot(boomerang,'-o'); +%! hold on +%! drawPoint(CoM,'xk;shape centroid;'); +%! drawPoint(Gcentroid,'xr;point centroid;'); +%! hold off +%! axis equal + +%!demo +%! Lshape = {[0.00000 0.76635; -0.67579 -0.24067]; ... +%! [0.77976 0.76635; 0.00000 -0.91646]; ... +%! [0.00000 1.54611; 0.38614 -0.91646]; ... +%! [-0.43813 1.54611; 0.00000 -0.53032]; ... +%! [0.00000 1.10798; 0.28965 -0.53032]; ... +%! [-0.34163 1.10798; 0.00000 -0.24067]};... +%! CoM = shapecentroid (Lshape) +%! Gcentroid = centroid (shape2polygon (Lshape)) +%! +%! shapeplot(Lshape,'-o'); +%! hold on +%! drawPoint(CoM,'xk;shape centroid;'); +%! drawPoint(Gcentroid,'xr;point centroid;'); +%! hold off +%! axis equal + +%!test +%! 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]}; +%! CoM = shapecentroid (square); +%! assert (CoM, [0 0], sqrt(eps)); + +%!test +%! 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]}; +%! square_t = shapetransform (square,[1;1]); +%! CoM = shapecentroid (square_t); +%! assert (CoM, [1 1], sqrt(eps)); + +%!test +%! circle = {[1.715729 -6.715729 0 5; ... +%! -1.715729 -1.568542 8.284271 0]; ... +%! [1.715729 1.568542 -8.284271 0; ... +%! 1.715729 -6.715729 0 5]; ... +%! [-1.715729 6.715729 0 -5; ... +%! 1.715729 1.568542 -8.284271 0]; ... +%! [-1.715729 -1.568542 8.284271 0; ... +%! -1.715729 6.715729 0 -5]}; +%! CoM = shapecentroid (circle); +%! assert (CoM , [0 0], 5e-3); + +%!shared shape +%! shape = {[-93.172 606.368 -476.054 291.429; ... +%! -431.196 637.253 11.085 163.791]; ... +%! [-75.3626 -253.2337 457.1678 328.5714; ... +%! 438.7659 -653.6278 -7.9953 380.9336]; ... +%! [-89.5841 344.9716 -275.3876 457.1429; ... +%! -170.3613 237.8858 1.0469 158.0765];... +%! [32.900 -298.704 145.804 437.143; ... +%! -243.903 369.597 -34.265 226.648]; ... +%! [-99.081 409.127 -352.903 317.143; ... +%! 55.289 -114.223 -26.781 318.076]; ... +%! [-342.231 191.266 168.108 274.286; ... +%! 58.870 -38.083 -89.358 232.362]}; + +%!test % x-Reflection +%! v = shapecentroid (shape)(:); +%! T = createLineReflection([0 0 1 0]); +%! nshape = shapetransform (shape, T); +%! vn = shapecentroid (nshape)(:); +%! assert(vn,T(1:2,1:2)*v); + +%!test % Rotation +%! v = shapecentroid (shape)(:); +%! T = createRotation(v.',pi/2); +%! nshape = shapetransform (shape, T); +%! vn = shapecentroid (nshape)(:); +%! assert(vn,v,1e-2); + +%!test % Translation +%! v = shapecentroid (shape)(:); +%! nshape = shapetransform (shape, -v); +%! vn = shapecentroid (nshape)(:); +%! assert(vn,[0; 0],1e-2);