如何拟合仅由缩放和平移组成的仿射变换？

Question

我正在尝试编写一个应用程序，允许用户通过单击地图上的航点来对地图图像进行地理配准，从而确定像素坐标x 和y，然后单击“真实的' 地图导致经度 lng 和纬度 lng。

我认为图像和“真实”地图都是南北方向的，因此它们之间没有旋转，而且南北和东方向只有一个比例-西。那就是：

lng(x) = scale * x + a
lat(y) = -scale * y + b

（减号的原因是y像素坐标从图像顶部到底部增加，而纬度lat从南到北增加。

我已将how to perform coordinates affine transformation using python? part 2 的答案修改如下：

import numpy as np

coords = [
    {"pixel": {"x": 610, "y": 1673}, "lnglat": {
        "lng": -119.66622566136141, "lat": 37.71690293889708}},
    {"pixel": {"x": 3616, "y": 948}, "lnglat": {
        "lng": -119.55987333997541, "lat": 37.739791632115}},
    {"pixel": {"x": 156, "y": 1582}, "lnglat": {
        "lng": -119.68242540789811, "lat": 37.719168689576634}},
    {"pixel": {"x": 1432, "y": 1079}, "lnglat": {
        "lng": -119.63773163590452, "lat": 37.733899511112554}},
    {"pixel": {"x": 1467, "y": 982}, "lnglat": {
        "lng": -119.6365899951677, "lat": 37.73740878429034}},
    {"pixel": {"x": 2045, "y": 464}, "lnglat": {
        "lng": -119.61643210247348, "lat": 37.75263501532096}},
    {"pixel": {"x": 2530, "y": 225}, "lnglat": {
        "lng": -119.59904847563081, "lat": 37.759640099263024}},
    {"pixel": {"x": 3611, "y": 217}, "lnglat": {
        "lng": -119.57440674003465, "lat": 37.769372182124215}},
    {"pixel": {"x": 4218, "y": 289}, "lnglat": {
        "lng": -119.53927620600871, "lat": 37.7590418448261}},
    {"pixel": {"x": 4972, "y": 819}, "lnglat": {
        "lng": -119.51283799895947, "lat": 37.7451015130886}},
    {"pixel": {"x": 4869, "y": 1178}, "lnglat": {
        "lng": -119.5150031101931, "lat": 37.73452849532761}},
    {"pixel": {"x": 4858, "y": 1268}, "lnglat": {
        "lng": -119.51537412573026, "lat": 37.731943969799104}},
    {"pixel": {"x": 4637, "y": 1307}, "lnglat": {
        "lng": -119.52293169964986, "lat": 37.730726899819345}},
    {"pixel": {"x": 4284, "y": 1599}, "lnglat": {
        "lng": -119.53520554208092, "lat": 37.72240153076238}},
    {"pixel": {"x": 4150, "y": 1676}, "lnglat": {
        "lng": -119.53996905111126, "lat": 37.71984653680312}},
    {"pixel": {"x": 3432, "y": 1989}, "lnglat": {
        "lng": -119.56520552108367, "lat": 37.70994983543632}},
    {"pixel": {"x": 2965, "y": 1408}, "lnglat": {
        "lng": -119.58234774459186, "lat": 37.72663636959598}},
    {"pixel": {"x": 2560, "y": 1921}, "lnglat": {
        "lng": -119.59584076119313, "lat": 37.712008849961066}},
    {"pixel": {"x": 1840, "y": 1593}, "lnglat": {
        "lng": -119.6231396666414, "lat": 37.72018991118786}},
    {"pixel": {"x": 1140, "y": 1590}, "lnglat": {
        "lng": -119.64782744839357, "lat": 37.71938854312988}},
]

pixel_coordinates = np.array([[coord['pixel']['x'], coord['pixel']['y']] for coord in coords])
lnglat_coordinates = np.array([[coord['lnglat']['lng'], coord['lnglat']['lat']] for coord in coords])


# Pad the data with ones, so that our transformation can do translations too
n = pixel_coordinates.shape[0]


def pad(x):
    return np.hstack([x, np.ones((x.shape[0], 1))])


def unpad(x):
    return x[:, :-1]


X = pad(pixel_coordinates)
Y = pad(lnglat_coordinates)

# Solve the least squares problem X * A = Y
# to find our transformation matrix A
A, res, rank, s = np.linalg.lstsq(X, Y, rcond=None)


def transform(x):
    return unpad(np.dot(pad(x), A))


print("Target:")
print(lnglat_coordinates)
print("Result:")
print(transform(pixel_coordinates))
print("Max error:", np.abs(lnglat_coordinates - transform(pixel_coordinates)).max())

print(A)

生成的矩阵A 如下所示：

[[ 3.55857577e-05  8.98377941e-07  2.81630741e-18]
 [ 3.43101520e-06 -2.97714115e-05 -8.56519716e-18]
 [-1.19693144e+02  3.77657997e+01  1.00000000e+00]]

我希望它有这种形式

[[scale 0      0]
 [0     -scale 0]
 [a     b      1]]

我注意到scale 的独立派生值实际上相差约 16%，而应该为零的非对角线列却没有。

（在实践中，我还注意到使用此算法计算的叠加层明显偏离，如下所示，部分不透明。注意道路似乎向西北移动）。

有没有办法将这些约束强加于仿射变换的最小二乘估计 - 即确保变换仅包含缩放和平移（没有倾斜或旋转）？

Answer 1

我按照UC Davis lecture notes 中描述的方法解决了这个问题；我还在related question on Math StackExchange 中描述了这种方法。

假设你有两个航路点，那么经纬度与像素坐标的关系可以表示为

（我还意识到除非您在赤道，否则经度和纬度的比例因子是不同的）。这是ax = b 在numpy.linalg.lstsq 中描述的形式，因此我们可以在最小二乘意义上直接估计x。

这是我用来实现它的代码：

from dataclasses import dataclass
from typing import List, Dict
import numpy as np
from cached_property import cached_property


@dataclass
class MapFitter:
    coords: List[Dict[str, object]]

    @cached_property
    def coeffs(self):
        lnglat = np.vstack(tuple(
            np.array([[coord['lnglat']['lng']], [coord['lnglat']['lat']]]) for coord in coords))

        xy = np.vstack(tuple(
            np.array([[coord['pixel']['x'], 0, 1, 0], [0, coord['pixel']['y'], 0, 1]]) for coord in coords))

        coefficients, residuals, rank, singular_values = np.linalg.lstsq(
            a=xy, b=lnglat, rcond=None)

        return coefficients

    def transform(self, xy_coord):
        x, y = xy_coord
        scale_lng, scale_lat, offset_lng, offset_lat = self.coeffs
        return [scale_lng * x + offset_lng, scale_lat * y + offset_lat]

例如，使用问题中给出的相同coords 数据结构，我可以这样做

map_fitter = MapFitter(coords=coords)

print(map_fitter.transform([0, 0]))
print(map_fitter.transform([5038, 2027]))

导致地图的角坐标：

[array([-119.68835968]), array([37.7690185])]
[array([-119.51026687]), array([37.70773367])]

我还发现这可以使叠加层与地图的对齐效果更好。