3

Hello I have an input image that has been imported as grayscale in float values. I am trying to use for loops to find the intensity this image. Using two for loops to loop over each row and column is not idea for larger image.

image = cv2.imread("W_A03.jpg", cv2.IMREAD_GRAYSCALE).astype(float)/255.0

intimage = np.zeros(image.shape, dtype=image.dtype)
output = np.zeros(image.shape, dtype=image.dtype)
w_Image, h_Image = image.shape[1], image.shape[0]

for h_Idx in range(h_Image):
        for w_Idx in range(w_Image):
            intimage[h_Idx, w_Idx] = image[0:h_Idx, 0:w_Idx].sum()

The resulting intimage should looks something like this:

intensity

Is there a shorter, cleaner, or more pythonic way to perform such tasks?

Edit: I am trying to use adaptive thresholding on a grayscale image. The pseudo code can be found here procedure AdaptiveThreshold

Kana
  • 43
  • 1
  • 5
  • Hi, interesting question. Could you please provide a full [MWE](https://stackoverflow.com/help/minimal-reproducible-example)? – mapf Oct 22 '21 at 21:54
  • 2
    Instead of `np.zeros(image.shape, dtype=image.dtype)`, use `np.zeros_like(image)` – Mad Physicist Oct 22 '21 at 21:55
  • Oh, Hi. I am following the pseudo code for adaptive thresholding like the one in https://stackoverflow.com/a/29599155/8520575 – Kana Oct 22 '21 at 22:05
  • @mapf. Posted a whole new answer. I think this makes more sense, since I actually read the question this time. – Mad Physicist Oct 22 '21 at 22:45

2 Answers2

5

You appear to be looking for a 2D cumulative sum. As it happens, np.cumsum is separable:

intimage = image.cumsum(0).cumsum(1)

Alternatively,

intimage = image.cumsum(1).cumsum(0)

The cumulative sum at a given index is the sum of the elements up to that point, which is exactly what your current loop is computing. This algorithm is O(n) with respect to the number of pixels, while yours is O(n2) because it recomputes from scratch for each pixel instead of reusing the prior computations.

If I were to rewrite your loop with the more efficient approach, it would look like this:

image = cv2.imread("W_A03.jpg", cv2.IMREAD_GRAYSCALE).astype(float)/255.0
intimage = np.zeros_like(image)
w_Image, h_Image = image.shape  # Use tuple unpacking to your advantage

# Fill in the first row
intimage[0, 0] = image0[0, 0]
for w_Idx in range(1, w_Image):
    intimage[0, w_Idx] = intimage[0, w_Idx - 1] + image[0, w_Idx]

for h_Idx in range(1, h_Image):
    # Fill in the first column
    intimage[h_Idx, 0] = image[h_Idx, 0]
    for w_Idx in range(1, w_Image):
        intimage[h_Idx, w_Idx] = intimage[h_Idx, w_Idx - 1] + intimage[h_Idx - 1, w_Idx] + image[h_Idx, w_Idx]
Mad Physicist
  • 107,652
  • 25
  • 181
  • 264
  • This is probably the most practical solution in using an optimized library function and is easy to parallelize. The problem more generally is the prefix sum. An O(n) algorithm can be found here:https://www.geeksforgeeks.org/prefix-sum-2d-array/. Based on observation that pixel intimage[I,j] only needs information from four cells in its neighborhood intimage[i-1:i,j-l:j] which would allow for reusing computations. The implementation in the link has for loops which could probably only be broken up into two orthogonal streams that need to be computed sequentially. – Mike Holcomb Oct 22 '21 at 22:57
  • 1
    @MikeHolcomb. I just posted a proper loop implementation that reflects your comment. That being said, `cumsum` is much easier to parallelize and not much slower. – Mad Physicist Oct 22 '21 at 22:59
  • Nice solution! Good that you went back and actually read the question this time :D – mapf Oct 23 '21 at 08:41
0

it has been a while since I posted this question. I ended up using cv2.integral2 to find the integral image and changed the function to accept the integral image as an input like the following:

def adaptiveThreshold(image, threshold, intimage):
    output = np.zeros(image.shape, dtype=image.dtype)
    # the image is in grayscale
    w_Image, h_Image = image.shape[1], image.shape[0]
    sample_window_1 = w_Image/12
    sample_window_2 = sample_window_1/2
    for h_Idx in range(h_Image):
        y0 = round(max(h_Idx - sample_window_2, 0))
        y1 = round(min(h_Idx + sample_window_2, h_Image - 1))
        for w_Idx in range(w_Image):
            x0 = round(max(w_Idx - sample_window_2, 0))
            x1 = round(min(w_Idx + sample_window_2, w_Image - 1))
            count = (y1-y0)*(x1-x0)
            sum_ = intimage[y1, x1]-intimage[y0, x1]-intimage[y1, x0] + intimage[y0, x0]
            if np.all(image[h_Idx, w_Idx]*count <= sum_*(100.-threshold)/100.):
                output[h_Idx, w_Idx] = 0.
            else:
                output[h_Idx, w_Idx] = 255.
   return output

read in the image file in grayscale and as type float

image = cv2.imread("file_name.jpg", cv2.IMREAD_GRAYSCALE).astype(float)/255.0

add one more row and column to the image, and then reshape

sum_in = np.zeros((image.shape[0]+1, image.shape[1]+1), image.dtype)
cv2.integral2(image, sum_in, -1)
sum_in = cv2.resize(sum_in, (image.shape[1], image.shape[0]))

define the mean as the threshold value

thresh_val = int(sum([sum(i) for i in image])*100/
                 sum(len(row) for row in image))

use the function defined above

output = adpativeThreshold(image, thresh_val, sum_in)
Kana
  • 43
  • 1
  • 5