4
\$\begingroup\$

I am solving a problem of increasing the ROI of a retinal image using the below algorithm-

First, the set of pixels of the exterior border of the ROI is determined, i.e., pixels that are outside the ROI and are neighbors (using four-neighbourhood) to pixels inside it. Then, each pixel value of this set is replaced with the mean value of its neighbors (this time using eight-neighbourhood) inside the ROI. Finally, the ROI is expanded by the inclusion of this altered set of pixels. This process is repeated and can be seen as artificially increasing the ROI.

My image is of following type- enter image description here

and the mask on which I am applying the algorithm is below- enter image description here

Pseudocode for the above approach which I got from here is below-

while there are pixels not in the ROI:
    border_pixels = []

    # find the border pixels
    for each pixel p=(i, j) in image
        if p is not in ROI and ((i+1, j) in ROI or (i-1, j) in ROI or (i, j+1) in ROI or (i, j-1) in ROI)):
            add p to border_pixels

    # calculate the averages
    for each pixel p in border_pixels:
        color_sum = (0, 0, 0)
        count = 0
        for each pixel n in 8-neighborhood of p:
            if n in ROI:
                color_sum += color(n)
                count += 1
        color(p) = color_sum / count

    # update the ROI
    for each pixel p=(i, j) in border_pixels:
        set p to be in ROI

and my implementation for the above pseudocode is below-

img = Image.open(path_dir)
pixelMap = img.load()
@jit
def roifun(img,pixelMap):
    roi = []
    for i in range(img.size[0]):
        for j in range(img.size[1]):
            if pixelMap[i,j] == 255:
                roi.append([i,j])

    return roi

roi= roifun(img,pixelMap)
notroi = img.size[0]*img.size[1] - len(roi)


@jit
def border_enhance(img,pixelMap,roi,notroi):
    while(notroi):
        border_pixels = []
        for i in range(img.size[0]):
            for j in range(img.size[1]):
                if [i,j] not in roi and ([i+1, j] in roi or [i-1, j] in roi or [i, j+1] in roi or [i, j-1] in roi):
                    border_pixels.append([i,j])

        for (each_i,each_j) in border_pixels:
            color_sum = 0
            count = 1
            eight_neighbourhood = [[each_i-1,each_j],[each_i+1,each_j],[each_i,each_j-1],[each_i,each_j+1],[each_i-1,each_j-1],[each_i-1,each_j+1],[each_i+1,each_j-1],[each_i+1,each_j+1]]
            for pix_i,pix_j in eight_neighbourhood:
                if (pix_i,pix_j) in roi:
                    color_sum+=pixelMap[pix_i,pix_j]
                    count+=1
            pixelMap[each_i,each_j]=(color_sum//count)

        for (each_i,each_j) in border_pixels:
            roi.append([each_i,each_j])
            border_pixels.remove([each_i,each_j])
            notroi = notroi-1
            print(notroi)

border_enhance(img,pixelMap,roi,notroi)

I run this code for image of dimension 50×50 and it is running correctly but for an image of larger size like 512×512, it is taking too long a time. I also tried modifying it using numba but it gave me plenty of warnings and then taking the same time as before.

\$\endgroup\$
  • \$\begingroup\$ Looping over pixels in Python will make your code very slow. Try Numba to compile it, or vectorize your code. \$\endgroup\$ – Cris Luengo Oct 4 at 2:18
  • \$\begingroup\$ @CrisLuengo I tried running it using Numba but it is taking the same time after giving warnings. \$\endgroup\$ – Mark Oct 11 at 18:20
  • \$\begingroup\$ I'm not sure I understand the code— you defined a border_enhance function, but it never gets called. \$\endgroup\$ – 200_success Oct 11 at 22:29
  • \$\begingroup\$ I just tried to post a minimal example, I edited it with by calling it now \$\endgroup\$ – Mark Oct 11 at 22:34
  • \$\begingroup\$ Okay posted it in the right place now! Let me know if you have any questions with the code. \$\endgroup\$ – CodeSurgeon Oct 21 at 0:42
2
\$\begingroup\$

The secret to performance is to choose the appropriate data structures which lay out the memory as raw arrays of bytes (or uint8_t values). This is done by creating cython typed memoryviews from numpy arrays or images loaded through PIL. This applies to the border list as well; that is, I create an array of (x, y) coordinates rather than a list of python tuples (which is slow, not contiguous in memory, and requires conversions to access the data from python to C or vice versa. Below is the preliminary code that I have so far, hope this helps.

from libc.stdint cimport *
from libc.string cimport *
cimport libc.math as c_math
import numpy as np
from PIL import Image

cpdef uint8_t[:, :, :] load_image(str image_path):
    cdef uint8_t[:, :, :] image_data
    image = Image.open(image_path).convert("RGBA")
    image_np = np.array(image)
    image_np.setflags(write=1)
    image_data = image_np
    return image_data

cpdef save_image(uint8_t[:, :, :] image_data, str image_path):
    image = Image.fromarray(np.array(image_data)).convert("RGBA")
    image.save(image_path)

cpdef border_enhance(uint8_t[:, :, :] image, uint8_t[:, :, :] mask):
    cdef:
        uint8_t[:, :, :] out
        uint8_t[:] color
        uint8_t[:, :] roi
        uint32_t[:, :] border
        size_t i, j, k
        size_t w = image.shape[0]
        size_t h = image.shape[1]
        size_t num_border
        size_t num_not_roi = 0
        float avg[4]
        size_t avg_count
        size_t x, y
        size_t sx, sy
        int a, b

    roi = np.zeros((w, h), dtype=np.uint8)
    border = np.zeros((w * h, 2), dtype=np.uint32)#assumes no product overflow
    out = image[:, :, :]

    #define ROI from mask; only white pixels are part of the mask?
    for i in range(w):
        for j in range(h):
            color = mask[i, j]
            if color[0] == 255 and color[1] == 255 and color[2] == 255 and color[3] == 255:
                roi[i, j] = True
            else:
                roi[i, j] = False
                num_not_roi += 1

    while num_not_roi:

        #Create border
        k = 0
        for i in range(w):
            for j in range(h):
                if not roi[i, j]:
                    #assumes edges wrap over
                    if roi[i-1, j] or roi[i+1, j] or roi[i, j-1] or roi[i, j+1]:
                        border[k, 0] = i
                        border[k, 1] = j
                        k += 1

        for i in range(k):
            x = border[i, 0]
            y = border[i, 1]
            avg_count = 0
            avg = [0, 0, 0, 0]
            for a in range(-1, 2):
                for b in range(-1, 2):
                    sx = x+a
                    sy = y+b
                    if roi[sx, sy]:
                        #print(np.array(image[sx, sy]))
                        avg[0] += image[sx, sy, 0]
                        avg[1] += image[sx, sy, 1]
                        avg[2] += image[sx, sy, 2]
                        avg[3] += image[sx, sy, 3]
                        avg_count += 1
            for a in range(4):
                avg[a] /= avg_count
                out[x, y, a] = <uint8_t>c_math.round(avg[a])

        for i in range(k):
            x = border[i, 0]
            y = border[i, 1]
            roi[x, y] = True
            num_not_roi -= 1
            #dprint(num_not_roi)

    save_image(out, "./images/out.png")

cdef:
    uint8_t[:, :, :] image
    uint8_t[:, :, :] mask
    size_t count

image = load_image("./images/image.png")
mask = load_image("./images/mask.png")
border_enhance(image, mask)

That being said, the results do not seem particularly correct. The intermediate border step generates the correct pixel results, but the final result is wrong. Here is my output for the out.png after the call to the border_enhance function:

enter image description here

\$\endgroup\$
  • 1
    \$\begingroup\$ damn, you wrote it in Cython. You are one badass programmer.I'll check it where the problem lies and will get back at you, though writing in cython is super dope. \$\endgroup\$ – Mark Oct 21 at 1:14
  • 1
    \$\begingroup\$ @Mark It is not too hard once you are familiar with how things are laid out in memory. I would take a look at the cython docs for how to use memoryviews. I would also look at this to get an idea as to what makes python so slow and what to avoid. Writing fast cython code is basically like writing C code with a nicer syntactic sugar, so taking some time to learn some C basics (e.g. for memory management and data structures) will pay off in the long run. \$\endgroup\$ – CodeSurgeon Oct 21 at 1:34
  • \$\begingroup\$ Thanks a lot for the links. I was finding it hard to wrap my head around cython, hope they will help me out. Thanks. \$\endgroup\$ – Mark Oct 21 at 15:08
  • 1
    \$\begingroup\$ @Mark Out of curiosity, were you able to figure out what was wrong with the code or get your issue solved? \$\endgroup\$ – CodeSurgeon Nov 7 at 1:41
  • \$\begingroup\$ There are some unusual edges which are present in the boundary which are responsible for the above unexpected increase in images, I got it from here-stackoverflow.com/questions/58512248/… but the thing is it is still slow, taking 149 sec for a single image so I am trying to convert the above solution into cython using your approach \$\endgroup\$ – Mark Nov 7 at 14:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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