function s = h3d_compute_cube_ht(p, s, points)

	n = size(points,2);

	for face=1:6
	for i=1:s.cube_ncells
		progress( (face-1)*s.cube_ncells + i, 6*s.cube_ncells);
		
	for j=1:s.cube_ncells
		normal = s.cube_normals{face,i,j};
	
		rhos = normal' * points;
		
		for k=1:n
		%	if p.ht_points_are_from_scan
		%		if rhos(k) <= 0
		%			continue
		%		end
		%	end
			
			[r,alpha] = h3d_rho2index(s, rhos(k));
			
			if r > 0
				s.cube_ht(face,i,j,r) = s.cube_ht(face,i,j,r) + 1-abs(alpha);
				
				if (alpha > 0) && (r < s.rho_ncells)
					s.cube_ht(face,i,j,r+1) = s.cube_ht(face,i,j,r+1) + abs(alpha);
				end
				
				if (alpha < 0) && (r > 1)
					s.cube_ht(face,i,j,r-1) = s.cube_ht(face,i,j,r-1) + abs(alpha);
				end
				
			end
		end
	
	end
	end
	end
		
	s.points = [s.points points];


	% compute spectrum

	for face=1:6
	for i=1:s.cube_ncells
		
	for j=1:s.cube_ncells

		normal = s.cube_normals{face,i,j};
	
	
		rhos = squeeze( s.cube_ht(face, i, j, :) );
		
		s.cube_hs(face, i, j) = sum(rhos.^2);
	end
	end
	end