1 ## Copyright (C) 2002 Etienne Grossmann <etienne@egdn.net>
3 ## This program is free software; you can redistribute it and/or modify it under
4 ## the terms of the GNU General Public License as published by the Free Software
5 ## Foundation; either version 3 of the License, or (at your option) any later
8 ## This program is distributed in the hope that it will be useful, but WITHOUT
9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 ## You should have received a copy of the GNU General Public License along with
14 ## this program; if not, see <http://www.gnu.org/licenses/>.
16 ## y = bound_convex(d,h,x,pad=0)
18 ## y : 3xQ : Corners that define the convex hull of the projection of x
19 ## in the plane d*y == v. The corners are sorted.
22 ## Author: Etienne Grossmann <etienne@egdn.net>
24 function y = bound_convex(d,h,x,pad)
26 if nargin<4, pad = 0 ; end
33 ## I don't really care if I'm not given coplanar points
35 # if prudent && any(abs(d*x-h)>10*eps),
36 # printf("bound_convex : Points are not on plane (max dist=%8.4g\n",...
41 ## Project x to plane { y | d'*y = h }
50 ## x = proplan (x,d,h);
52 c = mean(x')'*ones(1,P);
55 ext = zeros(1,P); # Extremal points
58 px = nuld'*xc; # Project on 2D
59 [chx,ich] = chull (px); # Find 2D convex hull
64 # [dum,jj] = max( xc(:,i)'*xc ) ;
66 # [dum,jj] = min( xc(:,i)'*xc ) ;
72 norms = sqrt( sum( y.^2 ) );
75 printf("bound_convex : Points project to line\n") ;
76 if sum( norms != 0 )!=2,
77 printf("bound_convex : Moreover the segment has more than 2 tips!!\n") ;
79 y = y(:,find(norms != 0)) ;
80 norms = norms(find(norms != 0));
83 ## Sort points so that they turn monotonously around the origin
84 d1 = y(:,1)'/norms(1) ;
85 if abs( d1*d1' - 1 )>10*eps,
86 error ("bound_convex : d1 ain't unit!\n");
89 norms = [1;1;1]*norms;
91 [dum,id2] = min( abs( d1*y./norms(1,:) ) ) ;
92 ## d2 = cross( d1, y(:,id2)' ) ;
96 [dum,iy] = sort( atan2( d2*y./norms(1,:), d1*y./norms(1,:) ) ) ;
102 ## foo = d2*y./norms(1,iy);
103 ## bar = d1*y./norms(1,iy);
108 ## Shift back y into place
109 y = y+c(:,1:size(y,2)) ;