function xb = rbfbnd(x1, x2, x, w)
%This is called by rbfbound.
%
%It calculate the boundary point crossing with a line

K = 4;

if nargin < 3
	x = [1     0     0    -1     0     0    -2; ...
 	     0     1    -1     0     2    -2     0];
	y = [-1   -1     -1     1     1     1     1];
	w = rbfweight(x, y);
end

y1 = rbf(x1, x, w);
y2 = rbf(x2, x, w);
if y1 * y2 >= 0
	disp('line not cross');
end;

while 1,
	d = x2 - x1;
	if norm(d) < 1e-5
		break;
	end;
	
	xt = repmat(x1, 1, K-1) + d * ((1:K-1)/K);
	yt = rbf(xt, x, w);
	if y1 < 0, 
		idx = find(yt > 0);
	else
		idx = find(yt < 0);
	end
	if isempty(idx),
		x1 = xt(:,K-1);
	elseif idx(1) > 1,
		x1 = xt(:,idx(1)-1);
		x2 = xt(:,idx(1));
	else
		x2 = xt(:, 1);
	end
end;
xb = (x1 + x2) / 2;