Initial commit

.gitignore

+poetry.lock
\ No newline at end of file

pyproject.toml

+[tool.poetry]
+name = "imagelab"
+version = "0.1.0"
+description = "This package is DTU Compute imagelab SDK for research in image analysis, computer vision, and computational imaging."
+authors = ["Søren K. S. Gregersen <sorgre@dtu.dk>"]
+license = "MIT"
+
+[tool.poetry.dependencies]
+python = "^3.7"
+opencv-contrib-python = "^4.4.0"
+numpy = "^1.17.2"
+scipy = "^1.5.2"
+
+[tool.poetry.dev-dependencies]
+
+[build-system]
+requires = ["poetry>=0.12"]
+build-backend = "poetry.masonry.api"

src/imagelab/__init__.py

+from . import io
+from . import phaseshift
+from . import construct
+from .camera import camera
+from . import utilities

src/imagelab/camera.py

+import numpy as np
+import cv2
+
+
+class camera:
+    def __init__(self, K, R=None, t=None, distortion=None, imshape=None):
+        self.K = K
+        self.R = R
+        self.t = t
+        self.distortion = distortion
+        self.imshape = imshape
+        if R is None:
+            self.R = np.eye(3, dtype=K.dtype)
+        if t is None:
+            self.t = np.zeros((3,), K.dtype)
+        if distortion is None:
+            self.distortion = np.array(tuple(), K.dtype)
+        if imshape is None:
+            self.imshape = np.array((0, 0), int)
+
+    @property
+    def position(self):
+        """position of the camera in *world space*."""
+        return - self.R.T.dot(self.t)
+
+    @property
+    def P(self):
+        """Projection matrix, (3, 4) array."""
+        return self.K.dot(np.c_[self.R, self.t])
+
+    @property
+    def focal_vector(self):
+        """Focal vector, (2,) view into array self.K."""
+        return np.diag(self.K)[:2]
+
+    @property
+    def principal_point(self):
+        """Principal point, (2,) view into array self.K."""
+        return self.K[:2, -1]
+
+    def relative_to(self, other):
+        R = self.R.dot(other.R.T)
+        t = R.dot(other.position - self.position)
+        return R, t
+
+    def __repr__(self):
+        def arr2str(s, A):
+            return s + np.array2string(A, precision=2, separator=',',
+                                       suppress_small=True,
+                                       prefix=s.strip('\n'))
+        return (arr2str("camera{   K: ", self.K) + ",\n" +
+                arr2str("          R: ", self.R) + ",\n" +
+                arr2str("          t: ", self.t) + ",\n" +
+                arr2str("    imshape: ", self.imshape) + ",\n" +
+                arr2str(" distortion: ", self.distortion) + "}")

src/imagelab/construct.py

+import numpy as np
+import plyfile
+import warnings
+from imagelab import utilities
+
+
+def fit_planes(points, mask=None):
+    barycenters = points.mean(axis=-2, keepdims=True)
+    baryvectors = (points - barycenters)
+    if mask is not None:
+        baryvectors[np.logical_not(mask)] *= 0
+    M = (baryvectors[..., None, :] * baryvectors[..., None]).sum(axis=-3)
+    eig_values = np.full((len(baryvectors), 3), np.nan)
+    eig_vectors = np.full((len(baryvectors), 3, 3), np.nan)
+    v = ~np.isnan(M).any(axis=-1).any(axis=-1)
+    eig_values[v], eig_vectors[v] = np.linalg.eigh(M[v])
+    i = tuple(np.arange(0, eig_values.shape[i], dtype=int)
+              for i in range(0, len(eig_values.shape) - 1))
+    _indices = np.zeros(len(eig_vectors), int)
+    _indices[v] = np.nanargmin(np.abs(eig_values[v]), axis=-1)
+    indices = (*i, slice(None), _indices)
+    return eig_vectors[indices]
+
+
+def normalize(v):
+    return v / np.linalg.norm(v, axis=-1, keepdims=True)
+
+
+def planefit_normals(camera, points3D, pixels):
+    imshape = tuple(camera.imshape[:2])
+    xyz = np.full(imshape + (3,), np.nan, points3D.dtype)
+    xyz_pixels = tuple(utilities.rint(pixels).T)
+    xyz[xyz_pixels] = points3D
+    mask = np.zeros(imshape + (3, 3), bool)
+    for i in [-1, 0, 1]:
+        for j in [-1, 0, 1]:
+            mask_pixels = tuple(utilities.rint(pixels + np.array([[i, j]])).T)
+            mask[(*mask_pixels, i, j)] = 1
+    mask = mask.reshape(imshape + (9,))
+
+    p3D = np.stack([xyz[mask[:, :, i]] for i in range(9)], axis=-2)
+    nmask = np.stack([mask[mask[:, :, i]][:, 4] for i in range(9)], axis=-1)
+
+    n = normalize(fit_planes(p3D, nmask))
+
+    nan_count = np.isnan(n).sum()
+    if nan_count > 0:
+        print('Warning: from_SVD() produced {} nans'.format(nan_count))
+
+    c = normalize(camera.position - xyz[mask[:, :, 4]])
+    deviation = (n * c).sum(axis=1)
+    n *= np.sign(deviation)[:, None]
+    return n, np.abs(deviation)
+
+
+class pointcloud(object):
+    """pointcloud encapsulates positions, normals, and colors.
+
+    The class can read and write Standford .ply files"""
+
+    def __init__(self, positions=None, colors=None, normals=None):
+        super(pointcloud, self).__init__()
+        self.positions = positions
+        self.colors = colors
+        self.normals = normals
+
+    def writePLY(self, filename, ascii=False):
+        dtype = []
+        N = -1
+        if self.positions is not None:
+            N = len(self.positions)
+            dtype += [('x', 'f4'), ('y', 'f4'), ('z', 'f4')]
+        if self.colors is not None:
+            N = len(self.colors) if N == -1 else N
+            dtype += [('red', 'u1'), ('green', 'u1'), ('blue', 'u1')]
+        if self.normals is not None:
+            N = len(self.normals) if N == -1 else N
+            dtype += [('nx', 'f4'), ('ny', 'f4'), ('nz', 'f4')]
+
+        error_msg = "Lengths of positions, colors, and normals must match."
+        if self.positions is not None and N != len(self.positions):
+            raise RuntimeError(error_msg)
+        if self.colors is not None and N != len(self.colors):
+            raise RuntimeError(error_msg)
+        if self.normals is not None and N != len(self.normals):
+            raise RuntimeError(error_msg)
+
+        vertex = np.zeros((N,), dtype=dtype)
+
+        if self.positions is not None:
+            vertex['x'] = self.positions[:, 0].astype('f4')
+            vertex['y'] = self.positions[:, 1].astype('f4')
+            vertex['z'] = self.positions[:, 2].astype('f4')
+
+        if self.colors is not None:
+            self.colors = np.array(self.colors)
+            if self.colors.shape == (N,) or self.colors.shape == (N, 1):
+                vertex['red'] = self.colors.squeeze().astype('u1')
+                vertex['green'] = self.colors.squeeze().astype('u1')
+                vertex['blue'] = self.colors.squeeze().astype('u1')
+            elif self.colors.shape == (N, 3):
+                vertex['red'] = self.colors[:, 0].astype('u1')
+                vertex['green'] = self.colors[:, 1].astype('u1')
+                vertex['blue'] = self.colors[:, 2].astype('u1')
+
+        if self.normals is not None:
+            vertex['nx'] = self.normals[:, 0].astype('f4')
+            vertex['ny'] = self.normals[:, 1].astype('f4')
+            vertex['nz'] = self.normals[:, 2].astype('f4')
+
+        vertex = plyfile.PlyElement.describe(vertex, 'vertex')
+
+        ext = filename.split('.')[-1]
+        if ext != "ply" and ext != "PLY":
+            filename = filename + '.ply'
+        plyfile.PlyData([vertex], text=ascii).write(filename)
+        return self
+
+    def readPLY(self, filename):
+        self.__init__()
+
+        vertex = plyfile.PlyData.read(filename)['vertex']
+
+        with warnings.catch_warnings():
+            # numpy does not like to .view() into structured array
+            warnings.simplefilter("ignore")
+
+            if all([p in vertex.data.dtype.names for p in ('x', 'y', 'z')]):
+                position_data = vertex.data[['x', 'y', 'z']]
+                N = len(position_data.dtype.names)
+                self.positions = position_data.view((position_data.dtype[0],
+                                                     N))
+
+            colored = all([p in vertex.data.dtype.names
+                           for p in ('red', 'green', 'blue')])
+            if colored:
+                color_data = vertex.data[['red', 'green', 'blue']]
+                N = len(color_data.dtype.names)
+                self.colors = color_data.view((color_data.dtype[0], N))
+
+            if all([p in vertex.data.dtype.names for p in ('nx', 'ny', 'nz')]):
+                normal_data = vertex.data[['nx', 'ny', 'nz']]
+                N = len(normal_data.dtype.names)
+                self.normals = normal_data.view((normal_data.dtype[0], N))
+        return self

src/imagelab/io.py

+import numpy as np
+import cv2
+import os
+import re
+
+
+def imread(path, color=True):
+    if not os.path.exists(path):
+        raise FileNotFoundError("No such file: '" + path + "'.")
+    if color:
+        opts = cv2.IMREAD_COLOR | cv2.IMREAD_ANYDEPTH
+        return cv2.imread(path, opts)[:, :, ::-1]  # BGR -> RGB
+    else:
+        opts = cv2.IMREAD_GRAYSCALE | cv2.IMREAD_ANYDEPTH
+        return cv2.imread(path, opts)[:, :, None]
+
+
+def imwrite(path, image):
+    image = image.reshape((*image.shape[:2], -1))  # Make at least 3D.
+    return cv2.imwrite(path, image[:, :, ::-1])  # RGB -> BGR
+
+
+def read_images(path, rexp=r'.*\.png', sort=None, filter=None, color=True):
+    files = [os.path.join(path, f)
+             for f in os.listdir(path)
+             if re.match(rexp, f)]
+    if filter is not None:
+        files = filter(files)
+    if sort is not None:
+        files = sort(files)
+    if files == []:
+        return []
+    image0 = imread(files[0], color)
+    images = np.empty((len(files), *image0.shape), dtype=image0.dtype)
+    images[0] = image0
+    for i, path in enumerate(files[1:], 1):
+        images[i] = imread(path, color)
+    return images
+
+
+def write_images(directory, images, format=None):
+    images = np.asarray(images)
+    shape = images.shape[:-3]
+    if format is None:
+        format = "image" + ("_{}" * len(shape)) + ".png"
+    path = os.path.join(directory, format)
+    for idx in np.ndindex(shape):
+        filepath = path.format(*idx)
+        directory, _ = os.path.split(filepath)
+        os.makedirs(directory, exist_ok=True)
+        imwrite(filepath, images[idx])

src/imagelab/phaseshift.py

+import numpy as np
+from scipy import ndimage
+from scipy.fftpack import fft
+from imagelab.utilities import ndtake, fint
+from imagelab import io  # for debugging purposes.
+
+
+def decode(data, DFT_channels=1):
+    data = np.squeeze(data, axis=-1)  # assume grayscale and squeese axis
+    shifting_axis = -3  # always assume lists of images
+    dft = fft(data, axis=shifting_axis)
+
+    N = data.shape[shifting_axis]  # N data points
+    sdft = 4 * dft / N
+    M = fint(N / 2)  # Nyquist frequency
+
+    sdft_ch = sdft[ndtake(DFT_channels, shifting_axis)]
+    phase = np.mod(-np.angle(sdft_ch), 2 * np.pi)
+    signal = np.absolute(sdft_ch)
+    sdft = sdft[ndtake(slice(1, M), shifting_axis)]  # keep until Nyquist
+    noise = np.absolute(sdft).sum(axis=shifting_axis)
+    noise = noise - signal  # noise is all minus the *individual* signals
+    return phase, signal, noise
+
+
+def unwrap_phase_with_cue(phase, cue, wave_count):
+    phase_cue = np.mod(cue - phase, 2 * np.pi)
+    phase_cue = np.round(((phase_cue * wave_count) - phase) / (2 * np.pi))
+    return (phase + (2 * np.pi * phase_cue)) / wave_count
+
+
+def stdmask(phase, dphase, signal_noise_list):
+    """ Standard masking function to use on the phase shift.
+
+        If you wish to change how it works or some parameters, then copy this
+        funtion and replace the relevant sections."""
+    mask = np.ones_like(phase, dtype=bool)
+    # Threshold on phase
+    mask &= np.logical_and(phase >= 0.0, phase <= 2 * np.pi)
+    for signal, noise in signal_noise_list:
+        # Threshold on signal trust region
+        mask &= np.logical_and(0.1 < signal, signal < 0.9)
+        # Threshold on signal to noise ratio
+        mask &= signal > 4 * noise
+    # Threshold on gradient and curvature of phase.
+    # TODO: Tune the tolarances for gradients and curvature.
+    dph = np.linalg.norm(dphase, axis=-1)
+    ddph = np.stack(np.gradient(dph, axis=(-2, -1)), axis=-1)
+    ddph = np.linalg.norm(ddph, axis=-1)
+    mask &= np.logical_and(1e-8 < dph, dph < 1e-2)
+    mask &= ddph < 1e-2
+    # Remove borders
+    mask[..., [0, -1]] = 0
+    mask[..., [0, -1], :] = 0
+    # Remove mask-pixels with masked neighbors
+    weights = np.array([[1, 1, 1], [1, 0, 1], [1, 1, 1]], dtype=int)
+    weights = np.broadcast_to(weights, mask.shape[:-2] + weights.shape)
+    neighbors = ndimage.convolve(mask.astype(int), weights, mode='constant')
+    mask = np.logical_and(mask, neighbors > 1)
+    return mask
+
+
+def decode_with_cue(primary, cue, wave_count, DFT_channels=1, maskfn=stdmask):
+    primary_phase, primary_signal, primary_noise = decode(primary)
+    cue_phase, cue_signal, cue_noise = decode(cue)
+    phase = unwrap_phase_with_cue(primary_phase, cue_phase, wave_count)
+    dphase = np.stack(np.gradient(phase, axis=(-2, -1)), axis=-1)
+    signal_noise_list = [[primary_signal, primary_noise],
+                         [cue_signal, cue_noise]]
+    mask = maskfn(phase, dphase, signal_noise_list)
+    return phase, dphase, mask
+
+
+def stereo_correspondances(L_enc, R_enc, L_mask, R_mask, keep_unsure=False):
+    if L_mask.sum() == 0 or R_mask.sum() == 0:
+        return np.empty((2, 0))
+    # To conserve computation size, we only consider a nonzero encoding patch.
+    L_nz_rows, L_nz_cols = L_mask.nonzero()
+    R_nz_rows, R_nz_cols = R_mask.nonzero()
+
+    # Notice that rows are shared between L and R.
+    rows = slice(min(L_nz_rows.min(), R_nz_rows.min()),
+                 max(L_nz_rows.max(), R_nz_rows.max()) + 1)
+    L_cols = slice(L_nz_cols.min(), L_nz_cols.max() + 1)
+    R_cols = slice(R_nz_cols.min(), R_nz_cols.max() + 1)
+
+    L_mask, R_mask = L_mask[rows, L_cols], R_mask[rows, R_cols]
+    L_enc, R_enc = L_enc[rows, L_cols], R_enc[rows, R_cols]
+
+    # both the left and right side of a match must exist
+    match = np.logical_and(L_mask[..., None], R_mask[..., None, :-1])
+    match &= R_mask[..., None, 1:]
+
+    # Value in left should lie between to neighbouring values in right
+    match[L_enc[..., :, None] < R_enc[..., None, :-1]] = 0
+    match[L_enc[..., :, None] >= R_enc[..., None, 1:]] = 0
+
+    if not keep_unsure:
+        match[match.sum(axis=-1) > 1] = 0
+
+    if match.sum() == 0:
+        return np.empty((2, 0))
+
+    if np.any(match.sum(axis=-1) > 1):
+        print('wrong match', (match.sum(axis=-1) > 1).sum())
+        errors = (match.sum(axis=-1) > 1)
+        for e in errors.nonzero()[0]:
+            print('i', e, 'phase', L_enc[..., e])
+            _e = match[e].nonzero()[0]
+            _a, _b = np.min(_e), np.max(_e) + 1
+            print('a,b', _a, _b)
+            print('left', R_enc[..., :-1][..., _a:_b])
+            print('right', R_enc[..., 1:][..., _a:_b])
+
+    r, cL, cR = tuple(match.nonzero())
+    step = R_enc[(r, cR + 1)] - R_enc[(r, cR)]
+    cRfrac = (L_enc[(r, cL)] - R_enc[(r, cR)]) / step
+
+    r += rows.start
+    cL += L_cols.start
+    cR = cR + R_cols.start + cRfrac
+    return np.array([[r, cL], [r, cR]]).swapaxes(1, 2)
+
+
+def stereo_triangulate_linearization(left_cam, right_cam):
+    dtype = left_cam.P.dtype
+    P0 = left_cam.P
+    P1 = right_cam.P
+    e = np.eye(4, dtype=dtype)
+    C = np.empty((4, 3, 3, 3), dtype)
+    for i in np.ndindex((4, 3, 3, 3)):
+        tmp = np.stack((P0[i[1]], P0[i[2]], P1[i[3]], e[i[0]]), axis=0)
+        C[i] = np.linalg.det(tmp.T)
+    C = C[..., None, None]
+    yx = np.mgrid[0:left_cam.imshape[0], 0:left_cam.imshape[1]].astype(dtype)
+    y, x = yx[None, 0, :, :], yx[None, 1, :, :]
+    offset = C[:, 0, 1, 0] - C[:, 2, 1, 0] * x - C[:, 0, 2, 0] * y
+    factor = -C[:, 0, 1, 2] + C[:, 2, 1, 2] * x + C[:, 0, 2, 2] * y
+    return offset, factor
+
+
+def stereo_triangulate(left, right, left_cam, right_cam):
+    offset, factor = stereo_triangulate_linearization(left_cam, right_cam)
+    idx = (slice(None), *(left + 0.5).astype(int).T)
+    xyzw = offset[idx] + factor[idx] * right[None, ..., 1]
+    return xyzw.T[..., :3] / xyzw.T[..., 3, None]

src/imagelab/pointcloud.py

+import numpy as np
+import plyfile
+import warnings
+
+
+class pointcloud(object):
+    """pointcloud encapsulates positions, normals, and colors.
+
+    The class can read and write Standford .ply files"""
+
+    def __init__(self, positions=None, colors=None, normals=None):
+        super(pointcloud, self).__init__()
+        self.positions = positions
+        self.colors = colors
+        self.normals = normals
+
+    def writePLY(self, filename, ascii=False):
+        dtype = []
+        N = -1
+        if self.positions is not None:
+            N = len(self.positions)
+            dtype += [('x', 'f4'), ('y', 'f4'), ('z', 'f4')]
+        if self.colors is not None:
+            N = len(self.colors) if N == -1 else N
+            dtype += [('red', 'u1'), ('green', 'u1'), ('blue', 'u1')]
+        if self.normals is not None:
+            N = len(self.normals) if N == -1 else N
+            dtype += [('nx', 'f4'), ('ny', 'f4'), ('nz', 'f4')]
+
+        error_msg = "Lengths of positions, colors, and normals must match."
+        if self.positions is not None and N != len(self.positions):
+            raise RuntimeError(error_msg)
+        if self.colors is not None and N != len(self.colors):
+            raise RuntimeError(error_msg)
+        if self.normals is not None and N != len(self.normals):
+            raise RuntimeError(error_msg)
+
+        vertex = np.zeros((N,), dtype=dtype)
+
+        if self.positions is not None:
+            vertex['x'] = self.positions[:, 0].astype('f4')
+            vertex['y'] = self.positions[:, 1].astype('f4')
+            vertex['z'] = self.positions[:, 2].astype('f4')
+
+        if self.colors is not None:
+            # assuming RGB format
+            vertex['red'] = self.colors[:, 0].astype('u1')
+            vertex['green'] = self.colors[:, 1].astype('u1')
+            vertex['blue'] = self.colors[:, 2].astype('u1')
+
+        if self.normals is not None:
+            vertex['nx'] = self.normals[:, 0].astype('f4')
+            vertex['ny'] = self.normals[:, 1].astype('f4')
+            vertex['nz'] = self.normals[:, 2].astype('f4')
+
+        vertex = plyfile.PlyElement.describe(vertex, 'vertex')
+
+        ext = filename.split('.')[-1]
+        if ext != "ply" and ext != "PLY":
+            filename = filename + '.ply'
+        plyfile.PlyData([vertex], text=ascii).write(filename)
+        return self
+
+    def readPLY(self, filename):
+        self.__init__()
+
+        vertex = plyfile.PlyData.read(filename)['vertex']
+
+        with warnings.catch_warnings():
+            # numpy does not like to .view() into structured array
+            warnings.simplefilter("ignore")
+
+            if all([p in vertex.data.dtype.names for p in ('x', 'y', 'z')]):
+                position_data = vertex.data[['x', 'y', 'z']]
+                N = len(position_data.dtype.names)
+                self.positions = position_data.view((position_data.dtype[0],
+                                                     N))
+
+            colored = all([p in vertex.data.dtype.names
+                           for p in ('red', 'green', 'blue')])
+            if colored:
+                color_data = vertex.data[['red', 'green', 'blue']]
+                N = len(color_data.dtype.names)
+                self.colors = color_data.view((color_data.dtype[0], N))
+
+            if all([p in vertex.data.dtype.names for p in ('nx', 'ny', 'nz')]):
+                normal_data = vertex.data[['nx', 'ny', 'nz']]
+                N = len(normal_data.dtype.names)
+                self.normals = normal_data.view((normal_data.dtype[0], N))
+        return self

src/imagelab/utilities.py

+import numpy as np
+import os
+import cv2
+
+
+def rint(array, dtype=int):
+    return np.rint(array).astype(dtype)
+
+
+def fint(array, dtype=int):
+    return np.floor(array).astype(dtype)
+
+
+def cint(array, dtype=int):
+    return np.ceil(array).astype(dtype)
+
+
+def ndtake(slicing, axis):  # simple slicing, much faster than np.take
+    if axis >= 0:
+        return (slice(None),) * axis + (slicing,)
+    if axis < 0:
+        return (Ellipsis, slicing,) + (slice(None),) * (-1 - axis)
+
+
+def ndsplit(array, widths, axis=0):  # split axis into len(widths) parts
+    array = np.asanyarray(array)
+    splits = (0,) + tuple(np.cumsum(widths))
+    return [array[ndtake(slice(s0, s1), axis)]
+            if s1 - s0 > 1 else array[ndtake(s0, axis)]
+            for s0, s1 in zip(splits[:-1], splits[1:])]
+
+
+def vectordot(u, v, *args, **kwargs):
+    """Specilization of the dot-operator. u and v are ndarrays of vectors"""
+    u, v = np.broadcast_arrays(u, v)
+    return np.einsum('...i,...i ->...', u, v).reshape(*u.shape[:-1], 1)
+
+
+def plot_correspondances(im_L, im_R, pixels_L, pixels_R, colors=None):
+    im_joined = np.concatenate([im_L, im_R], axis=1)
+    if len(im_joined.shape) == 2:
+        im_joined = im_joined[:, :, None]
+
+    pixels_L = np.array(pixels_L)
+    pixels_R = np.array(pixels_R)
+    if len(pixels_L) != len(pixels_R):
+        raise RuntimeError("Pixel lists must be same length")
+    elif pixels_L.shape[1] != 2 or pixels_R.shape[1] != 2:
+        raise RuntimeError("Pixel must have 2 coordinates")
+
+    if colors is None:
+        if im_joined.dtype == np.uint8:
+            colors = 255
+        else:
+            colors = 1
+    colors = np.array(colors)
+    if colors.ndim == 0 or colors.shape == (1,):
+        colors = [colors.item()] * im_joined.shape[2]
+    if len(colors) == len(pixels_L):
+        if colors.ndim == 1:
+            colors = colors[:, None]
+        elif colors.shape[1] != im_joined.shape[2]:
+            raise RuntimeError("Color(s) do not match channel depth")
+    if colors.shape == (im_joined.shape[2],):
+        colors = [colors] * len(pixels_L)
+    else:
+        raise RuntimeError("Color(s) do not match channel depth")
+    colors = np.array(colors)
+
+    start = np.array(pixels_L, dtype=np.float64)
+    end = np.array(pixels_R, dtype=np.float64) + np.array([[0, im_L.shape[1]]])
+    v = end - start
+    dist = np.linalg.norm(v, axis=1)
+    v /= dist[:, None]
+    dist -= 0.25
+
+    nx = np.array(start)
+    crit = np.linalg.norm(nx - start, axis=1) < dist
+    while np.any(crit):
+        idx = tuple(rint(nx[crit]).T)
+        im_joined[idx] = colors[crit]
+        nx += v
+        crit = np.linalg.norm(nx - start, axis=1) < dist
+
+    idx = tuple(rint(pixels_L).T)
+    im_joined[idx] = colors
+
+    return im_joined
+
+
+# Precomputed LUTs for combine_HDR
+triangle_weights_half = np.arange(1, 129)
+triangle_weights = np.empty((256,))
+triangle_weights[:128] = triangle_weights_half
+triangle_weights[128:] = triangle_weights_half[::-1]
+del triangle_weights_half
+linear_response = np.arange(0, 256, dtype=np.float32)
+linear_response[0] = 1
+log_linear_response = np.log(linear_response)
+del linear_response
+
+
+def combine_HDR_vectorized(frames, exposure_times=None, out=None, weight=None):
+    # debevec method (very memory intensive and slow. Dont know why.)
+    if exposure_times is None:
+        exposures = np.arange(0, frames.shape[0], dtype=np.float32) / 9
+    else:
+        exposures = -np.log(exposure_times)
+    if out is None:
+        out = np.zeros((frames.shape[1:]), dtype=np.float32)
+    if weight is None:
+        weight = np.zeros((*frames.shape[1:-1], 1), dtype=np.float32)
+    for i, exposure in enumerate(exposures):
+        frame = frames[i]
+        response = log_linear_response[frame] + exposure
+        weights = np.sum(triangle_weights[frame], axis=-1, keepdims=True)
+        weight += weights
+        out += weights * response
+    out /= weight
+    return np.exp(out)
+
+
+def combine_HDR(frames, exposure_times=None, step=1):
+    list_shape = frames.shape[1:-3]
+    if len(list_shape) == 0:
+        return combine_HDR_vectorized(frames, exposure_times)
+    if exposure_times is None:
+        exposure_times = np.arange(0, frames.shape[0], dtype=np.float32)
+        exposure_times = np.power(2, exposure_times) / 2**9
+    result = np.zeros((frames.shape[1:]), dtype=np.float32)
+    tmp_weights = np.zeros((step, *frames.shape[-3:-1], 1), dtype=np.float32)
+    iter = np.ndindex(list_shape)
+    for idx in iter:
+        [iter.next() for i in range(step - 1)]
+        idcs = idx[:-1] + (slice(idx[-1], idx[-1] + step),)
+        ith_frames = frames[(slice(None),) + idcs]
+        ith_tmp_weights = tmp_weights[:min(frames.shape[-4] - idx[-1], step)]
+        result[idcs] = combine_HDR_vectorized(ith_frames, exposure_times,
+                                              out=result[idcs],
+                                              weight=ith_tmp_weights)
+    return result
+
+
+def combine_HDR_using_opencv(frames, exposure_times=None):
+    list_shape = frames.shape[1:-3]
+    if len(list_shape) == 0:
+        return combine_HDR_using_opencv(frames[:, None], exposure_times)[0]
+    if exposure_times is None:
+        exposure_times = np.arange(0, frames.shape[0], dtype=np.float32)
+        exposure_times = 255 * np.power(2, exposure_times)
+    result = np.empty((frames.shape[1:]), dtype=np.float32)
+    merge_debvec = cv2.createMergeDebevec()
+    for idx in np.ndindex(list_shape):
+        frame = frames[(slice(None),) + idx]
+        result[idx] = merge_debvec.process(frame, times=exposure_times.copy())
+    return result