voxpts

 voxpts method - voxels within a sphere in given space

 $Id$

0001 function [pts, vals] = voxpts(obj, sp)
0002 % voxpts method - voxels within a sphere in given space
0003 %
0004 % $Id$
0005 
0006 if nargin < 2
0007   error('Need object and space arguments');
0008 end
0009 sp = mars_space(sp);
0010 
0011 mc = obj.centre;
0012 mr = obj.radius;
0013 
0014 % find voxels that are within sphere
0015 vox = sqrt(sum(sp.mat(1:3,1:3)'.^2));
0016 vr = mr ./ vox;
0017 vc = sp.mat \ [mc(1:3) 1]';
0018 vc = vc(1:3)';
0019 blim = [max([1 1 1; ceil(vc-vr)]); min([sp.dim(1:3); floor(vc+vr)])];
0020 [R,C,P]=ndgrid(...
0021     blim(1,1):blim(2,1),blim(1,2):blim(2,2),blim(1,3):blim(2,3));
0022 vXYZ = [R(:)';C(:)';P(:)'];
0023 o = ones(1, size(vXYZ, 2));
0024 pts = find(sqrt(sum(((vXYZ-(vc'*o)) .* (vox'*o)).^2)) <= mr);
0025 pts = vXYZ(:,pts);
0026 vals = ones(1, size(pts, 2));