I am implementing ambient occlusion using a cosine weighted hemisphere. So far I assume that the hemisphere is centered about normal (0,0,1), which is pointing in the positive z direction. I cast rays about this hemisphere using zenith and azimuth angles $$0 < \theta < \pi/2, 0 < \phi < 2\pi$$, where the direction of the ray is computed as (sin theta cos phi, sin theta sin phi, cos theta).
How could I generalize this to work for an arbitrary normal?
Find a local coordinate system of the surface.
I.e. let the z-axis be the normal. Then, find two orthogonal axes for x and y. This can be achieved by first assuming (0, 1, 0) as the y-axis (or an alternate axis if this is the normal). Then calculate the x-axis as x = cross(y, z) and refine the y-axis as y = cross(z, x).
Then you can use your point calculation (p) as weights on this coordinate system:
dir = p.x * x + p.y * y + p.z * z
