function s = h3d_scan_sample_points(s, ndes)
	% 
	% assume s has already been passed to h3d_scan_convert_points
	
	nrays = numel(s.thetas_deg);
	nscans = numel(s.tilt_deg);
	
	s.area = zeros(nscans,nrays)
	
	for k=2:nscans-1
	for i=2:nrays-1
		rho = s.readings(k,i);
		if not(rho > 0) 
			s.area(k,i) = 0;
			continue
		end
		
		if not(s.valid(k+1,i) && s.valid(k-1,i))
			side1 = 0;
		else
			v1 =  s.p3d{k+1,i} - s.p3d{k-1,i};
			side1 = norm(v1);
		end
		
		if not(s.valid(k,i+1) && s.valid(k,i-1))
			side2 = 0;
		else
			v2 = s.p3d{k,i+1} - s.p3d{k,i-1};
			side2 = norm( v2); 
		end
		
		
		s.area(k,i) = side1 * side2;
		
	end
	end
	
	tot_area = sum(sum(s.area));
	
	indexes = randperm(nscans*nrays);
%	indexes = 1:nscans*nrays;
	
	prob = zeros(1,nscans*nrays);
	for a=1:numel(indexes)
		i = mod(indexes(a),nrays) + 1;
		k = floor((indexes(a)-1)/nrays) + 1;
		prob(a) = s.area(k,i) / tot_area;
	end

	
	chosen = zeros(1,numel(indexes));
	tot = 0;
	for a=1:numel(indexes)
		tot = tot + prob(a);
		if tot > 1/ndes
			chosen(a) = 1;
			tot = 0;
		end
	end
	

	fprintf('Nchosen: %d\n',sum(chosen));

	s.sampled = zeros(nscans,nrays);
	s.sampled_points = zeros(3, sum(chosen));
	ns = 1;
	
	for a=1:numel(indexes)
		if chosen(a)
			i = mod(indexes(a),nrays) + 1;
			k = floor((indexes(a)-1)/nrays) + 1;
	
			s.sampled(k,i) = 1;
			s.sampled_points(1:3,ns) = s.p3d{k,i};
			ns = ns + 1;
		end
	end