3
\$\begingroup\$

Within Matlab there exist a function to create a structuringselement in the form of a line.
SE = strel('line',len,deg) creates a linear structuring element that is symmetric with respect to the neighborhood center, with approximate length len and angle deg. For a couple of morphological operations I need to filter out horizontal-, vertical and other angles of lines. Within Python's CV2 a similiar function does not exist. Therefore, I wanted to recreate this function in Python

My code:

import math
import numpy as np

# bresenham function is the accepted answer of SO's post https://stackoverflow.com/questions/23930274/list-of-coordinates-between-irregular-points-in-python
def bresenham(x0, y0, x1, y1):
   points = []
   dx = abs(x1 - x0)
   dy = abs(y1 - y0)
   x, y = x0, y0
   sx = -1 if x0 > x1 else 1
   sy = -1 if y0 > y1 else 1
   if dx > dy:
      err = dx / 2.0
      while x != x1:
         points.append((x, y))
         err -= dy
         if err < 0:
            y += sy
            err += dx
         x += sx
   else:
      err = dy / 2.0
      while y != y1:
         points.append((x, y))
         err -= dx
         if err < 0:
            x += sx
            err += dy
         y += sy
   points.append((x, y))

   return points


def strel_line(length, degrees):
   if length >= 1:
      theta = degrees * np.pi / 180
      x = round((length - 1) / 2 * np.cos(theta))
      y = -round((length - 1) / 2 * np.sin(theta))
      points = bresenham(-x, -y, x, y)
      points_x = [point[0] for point in points]
      points_y = [point[1] for point in points]
      n_rows = int(2 * max([abs(point_y) for point_y in points_y]) + 1)
      n_columns = int(2 * max([abs(point_x) for point_x in points_x]) + 1)
      strel = np.zeros((n_rows, n_columns))
      rows = ([point_y + max([abs(point_y) for point_y in points_y]) for point_y in points_y])
      columns = ([point_x + max([abs(point_x) for point_x in points_x]) for point_x in points_x])
      idx = []
      for x in zip(rows, columns):
         idx.append(np.ravel_multi_index((int(x[0]), int(x[1])), (n_rows, n_columns)))
      strel.reshape(-1)[idx] = 1

   return strel


if __name__=='__main__':
   strel = strel_line(15, 105)
   print(strel)
   print(strel.shape)

Output:

[[1. 0. 0. 0. 0.]
 [1. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0.]
 [0. 1. 0. 0. 0.]
 [0. 1. 0. 0. 0.]
 [0. 1. 0. 0. 0.]
 [0. 0. 1. 0. 0.]
 [0. 0. 1. 0. 0.]
 [0. 0. 1. 0. 0.]
 [0. 0. 0. 1. 0.]
 [0. 0. 0. 1. 0.]
 [0. 0. 0. 1. 0.]
 [0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 1.]
 [0. 0. 0. 0. 1.]]

Do you see any ways to write my function more efficiently?

\$\endgroup\$
2
\$\begingroup\$

Type hints

def bresenham(x0, y0, x1, y1):

should be, at a guess,

def bresenham(x0: float, y0: float, x1: float, y1: float) -> List[float]:

or maybe int, but you get the idea.

Type promotion

You don't need decimals here:

err = dx / 2.0

/2 will accomplish the same thing.

Termination

What happens if you call bresenham(0, 0, 9.5, 9.5)? I suspect that it will iterate infinitely and hang, due to

  while x != x1:

One way to fix this, if your algorithm permits, is - rather than making an sx dependent on the comparison of x0 and x1, simply swap x0 and x1 if they're in the wrong order. Then, rather than !=, use <, and you don't have to calculate sx.

Use a generator

For some use cases this will dramatically decrease memory consumption. Instead of building up points in memory, just yield whenever you have a point. Allow the caller to decide whether they want to cast to a list, a tuple, or just iterate over the points.

Int casting

  n_rows = int(2 * max([abs(point_y) for point_y in points_y]) + 1)

is overambitious in application of the int cast. The multiplication and addition don't need to be done in floating-point math, so you can instead do

  n_rows = 2 * int(max(abs(point_y) for point_y in points_y)) + 1

Also note that you shouldn't pass a list to max - just pass the generator directly.

\$\endgroup\$

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.