Any idea to optimise this algorithm?

Question

Let's imagine I have a binary 40*40 matrix. In this matrix, values can be either ones or zeros.

I need to parse the whole matrix, for any value == 1, apply the following :

If the following condition is met, keep the value = 1, else modify the value back to 0:

Condition: in a square of N*N (centered on the currently evaluated value), I can count at least M values == 1.

N and M are parameters that can be set for the algorithms, naturally N<20 (in this case) and M<(N*N - 1).

The obvious algorithm is to iterate over the whole image and then each time a value == 1. Perform another search around this value. It would make an O^3 complex algorithm. Any idea to make this more efficient?

Edit: Some code to make this more easy to understand.

Let's create a randomly initialized 40*40 matrix of 1s and 0s. 5% (arbitrarily chosen) of the values are 1s, 95% are 0s.

import matplotlib.pyplot as plt
import numpy as np

%matplotlib inline

def display_array(image):
    image_display_ready = image * 255

    print(image_display_ready)

    plt.imshow(image_display_ready, cmap='gray')
    plt.show()

image = np.zeros([40,40])

for _ in range(80):  # 5% of the pixels are == 1
    i, j = np.random.randint(40, size=2)
    image[i, j] = 1

# Image displayed on left below before preprocessing    
display_array(image)   

def cleaning_algorithm_v1(src_image, FAT, LR, FLAG, verbose=False):
    """
    FAT  = 4 # False Alarm Threshold (number of surrounding 1s pixel values)
    LR   = 8 # Lookup Range +/- i/j value
    FLAG = 2 # 1s pixels kept as 1s after processing are flag with this value.
    """

    rslt_img = np.copy(src_image)

    for i in range(rslt_img.shape[0]):
        for j in  range(rslt_img.shape[1]):

            if rslt_img[i, j] >= 1:

                surrounding_abnormal_pixels = list()

                lower_i = max(i - LR, 0)
                upper_i = min(i + LR + 1, rslt_img.shape[0])

                lower_j = max(j - LR, 0)
                upper_j = min(j + LR + 1, rslt_img.shape[1])

                abnormal_pixel_count = 0

                for i_k in range(lower_i, upper_i):
                    for j_k in range(lower_j, upper_j):

                        if i_k == i and j_k == j:
                            continue

                        pixel_value = rslt_img[i_k, j_k]

                        if pixel_value >= 1:
                            abnormal_pixel_count += 1

                if abnormal_pixel_count >= FAT:
                    rslt_img[i, j] = FLAG

    rslt_img[rslt_img != FLAG] = 0
    rslt_img[rslt_img == FLAG] = 1

    return rslt_img

# Image displayed on right below after preprocessing  
display_array(cleaning_algorithm_v1(image, FAT=10, LR=6, FLAG=2))

Which gives the following:

Actually, the complexity of the algorithm you suggest is O((M^2)*(N^2)). You parse through a matrix of M^2 places and, whenever a value of 1 is found, you parse through another matrix of N^2 places. Worst case, you can find M^2 ones (MxM matrix is all ones). — RenanB, Jul 18 '18 at 02:33
Is the matrix 40x40? Time complexity only makes sense in the context of an asymptotic limit... — unutbu, Jul 18 '18 at 02:37
What is the condition when the location being evaluated is close to the edge? (Does the square get wrapped around or truncated or....?) — unutbu, Jul 18 '18 at 02:39
Well I don't know about reducing the complexity, but given that you know N and M at runtime, you could generate a matrix for each point (which would be 40^2) of size N^2 and compute the value for each matrix. ***However***, like @unutbu pointed out, you did not mention the corner cases. Also, what happens to the future points if you toggle one of the values? Do you use the old one or the new toggled one? If the old one, my suggestion of making it embarrassingly parallel works, otherwise it breaks. So you may want to state that in your problem question. — Haris Nadeem, Jul 18 '18 at 03:43
@haris-nadeem I use only the original matrix as source not the modified version — Jonathan DEKHTIAR, Jul 18 '18 at 12:50
@unutbu the objective is to try to optimize this algorithm ;) I just look for new ideas to make it faster ;) When the location is close to the edge, it is truncated. — Jonathan DEKHTIAR, Jul 18 '18 at 12:50

William Yolland · Accepted Answer · 2018-07-18T03:34:35.677

4

What about using convolution?

Your kernel would be an NxN window of 1's. In this case the kernel is separable so you can process the convolution as 2 1-D convolutions. You could do something like:

import numpy as np
from scipy.ndimage.filters import convolve1d
from time import time

mat = np.random.random_integers(0, 1, (40, 40))

N = 5
M = 15

window = np.ones((N, ), dtype=np.int)

start = time()

interm = convolve1d(mat, window, axis=0)
result = convolve1d(interm, window, axis=1)

[rows, cols] = np.where(result >= M)
result[:, :] = 0
result[(rows, cols)] = 1

end = time()
print "{} seconds".format(end - start)
0.00155591964722 seconds

Not sure how the complexity compares, but convolution is pretty well optimised in a variety of python deep learning libraries.

edited Jul 18 '18 at 03:34

answered Jul 18 '18 at 02:52

William Yolland

141
2

Since your 1d convolution kernels are both of the form `[1, 1, 1, 1, ..., 1]`, you can compute them incrementally in O(N) time (when applied to a row or column of length N). I guess numpy can't perform this optimization automatically. – Paul Hankin Jul 18 '18 at 12:54
I absolutely haven't thought about using convolution, that's a terrific idea ! Let's see how it compares ! – Jonathan DEKHTIAR Jul 18 '18 at 12:54
1

Damned, for information @william-yolland, your proposal is 10 TIMES faster :D – Jonathan DEKHTIAR Jul 18 '18 at 13:29

Any idea to optimise this algorithm?

Edit: Some code to make this more easy to understand.

1 Answers1