make sample_point_in_half_sphere_shell more efficient

[soskek]

Thank you for your great work.

sample_point_in_half_sphere_shell function can be extremely slow if the region between inner_radius and outer_radius is narrow because many rejections happen. For example, if they are 0.99999 and 1.00001, it takes 0.44 sec/iter.

This PR implements a sampling function, which has the equivalent feature but is much faster (e.g., 0.000015 sec/iter in the same condition).
See https://math.stackexchange.com/questions/1111894/random-points-in-spherical-shell in detail of the algorithm.

METHOD   sec/iter 	                      inner_radius  outer_radius

CURRENT   6.361722946166992e-05 sec/iter for	 0.9 1.1
EFFICIENT 1.516876220703125e-05 sec/iter for	 0.9 1.1
CURRENT   0.0004531641721725464 sec/iter for	 0.99 1.01
EFFICIENT 1.4918088912963867e-05 sec/iter for	 0.99 1.01
CURRENT   0.004436958122253418 sec/iter for	 0.999 1.001
EFFICIENT 1.507430076599121e-05 sec/iter for	 0.999 1.001
CURRENT   0.04548759651184082 sec/iter for	 0.9999 1.0001
EFFICIENT 1.4777421951293946e-05 sec/iter for	 0.9999 1.0001
CURRENT   0.4282598567008972 sec/iter for	 0.99999 1.00001
EFFICIENT 1.5428066253662108e-05 sec/iter for	 0.99999 1.00001

---

profile snippet:

if __name__ == "__main__":
    import kubric as kb
    import numpy as np
    import time

    N = 10000
    def check_current_version(inner_radius, outer_radius):
        start = time.time()
        outs = []
        for frame in range(N):
            pos = kb.sample_point_in_half_sphere_shell(inner_radius, outer_radius, 0.0)
            outs.append(pos)
        print("CURRENT  ", (time.time() - start) / N, "sec/iter for\t", inner_radius, outer_radius)
        return np.array(outs)

    def check_efficient_version(inner_radius, outer_radius):
        start = time.time()
        outs = []
        for frame in range(N):
            pos = efficient__sample_point_in_half_sphere_shell(inner_radius, outer_radius, 0.0)
            outs.append(pos)
        print("EFFICIENT", (time.time() - start) / N, "sec/iter for\t", inner_radius, outer_radius)
        return np.array(outs)

    # check time
    check_current_version(inner_radius=0.9, outer_radius=1.1)
    check_efficient_version(inner_radius=0.9, outer_radius=1.1)
    check_current_version(inner_radius=0.99, outer_radius=1.01)
    check_efficient_version(inner_radius=0.99, outer_radius=1.01)
    check_current_version(inner_radius=0.999, outer_radius=1.001)
    check_efficient_version(inner_radius=0.999, outer_radius=1.001)
    check_current_version(inner_radius=0.9999, outer_radius=1.0001)
    check_efficient_version(inner_radius=0.9999, outer_radius=1.0001)
    check_current_version(inner_radius=0.99999, outer_radius=1.00001)
    check_efficient_version(inner_radius=0.99999, outer_radius=1.00001)

[soskek]

Visual comparison of samples from the current algorithm (red) and the new one (blue); they are matched.

inner_radius=0.5, outer_radius=1.5

[Qwlouse]

That looks great! Thanks for the upgrade.