2
\$\begingroup\$

As part of an assignment, I need to write a program which loads a colored image to the screen, converts it to HSI and displays each component in separate grayscale image.

My program:

#!/usr/bin/env python3

import cv2
import numpy as np
import math
from PIL import Image
from io import BytesIO
from matplotlib import pyplot as plt


IMAGE_LOCATION = "lena.png"

class RGB2HSI_Coverter:
    def __init__(self, img):
        self.img = img

    # Method calculates the intensity
    # I = (R + G + B) / 3
    @staticmethod
    def calc_intensity(split):
        (R, G, B) = split
        return np.divide(R + G + B, 3)

    # Method calculates the saturation
    # S = 1 - (3 * min[R,G,B]) / (R + G + B)
    @staticmethod
    def calc_saturation(split):
        (R, G, B) = split
        min_value = np.minimum(np.minimum(R, G), B)
        return 1 - np.divide(3 * min_value, (R + G + B))

    # Method calculates the hue
    # H = theta if B <= G, otherwise 360 - theta
    # Where theta = arccos(( 0.5 * [(R - G) + (R - B)] ) / ([(R - G)^2 + (R - B)(G - B)]^0.5))
    @staticmethod
    def calc_hue(split):
        (R, G, B) = split
        result_hue = np.copy(R)

        for i in range(0, B.shape[0]):
            for j in range(0, B.shape[1]):
                # Calculate numerator = ( 0.5 * [(R - G) + (R - B)] )
                numerator = 0.5 * ((R[i][j] - G[i][j]) + (R[i][j] - B[i][j]))
                
                # Calculate denominator = ([(R - G)^2 + (R - B)(G - B)]^0.5)
                denominator = math.sqrt((R[i][j] - G[i][j])**2 + ((R[i][j] - B[i][j]) * (G[i][j] - B[i][j])))
                
                # Calculate divistion = numerator / denominator
                divistion = np.divide(numerator, denominator)
                
                # Calculate theta = arccos(divistion)
                result_hue[i][j] = math.acos(divistion)

                # If B > G then H = 360 - theta
                if B[i][j] > G[i][j]:
                    result_hue[i][j] = ((360 * math.pi) / 180.0) - result_hue[i][j]

        return result_hue
        
    def convert_R2G_to_HSI(self):
        with np.errstate(divide='ignore', invalid='ignore'):
            # Split into channels
            block = np.float32(self.img) / 255
            split = (block[:,:,2], block[:,:,1], block[:,:,0])
            
            # Calculate saturation
            saturation = self.calc_saturation(split) 
            
            # Calculate intensity
            intensity = self.calc_intensity(split)
            
            # Calculate hue
            hue = self.calc_hue(split)
            
            # Merge channels into one image
            hsi_img = cv2.merge((hue, saturation, intensity))

            return hsi_img

def main():
    # Load RGB image
    img = Image.open(IMAGE_LOCATION).convert("RGB")
    
    # Display RGB image
    plt.imshow(img)
    plt.show()

    # Convert RGB image to HSI
    converter = RGB2HSI_Coverter(img)
    hsi = converter.convert_R2G_to_HSI()

    # Display HSI image
    plt.imshow(hsi)
    plt.show()

if __name__ == "__main__":
    main()

The methods that are shown in my textbook and that I need to use to convert to HSI are:

enter image description here

What do you think about my code? I'm more interested if my code does what expected to be done and if I handled all of the cases (less interested in code design, but still will be nice). Just looking for additional pair of eyes, to check the correctness of my code.

\$\endgroup\$
4
  • \$\begingroup\$ Seems right, except that you compute H in radian instead of degrees. Usually degrees are used so that it can be stuck in an integer array. Also, you can use array operations to compute H, which will be a lot faster. Finally, you should reuse the I result to compute S, no need to repeat the same computations twice. (Edit: Sorry, didn’t realize I was on Code Review, I shouldn’t have posted this as a comment. I don’t have time to write a full answer right now.) \$\endgroup\$ Nov 6, 2021 at 18:41
  • \$\begingroup\$ @CrisLuengo Using degrees rather than radians just because they can be represented as integers assumes that you don't care about the (big) hit to resolution. This assumption may or may not hold. \$\endgroup\$
    – Reinderien
    Nov 6, 2021 at 20:27
  • \$\begingroup\$ @Reinderien I didn’t say OP should round. I said that the common definition uses degrees, and I explained why this is the common definition. OpenCV even divides the degrees by two so that they fit in an 8-bit integer. It’s awful, but it’s also surprisingly hard to convince people that an image doesn’t need to be represented by 8-bit integers. I much prefer using single floats, but I’m in a minority. \$\endgroup\$ Nov 6, 2021 at 22:08
  • \$\begingroup\$ Why did you put everything in the class RGB2HSI_Coverter? I don't think that class needs to exist. IMO, it would be cleaner if it was living directly in the module. Then calc_hue and the others would be pure functions, taking the image as an input. \$\endgroup\$ Nov 7, 2021 at 16:49

1 Answer 1

1
\$\begingroup\$

It would help future users to either add a hint comment that cv2 comes from the pip package opencv-contrib-python-headless, or better yet write an appropriate requirements.txt.

We need to lose Lena. Particularly when there are so many other excellent test pictures around.

Your class structure is a little odd, betrayed partly by there being so many static methods. You can store an intermediate rgb array during construction, and then just list off a bunch of @propertys, one per channel conversion you're interested in. This will convert all of your statics.

In many cases, it's not actually useful to split your RGB array to three channels. Many of your expressions are not properly vectorised, and after vectorisation, the entire array can be operated on. Your hue calculation is especially amenable to vectorisation improvements, since both of the loops should go away, there should be no copy, and no calls to math.

divistion isn't a word; did you mean "division" or "quotient"?

Do not call show() twice. The nicer display is to show the original and converted images side-by-side.

Importantly,

displays each component in separate grayscale image

is not what your code does. Your code was interpreting the HSI space as RGB space on display. The suggested code below technically does not use grayscale, as the default viridis colour scale is better for most purposes, including IMO this one.

Suggested

#!/usr/bin/env python3
import numpy as np
from PIL import Image
from cv2.cv2 import merge  # pip: opencv-contrib-python-headless
from matplotlib import pyplot as plt
import matplotlib.image

# https://www.gannett-cdn.com/-mm-/6830b11dd7334714ee859c4104be6bcc69f3abec/c=0-323-2833-1924/local/-/media/2017/03/21/USATODAY/USATODAY/636257089475691215-USP-ENTERTAINMENT-WEIRD-AL-YANKOVIC-82386343.JPG?width=2833&height=1601&fit=crop&format=pjpg&auto=webp
from mpl_toolkits.axes_grid1 import make_axes_locatable
from mpl_toolkits.axes_grid1.axes_divider import AxesDivider

IMAGE_LOCATION = "weirdal.webp"


class RGB2HSI_Converter:
    def __init__(self, image: Image) -> None:
        self.image = image

        # To a three-dimensional ndarray with values 0-1
        block = np.float32(self.image) / 255

        # From BGR to RGB, last axis to first
        self.rgb = np.moveaxis(block[:, :, ::-1], source=2, destination=0)

    @property
    def intensity(self) -> np.ndarray:
        # I = (R + G + B) / 3
        return np.sum(self.rgb, axis=0) / 3

    @property
    def saturation(self) -> np.ndarray:
        # S = 1 - (3 * min[R,G,B]) / (R + G + B)
        min_values = np.min(self.rgb, axis=0)
        sum_values = np.sum(self.rgb, axis=0)
        return 1 - 3 * min_values / sum_values

    @property
    def hue(self) -> np.ndarray:
        """
        H = theta if B <= G, otherwise 360 - theta
        Where theta = arccos(( 0.5 * [(R - G) + (R - B)] ) / ([(R - G)^2 + (R - B)(G - B)]^0.5))
        """
        R, G, B = self.rgb

        # Calculate numerator = ( 0.5 * [(R - G) + (R - B)] )
        numerator = 0.5 * (2*R - G - B)
        # Calculate denominator = ([(R - G)^2 + (R - B)(G - B)]^0.5)
        denominator = np.sqrt((R - G)**2 + (R - B)*(G - B))

        # Calculate theta = arccos(division)
        result_hue = np.arccos(numerator / denominator)

        # If B > G then H = 360 - theta
        np.putmask(result_hue, B > G, 2*np.pi - result_hue)

        return result_hue

    @property
    def hsi(self) -> Image:
        with np.errstate(divide='ignore', invalid='ignore'):
            return merge((self.hue, self.saturation, self.intensity))


def imshow_with_bar(fig: plt.Figure, ax: plt.Axes, data: np.ndarray) -> None:
    divider: AxesDivider = make_axes_locatable(ax)
    cax: plt.Axes = divider.append_axes(position='right', size='5%', pad=0.05)
    image: matplotlib.image.AxesImage = ax.imshow(data)
    fig.colorbar(mappable=image, cax=cax, orientation='vertical')


def main():
    img = Image.open(IMAGE_LOCATION).convert("RGB")
    converter = RGB2HSI_Converter(img)

    fig, (
        (top_left, top_right),
        (bottom_left, bottom_right),
    ) = plt.subplots(nrows=2, ncols=2)

    top_right.set_title('Hue')
    imshow_with_bar(fig, top_right, converter.hue)

    bottom_left.set_title('Saturation')
    imshow_with_bar(fig, bottom_left, converter.saturation)

    bottom_right.set_title('Intensity')
    imshow_with_bar(fig, bottom_right, converter.intensity)

    top_left.set_title('Original (RGB)')
    top_left.imshow(img)

    plt.show()


if __name__ == "__main__":
    main()

images by channel

\$\endgroup\$
3
  • \$\begingroup\$ Thank you! Part of my assignment, I was asked to use grayscale. Since I didn't use it myself, could you please explain how should it be added to your suggested code? \$\endgroup\$
    – vesii
    Nov 7, 2021 at 15:10
  • \$\begingroup\$ @vesii Use Greys \$\endgroup\$
    – Reinderien
    Nov 7, 2021 at 16:39
  • \$\begingroup\$ Does it mean I need to just use .convert('L')? \$\endgroup\$
    – vesii
    Nov 7, 2021 at 21:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.