function a = createUniformWithinBoundary(npoints, areaParams)
% Random points, guaranteed to be within a user-defined geometrical
% volume. (Not very well thought-of.. just needed this kind of thing
% the other day).
%
%
% The numbers are drawn from a uniform random distribution, and filtered to
% be 
%   (1) bounded by a bounding box,
%       given as two vectors that contain lower and upper limits of 
%       coordinate values, 
%
% and within the union of 
%
%   (2) any user-given hyperspheres and 
%
%   (3) any user-given "hypercylinders".
%

tbd = repmat(true, [1,npoints]);

hyperc = areaParams.hyperc;   % bounding hypercube
center = areaParams.center;   % centers of bounding balls
radius = areaParams.radius;   % radii of bounding balls

tubeA = areaParams.tubeA;     % beginnings of 'cylinder' centers
tubeB = areaParams.tubeB;     % endings of 'cylinder' centers
tubeR = areaParams.tubeR;     % radii of 'cylinders'

dim = size(hyperc,1);

while(true)
    % Re-draw those who weren't yet inside limits:
    a(:,tbd) = bsxfun(@times, hyperc(:,2)-hyperc(:,1), rand(dim,sum(tbd)));
    a(:,tbd) = bsxfun(@plus, hyperc(:,1), a(:,tbd));

    % Filter with geometric area:
    ok = repmat(false,[1 sum(tbd)]);
    % Hyperballs
    for ii = 1:size(center,2)
      c = center(:,ii);
      r = sum( (a(:,tbd) - repmat(c, [1, sum(tbd)] ) ) .^ 2, 1);
      ok = ok | (r < radius(ii).^2);
    end
    % "Hypertubes"
    for ii = 1:size(tubeA,2)
        tA = tubeA(:,ii);
        tB = tubeB(:,ii);
        tV = tB-tA; % Tube vector
        tL = sqrt(tV'*tV); % Tube length
        tD = tV ./ tL;  % Tube direction
        
        % We calculate (point-tubeA)
        trP = bsxfun(@minus, a(:,tbd), tubeA);
        alongDist = tD' * trP;
        crossDist = sqrt(sum(bsxfun(@minus, trP, bsxfun(@times, tD, alongDist)) .^2));
        
        ok = ok | ((0<alongDist) & (alongDist < tL)&(crossDist < tubeR(ii)));
    end
    
    tbd(tbd) = ~ok;
    %tbd(tbd) = ~tbd(tbd);
    if sum(tbd)==0, break; end;
end
