function s = h3d_scan_find_normals(s)
	% s = h3d_scan_find_normals(s)
	% assume s has already been passed to h3d_scan_convert_points
	
	nrays = numel(s.thetas_deg);
	nscans = numel(s.tilt_deg);
	
	s.alpha = cell(nscans, nrays);
	s.alpha_valid = zeros(nscans, nrays);
	s.area = zeros(nscans, nrays);
	
	for k=2:nscans-1
	for i=2:nrays-1
		rho = s.readings(k,i);

		if not(rho > 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
		
		if (side2>0) && (side1>0)
		
			n = cross(v1,v2);
			n = n / norm(n);
		
			if n'*s.p3d{k,i} > 0
				n = -n;
			end
		
			s.alpha{k,i} = n;
			s.alpha_valid(k,i) = 1;

			s.area(k,i) = side1*side2;
		end
	end
	end