Create a structuringselement in the form of a line with a certain degree and length

Question

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?

Answer 1

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.