bound_convex

(tmp_for_tar/vrml.doc/bound_convex.m)

Function Synopsis

y = bound_convex(d,v,x,pad) 

Help text

 y = bound_convex(d,v,x,pad=0) 

  y : 3xQ : Corners that define the convex hull of the projection of x
            in the plane d*y == v. The corners are sorted.

Listing of function file tmp_for_tar/vrml.doc/bound_convex.m

## y = bound_convex(d,v,x,pad=0) 
##
##  y : 3xQ : Corners that define the convex hull of the projection of x
##            in the plane d*y == v. The corners are sorted.
##

## Author:        Etienne Grossmann  <etienne@isr.ist.utl.pt>

function y = bound_convex(d,v,x,pad) 

if  nargin<4, pad = 0 ; end

d = d(:)' ;

P = size(x,2) ;
prudent = 1;

## I don't really care if I'm not given coplanar points
## 
#  if prudent && any(abs(d*x-v)>10*eps),
#    printf("bound_convex : Points are not on plane (max dist=%8.4g\n",...
#  	 max(abs(d*x-v)));
#    keyboard
#  end

x = proplan (x,d,v);

c = mean(x')'*ones(1,P);
xc = x-c ;

ext = zeros(1,P);		# Extremal points

nuld = null(d);
px = nuld'*xc; 			# Project on 2D 
[chx,ich] = chull (px);		# Find 2D convex hull
ext(ich) = 1;
## keyboard
#  for i = 1:P,

#    [dum,jj] = max( xc(:,i)'*xc ) ;
#    ext(jj) = 1 ;
#    [dum,jj] = min( xc(:,i)'*xc ) ;
#    ext(jj) = 1 ;
#  end

y = xc(:,find(ext)) ;

norms = sqrt( sum( y.^2 ) );

if any(norms==0),
  printf("bound_convex : Points project to line\n") ;
  if sum( norms != 0 )!=2,
    printf("bound_convex : Moreover the segment has more than 2 tips!!\n") ;
  end
  y = y(:,find(norms != 0)) ;
  norms = norms(find(norms != 0));
  return
end
## Sort points so that they turn monotonously around the origin
d1 = y(:,1)'/norms(1) ;
if abs( d1*d1' - 1 )>10*eps,
  printf("bound_convex : d1 ain't unit!\n");
  keyboard
end
norms = [1;1;1]*norms;

[dum,id2] = min( abs( d1*y./norms(1,:) ) ) ;
## d2 = cross( d1, y(:,id2)' ) ;
d2 = cross( d1, d ) ;
d2 = d2/norm(d2) ;

[dum,iy] = sort( atan2( d2*y./norms(1,:), d1*y./norms(1,:) ) ) ;

y = y(:,iy) ;



## foo = d2*y./norms(1,iy);
## bar = d1*y./norms(1,iy);
## plot(bar,foo) ;

## keyboard

## Shift back y into place 
y = y+c(:,1:size(y,2)) ;