blinn-phong shading with numpy

Question

I am trying to implement blinn-phong shading in numpy for educational purposes. However I am stuck at debugging what parameters are doing for several days now.

My general idea was the following. Since the equation was given for a channel. I apply the model to each color channel to get the relative pixel intensities in the channel, then regroup the channels back togather to have all the image.

My lambertian coefficiant does not seem to take into account the light position changes, but it does change the pixel intensity but other parameters have almost no effect on the output. Any help would be appreciated. Here are the relative bits of the code (full code is here for anyone interested):

def normalize_1d_array(arr):
    "Normalize 1d array"
    assert arr.ndim == 1
    result = None
    if np.linalg.norm(arr) == 0:
        result = arr
    else:
        result = arr / np.linalg.norm(arr)
    return result


def normalize_3col_array(arr):
    "Normalize 3 column array"
    assert arr.shape[1] == 3
    assert arr.ndim == 2
    normal = np.copy(arr)
    normal[:, 0] = normalize_1d_array(normal[:, 0])
    normal[:, 1] = normalize_1d_array(normal[:, 1])
    normal[:, 2] = normalize_1d_array(normal[:, 2])
    return normal


def get_vector_dot(arr1, arr2):
    "Get vector dot product for 2 matrices"
    assert arr1.shape == arr2.shape
    newarr = np.sum(arr1 * arr2, axis=1, dtype=np.float32)
    return newarr


class LightSource:
    "Simple implementation of a light source"

    def __init__(self,
                 x=10.0,  # x
                 y=5.0,  # y
                 z=0.0,  # light source at infinity
                 intensity=1.0,  # I_p
                 ambient_intensity=1.0,  # I_a
                 ambient_coefficient=0.1,  # k_a
                 light_power=80.0):
        "light source"
        self.x = x
        self.y = y
        if z is not None:
            assert isinstance(z, float)
        self.z = z
        self.intensity = intensity
        self.power = light_power
        self.ambient_intensity = ambient_intensity  # I_a
        self.ambient_coefficient = ambient_coefficient  # k_a
        # k_a can be tuned if the material is known

    def copy(self):
        "copy self"
        return LightSource(x=self.x,
                           y=self.y,
                           z=self.z,
                           intensity=self.intensity,
                           light_power=self.power)


class ChannelShader:
    "Shades channels"

    def __init__(self,
                 coordarr: np.ndarray,
                 light_source: LightSource,  # has I_a, I_p, k_a
                 surface_normal: np.ndarray,
                 imagesize: (int, int),
                 color: np.ndarray,  # they are assumed to be O_d and O_s
                 spec_coeff=0.5,  # k_s
                 screen_gamma=2.2,
                 diffuse_coeff=0.9,  # k_d
                 attenuation_c1=2.0,  # f_attr c1
                 attenuation_c2=0.0,  # f_attr c2 d_L coefficient
                 attenuation_c3=0.0,  # f_attr c3 d_L^2 coefficient
                 shininess=270.0  # n
                 ):
        self.light_source = light_source
        self.light_intensity = self.light_source.intensity  # I_p
        self.ambient_coefficient = self.light_source.ambient_coefficient  # k_a
        self.ambient_intensity = self.light_source.ambient_intensity  # I_a
        self.coordarr = coordarr
        self.surface_normal = np.copy(surface_normal)
        self.screen_gamma = screen_gamma
        self.shininess = shininess
        self.diffuse_coeff = diffuse_coeff  # k_d
        self.diffuse_color = normalize_1d_array(color)  # O_d: object diffuse color
        self.spec_color = normalize_1d_array(color)  # O_s: object specular color
        self.spec_coeff = spec_coeff  # k_s: specular coefficient
        self.imsize = imagesize
        self.att_c1 = attenuation_c1
        self.att_c2 = attenuation_c2
        self.att_c3 = attenuation_c3

    def copy(self):
        return ChannelShader(coordarr=np.copy(self.coordarr),
                             light_source=self.light_source.copy(),
                             surface_normal=np.copy(self.surface_normal),
                             color=np.copy(self.diffuse_coeff) * 255.0)

    @property
    def distance(self):
        yarr = self.coordarr[:, 0]  # row nb
        xarr = self.coordarr[:, 1]  # col nb
        xdist = (self.light_source.x - xarr)**2
        ydist = (self.light_source.y - yarr)**2
        return xdist + ydist

    @property
    def distance_factor(self):
        resx = self.imsize[1]
        factor = self.distance / self.light_source.z * resx
        return 1.0 - factor

    @property
    def light_direction(self):
        "get light direction matrix (-1, 3)"
        yarr = self.coordarr[:, 0]
        xarr = self.coordarr[:, 1]
        xdiff = self.light_source.x - xarr
        ydiff = self.light_source.y - yarr
        light_matrix = np.zeros((self.coordarr.shape[0], 3))
        light_matrix[:, 0] = ydiff
        light_matrix[:, 1] = xdiff
        light_matrix[:, 2] = self.light_source.z
        # light_matrix[:, 2] = 0.0
        return light_matrix

    @property
    def light_attenuation(self):
        """
        Implementing from Foley JD 1996, p. 726

        f_att : light source attenuation function:
        f_att = min(\frac{1}{c_1 + c_2{\times}d_L + c_3{\times}d^2_{L}} , 1)
        """
        second = self.att_c2 * self.distance
        third = self.att_c3 * self.distance * self.distance
        result = self.att_c1 + second + third
        result = 1 / result
        return np.where(result < 1, result, 1)

    @property
    def normalized_light_direction(self):
        "Light Direction matrix normalized"
        return normalize_3col_array(self.light_direction)

    @property
    def normalized_surface_normal(self):
        return normalize_3col_array(self.surface_normal)

    @property
    def costheta(self):
        "set costheta"
        # pdb.set_trace()
        costheta = get_vector_dot(
            arr1=self.normalized_light_direction,
            arr2=self.normalized_surface_normal)
        # products of vectors
        costheta = np.abs(costheta)  # as per (Foley J.D, et.al. 1996, p. 724)
        return costheta

    @property
    def ambient_term(self):
        "Get the ambient term I_a * k_a * O_d"
        term = self.ambient_coefficient * self.ambient_intensity
        return term * self.diffuse_color

    @property
    def view_direction(self):
        "Get view direction"
        # pdb.set_trace()
        cshape = self.coordarr.shape
        coord = np.zeros((cshape[0], 3))  # x, y, z
        coord[:, :2] = -self.coordarr
        coord[:, 2] = 0.0  # viewer at infinity
        coord = normalize_3col_array(coord)
        return coord

    @property
    def half_direction(self):
        "get half direction"
        # pdb.set_trace()
        arr = self.view_direction + self.normalized_light_direction
        return normalize_3col_array(arr)

    @property
    def spec_angle(self):
        "get spec angle"
        specAngle = get_vector_dot(
            arr1=self.half_direction,
            arr2=self.normalized_surface_normal)
        return np.where(specAngle > 0.0, specAngle, 0.0)

    @property
    def specular(self):
        return self.spec_angle ** self.shininess

    @property
    def channel_color_blinn_phong(self):
        """compute new channel color intensities
        Implements: Foley J.D. 1996 p. 730 - 731, variation on equation 16.15
        """
        second = 1.0  # added for structuring code in this fashion, makes
        # debugging easier
        # lambertian terms
        second *= self.diffuse_coeff  # k_d
        second *= self.costheta  # (N \cdot L)
        second *= self.light_intensity  # I_p
        # adding phong terms
        second *= self.light_attenuation  # f_attr
        second *= self.diffuse_color  # O_d
        third = 1.0
        third *= self.spec_color  # O_s
        third *= self.specular  # (N \cdot H)^n
        third *= self.spec_coeff  # k_s
        result = 0.0
        result += self.ambient_term  # I_a × k_a × O_d
        result += second
        result += third
        pdb.set_trace()
        return result

Thanks

Answer 1

Well after all, the implementation did not have a lot of problems, however the images I was working with required very weird values for parameters due to their specific condition in which they are produced.

Most images I have used contained rough surfaces with unglazed clay as material and the images were taken in a controlled environment with a single source of light, contrary to real world environments where objects are illuminated from multiple light spots.

So most of the parameters about ambient illumination and specular reflection did not make much sense in usage.

I am putting here the relative parts of my implementation as a reference for future users, be sure NOT to use default values.

Some details about the implementation:

• It largely follows the equation specified in Foley J.D. et.al., 1996, p. 730 - 731, no 16.15 with the variation of adding a halfway vector, required for blinn-phong.

• ChannelShader expects the following:

  • Coordinates of channel pixels of shape: (-1, 2)
  • Surface normals of shape: (-1, 3)
  • Channel colors of shape: (-1,)
  • A light source of the type LightSource
  • As stated above, please do change the default values before proceeding with experimentation.

You can use same surface normals for each channel, if you are shading multiple channels.

Final word of caution, it is noticably slow, even with numpy. The proper way to render shading is gpu based libraries like pyopengl etc. I have not tested it with gpu ports of numpy, like cupy, nor with other libraries like numba etc though:

def normalize_1d_array(arr):
    "Normalize 1d array"
    assert arr.ndim == 1
    result = None
    if np.linalg.norm(arr) == 0:
        result = arr
    else:
        result = arr / np.linalg.norm(arr)
    return result

def normalize_3col_array(arr):
    "Normalize 3 column array"
    assert arr.shape[1] == 3
    assert arr.ndim == 2
    normal = np.copy(arr)
    normal[:, 0] = normalize_1d_array(normal[:, 0])
    normal[:, 1] = normalize_1d_array(normal[:, 1])
    normal[:, 2] = normalize_1d_array(normal[:, 2])
    return normal


def get_vector_dot(arr1, arr2):
    "Get vector dot product for 2 matrices"
    assert arr1.shape == arr2.shape
    newarr = np.sum(arr1 * arr2, axis=1, dtype=np.float32)
    return newarr


class ImageArray:
    "Image array have some additional properties besides np.ndarray"

    def __init__(self, image: np.ndarray):
        assert isinstance(image, np.ndarray)
        self.image = image

    @property
    def norm_coordinates(self):
        "Get normalized coordinates of the image pixels"
        # pdb.set_trace()
        rownb, colnb = self.image.shape[0], self.image.shape[1]
        norm = np.empty_like(self.coordinates, dtype=np.float32)
        norm[:, 0] = self.coordinates[:, 0] / rownb
        norm[:, 1] = self.coordinates[:, 1] / colnb
        return norm

    @property
    def norm_image(self):
        "Get normalized image with pixel values divided by 255"
        return self.image / 255.0

    @property
    def coordinates(self):
        "Coordinates of the image pixels"
        rownb, colnb = self.image.shape[:2]
        coords = [[(row, col) for col in range(colnb)] for row in range(rownb)]
        coordarray = np.array(coords)
        return coordarray.reshape((-1, 2))

    @property
    def arrshape(self):
        "get array shape"
        return self.image.shape

    @property
    def flatarr(self):
        "get flattened array"
        return self.image.flatten()


def interpolateImage(imarr: ImageArray):
    "Interpolate image array"
    imshape = imarr.image.shape
    newimage = imarr.image.flatten()
    newimage = np.uint8(np.interp(newimage,
                                  (newimage.min(),
                                   newimage.max()),
                                  (0, 255))
                        )
    newimage = newimage.reshape(imshape)
    return ImageArray(newimage)


class LightSource:
    "Simple implementation of a light source"

    def __init__(self,
                 x=1.0,  # x
                 y=1.0,  # y
                 z=20.0,  # light source distance: 0 to make it at infinity
                 intensity=1.0,  # I_p
                 ambient_intensity=1.0,  # I_a
                 ambient_coefficient=0.000000002,  # k_a
                 ):
        "light source"
        self.x = x
        self.y = y
        if z is not None:
            assert isinstance(z, float)
        self.z = z
        self.intensity = intensity
        self.ambient_intensity = ambient_intensity  # I_a
        self.ambient_coefficient = ambient_coefficient  # k_a
        # k_a can be tuned if the material is known

    def copy(self):
        "copy self"
        return LightSource(x=self.x,
                           y=self.y,
                           z=self.z,
                           intensity=self.intensity,
                           light_power=self.power)


class ChannelShader:
    "Shades channels"

    def __init__(self,
                 coordarr: np.ndarray,
                 light_source: LightSource,  # has I_a, I_p, k_a
                 surface_normal: np.ndarray,
                 color: np.ndarray,  # they are assumed to be O_d and O_s
                 spec_coeff=0.1,  # k_s
                 spec_color=1.0,  # O_s: obj specular color. It can be
                 # optimized with respect to surface material
                 screen_gamma=2.2,
                 diffuse_coeff=0.008,  # k_d
                 # a good value is between 0.007 and 0.1
                 attenuation_c1=1.0,  # f_attr c1
                 attenuation_c2=0.0,  # f_attr c2 d_L coefficient
                 attenuation_c3=0.0,  # f_attr c3 d_L^2 coefficient
                 shininess=20.0  # n
                 ):
        self.light_source = light_source
        self.light_intensity = self.light_source.intensity  # I_p
        self.ambient_coefficient = self.light_source.ambient_coefficient  # k_a
        self.ambient_intensity = self.light_source.ambient_intensity  # I_a
        self.coordarr = coordarr
        self.surface_normal = np.copy(surface_normal)
        self.screen_gamma = screen_gamma
        self.shininess = shininess
        self.diffuse_coeff = diffuse_coeff  # k_d
        # self.diffuse_color = normalize_1d_array(color)  # O_d: obj diffuse color
        self.diffuse_color = color  # O_d: obj diffuse color
        self.spec_color = spec_color  # O_s
        self.spec_coeff = spec_coeff  # k_s: specular coefficient
        self.att_c1 = attenuation_c1
        self.att_c2 = attenuation_c2
        self.att_c3 = attenuation_c3

    def copy(self):
        return ChannelShader(coordarr=np.copy(self.coordarr),
                             light_source=self.light_source.copy(),
                             surface_normal=np.copy(self.surface_normal),
                             color=np.copy(self.diffuse_color))

    @property
    def distance(self):
        yarr = self.coordarr[:, 0]  # row nb
        xarr = self.coordarr[:, 1]  # col nb
        xdist = (self.light_source.x - xarr)**2
        ydist = (self.light_source.y - yarr)**2
        return xdist + ydist

    @property
    def light_direction(self):
        "get light direction matrix (-1, 3)"
        yarr = self.coordarr[:, 0]
        xarr = self.coordarr[:, 1]
        xdiff = self.light_source.x - xarr
        ydiff = self.light_source.y - yarr
        light_matrix = np.zeros((self.coordarr.shape[0], 3))
        light_matrix[:, 0] = ydiff
        light_matrix[:, 1] = xdiff
        light_matrix[:, 2] = self.light_source.z
        # light_matrix[:, 2] = 0.0
        # pdb.set_trace()
        return light_matrix

    @property
    def light_attenuation(self):
        """
        Implementing from Foley JD 1996, p. 726

        f_att : light source attenuation function:
        f_att = min(\frac{1}{c_1 + c_2{\times}d_L + c_3{\times}d^2_{L}} , 1)
        """
        second = self.att_c2 * self.distance
        third = self.att_c3 * self.distance * self.distance
        result = self.att_c1 + second + third
        result = 1 / result
        return np.where(result < 1, result, 1)

    @property
    def normalized_light_direction(self):
        "Light Direction matrix normalized"
        return normalize_3col_array(self.light_direction)

    @property
    def normalized_surface_normal(self):
        return normalize_3col_array(self.surface_normal)

    @property
    def costheta(self):
        "set costheta"
        # pdb.set_trace()
        costheta = get_vector_dot(
            arr1=self.normalized_light_direction,
            arr2=self.normalized_surface_normal)
        # products of vectors
        # costheta = np.abs(costheta)  # as per (Foley J.D, et.al. 1996, p. 724)
        costheta = np.where(costheta > 0, costheta, 0)
        return costheta

    @property
    def ambient_term(self):
        "Get the ambient term I_a * k_a * O_d"
        term = self.ambient_coefficient * self.ambient_intensity
        term *= self.diffuse_color
        # pdb.set_trace()
        return term

    @property
    def view_direction(self):
        "Get view direction"
        # pdb.set_trace()
        cshape = self.coordarr.shape
        coord = np.zeros((cshape[0], 3))  # x, y, z
        coord[:, :2] = -self.coordarr
        coord[:, 2] = 0.0  # viewer at infinity
        coord = normalize_3col_array(coord)
        return coord

    @property
    def half_direction(self):
        "get half direction"
        # pdb.set_trace()
        arr = self.view_direction + self.normalized_light_direction
        return normalize_3col_array(arr)

    @property
    def spec_angle(self):
        "get spec angle"
        specAngle = get_vector_dot(
            arr1=self.half_direction,
            arr2=self.normalized_surface_normal)
        return np.where(specAngle > 0.0, specAngle, 0.0)

    @property
    def specular(self):
        return self.spec_angle ** self.shininess

    @property
    def channel_color_blinn_phong(self):
        """compute new channel color intensities
        Implements: Foley J.D. 1996 p. 730 - 731, variation on equation 16.15
        """
        second = 1.0  # added for structuring code in this fashion, makes
        # debugging easier
        # lambertian terms
        second *= self.diffuse_coeff  # k_d
        second *= self.costheta  # (N \cdot L)
        second *= self.light_intensity  # I_p
        # adding phong terms
        second *= self.light_attenuation  # f_attr
        second *= self.diffuse_color  # O_d
        third = 1.0
        third *= self.spec_color  # O_s
        third *= self.specular  # (N \cdot H)^n
        third *= self.spec_coeff  # k_s
        result = 0.0
        #
        result += self.ambient_term  # I_a × k_a × O_d
        result += second
        result += third
        # pdb.set_trace()
        return result