快速实现1参数分割失真模型

Question

我已经实现了一些代码，允许我估计图像的畸变中心以及试图消除图像中存在的任何径向畸变的畸变参数。它是根据 1 参数除法模型建模的：

xu = xd / (1 + K1 * r^2)

yu = yd / (1 + K1 * r^2)

我想知道尽快处理此转换的最佳方法是什么。我在想应该有一种方法可以使用 numpy 矢量化来为给定图像创建映射，我可以用它来执行不失真过程。我有 OpenCV 的经验（un Distor()、initUn DistorifyMap() 等），但是，这些方法需要估计相机焦点属性 (

fx, fy

)，而我没有。

我可以实现一个嵌套的 for 循环结构来计算给定输入图像大小的所有未失真目标，但是，这效率非常低。有没有一种方法可以用 numpy 来构建这个映射？

这是我如何将单个点从扭曲状态转换为未扭曲状态。

def get_undistorted(pd, dc, k):
    xd, yd = pd
    dcx, dcy = dc
    r2 = (dcx - xd)**2 + (dcy - yd)**2
    xu = dcx + (xd - dcx)/(1 + k*r2)
    yu = dcy + (yd - dcy)/(1 + k*r2)
    return xu, yu

Answer 1

我不能保证这是应用变换的最有效方法，但它是完全矢量化的，并且插值是可配置的。坐标空间是预先计算的，因此每帧只需要执行一次插值操作。

from typing import NamedTuple

import numpy as np
from PIL import Image
from matplotlib import pyplot as plt
from scipy.interpolate import RegularGridInterpolator


class Transform(NamedTuple):
    mn_high: tuple[int, int]
    interp_grid: tuple[np.ndarray, np.ndarray]
    yxdistort: np.ndarray
    yxfixed: np.ndarray
    rect_axes: tuple[np.ndarray, np.ndarray]

    @classmethod
    def setup(
        cls,
        dc: np.ndarray,
        k: float,
        mn: tuple[int, int],
        scale: float,
    ) -> 'Transform':
        ylong = np.arange(0, mn[0], scale)
        xlong = np.arange(0, mn[1], scale)
        yxrect_d3 = np.stack(
            np.meshgrid(xlong, ylong), axis=-1,
        )
        yxrect = yxrect_d3.reshape((-1, 2))[:, ::-1]

        interp_grid = tuple(
            dim.ravel()
            for dim in np.indices(dimensions=mn, sparse=True)
        )

        return cls(
            mn_high=yxrect_d3.shape[:2],
            interp_grid=interp_grid,
            yxdistort=redistort(dc=dc, k=k, vu=yxrect),
            yxfixed=undistort(dc=dc, k=k, yx=yxrect),
            rect_axes=(ylong, xlong),
        )

    def distort_image(self, image: np.ndarray) -> np.ndarray:
        distort_interp = RegularGridInterpolator(
            points=self.interp_grid, values=image, bounds_error=False,
        )
        return distort_interp(self.yxdistort).reshape(
            (*self.mn_high, -1)
        )

    def undistort_image(self, image: np.ndarray) -> np.ndarray:
        undistort_interp = RegularGridInterpolator(
            points=self.rect_axes, values=image, bounds_error=False,
        )
        return undistort_interp(self.yxfixed).reshape(
            (*self.mn_high, -1)
        )


def undistort(
    dc: np.ndarray,  # centre point, *2
    k: float,        # distort parameter aka. lambda
    yx: np.ndarray,  # distorted, n * 2
) -> np.ndarray:     # n * 2
    r2 = ((yx - dc)**2).sum(axis=1, keepdims=True)
    vu = dc + (yx - dc)/(1 + k*r2)
    return vu


def redistort(
    dc: np.ndarray,  # centre point, *2
    k: float,        # distort parameter aka. lambda
    vu: np.ndarray,  # non-distorted, n * 2
) -> np.ndarray:     # n * 2
    inner = k*((vu - dc)**2).sum(axis=1, keepdims=True)

    return (
        k*(
            dc*(
                + (vu**2).sum(axis=1, keepdims=True)
                - 2*dc*vu
                + dc**2
                + dc[::-1]**2
            )
            - 2*dc.prod()*vu[:, ::-1]
        )
        + 0.5*(vu - dc)*(
            1 - np.sqrt(1 - 4*inner)
        )
    )/inner


def regression_test() -> None:
    dcx = 320/2
    dcy = 240/2
    xd = 5
    yd = 17
    lambd = 1e-6

    r2 = (dcx - xd)**2 + (dcy - yd)**2
    xu = dcx + (xd - dcx)/(1 + lambd*r2)
    yu = dcy + (yd - dcy)/(1 + lambd*r2)

    dc = np.array((dcy, dcx))
    actual = undistort(
        dc=dc, k=lambd, yx=np.array((
            (yd, xd),
        )),
    )
    assert np.allclose(actual, (yu, xu))

    vu = np.array((
        (yu, xu),
        (40, 30),
        (25, 20),
    ))
    expected = np.array((
        (         yd,          xd),
        (38.04372287, 26.82104966),
        (22.11282237, 15.74521192),
    ))
    actual = redistort(dc, lambd, vu)
    assert np.allclose(expected, actual, rtol=0, atol=1e-8)


def estimate_resolution_loss(
    mn: tuple[int, int],
    dc: np.ndarray,
    k: float,
) -> float:
    # Pretty bad! This is not efficient.
    n = 200
    middle_strip = np.empty((n, 2))
    middle_strip[:, 0] = np.linspace(0, mn[0]-1, n)
    middle_strip[:, 1] = mn[1]/2

    distorted = redistort(dc=dc, k=k, vu=middle_strip)
    y, x = distorted.T
    y = y[np.isfinite(y)]

    res_reduce = mn[0]/(y.max() - y.min())
    return res_reduce


def roundtrip_demo(filename: str) -> None:
    # Once per parameter+resolution set
    imorig = np.array(Image.open(filename))
    mn = imorig.shape[:2]
    dc = 0.5*np.array(mn) + 0.1  # avoid a divide-by-zero
    k = 1e-6
    res_reduce = estimate_resolution_loss(mn=mn, dc=dc, k=k)
    transform = Transform.setup(mn=mn, dc=dc, k=k, scale=res_reduce)

    # Once every frame
    imdistort = transform.distort_image(imorig/255)
    imfixed = transform.undistort_image(imdistort)

    fig, (ax_orig, ax_distort, ax_fixed) = plt.subplots(ncols=3)
    ax_orig.imshow(imorig)
    ax_distort.imshow(imdistort)
    ax_fixed.imshow(imfixed)
    ax_orig.set_title('Original')
    ax_distort.set_title('Distorted')
    ax_fixed.set_title('Corrected')
    plt.show()


if __name__ == '__main__':
    regression_test()
    roundtrip_demo('baby goat.jpg')

可以采取更多措施更好地处理边缘等