关于根据鱼眼相机的uv坐标反推相机在3维场景下的入射向量

之前我写了一篇关于针孔相机和鱼眼相机从3维坐标转换为2d相机的uv坐标的公式步骤.
同样的,当我们已知图片的uv坐标时,是可以反推出在3d场景下以相机位置为起点的一条射线的(也就是得到一个向量值)
原理就是根据正向的公式进行反推,反推过程中也遇到一些坑,这里记录一下步骤和关键点
注意这里只是鱼眼相机的反推,针孔的畸变处理大同小异,实际项目里偏差很小,并未实现,将来如果实现了再补充.

1. 根据(u,v,1)乘以相机投影矩阵的逆矩阵得到3d场景的相机坐标系下的归一化坐标(x,y,1)
2. 根据x,y得到半径r_d=sqrt(x^2+y^2)
3. 根据牛顿迭代法根据r_d得到实际r的值
4. 得到r_d和r的比值scale = r_d / r
5. 根据比值scale得到真实的(x',y')=(x/scale,y/scale)
6. 根据真实的归一化坐标(x',y',1,1)乘以相机的平移旋转矩阵的逆矩阵得到世界坐标

这里最复杂的处理要数第3步根据r_d得到实际的r值,本质就是因为鱼眼相机导致的真实入射角偏移.
https://blog.csdn.net/qq_16137569/article/details/112398976

这篇博客的最后有提到是求thera角的一元高次方程,不过具体的实现公式并没有详细给出,还是稍微有一点坑需要注意的:
因为我们前面把相机3维相机归一化,所以这边根据半径的反正切值就是thera角度,具体公式为thera=arctan(r)
之前讲过opencv里面鱼眼相机的畸变矫正
r_d = theta(1+k1*theta^2+k2*theta^4+k3*theta^6+k4*theta^8)
所以得出r_d和r的关系是
r_d = arctan(r)+k1*arctan(r)^3+k2*arctan(r)^5+k3*arctan(r)^7+k4*arctan(r)^9
接下来要根据牛顿迭代法来求r的近似解,
牛顿迭代法公式为f(xn+1) = xn - f(xn)/f'(xn)
根据上面的公式,可以得到
f(r) = arctan(r)+k1*arctan(r)^3+k2*arctan(r)^5+k3*arctan(r)^7+k4*arctan(r)^9-r_d = 0
arctan(r)的导数为1/(1+x^2);x^n的导数为n*x^(n-1)
推出f'(r) = 1/(1+r^2)+3*k1*arctan(r)^2/(1+r^2)+5*k2*arctan(r)^4/(1+r^2)+7*k3*arctan(r)^6/(1+r^2)+9*k4*arctan(r)^8/(1+r^2)
然后第一步令 r0=r_d
计算r1=r0-f(r0)/f'(r0)
然后计算r1-r0的绝对值,如果值大于阈值(如0.000001),则计算r2=r1-f(r1)/f'(r1)
以此类推一直迭代到r(n)-r(n-1)的绝对值在阈值以下(非常接近0)或者迭代了很多次(如15次)

实现的关键代码

  const r_d = Math.sqrt(Math.pow(sx, 2) + Math.pow(sy, 2))

  const e = 1-e6
  let c = 0
  let x = r_d

  let y = getNextX(r_d, x, k1, k2, k3, k4)

  while (Math.abs(x - y) > e && c < 10) {
    x = y
    y = getNextX(r_d, x, k1, k2, k3, k4)
    c++
  }

  const r = x
  const s = r_d / r
  // 修正值
  const rx = sx / s
  const ry = sy / s



function getNextX(r_d: number, x: number, k1: number, k2: number, k3: number, k4: number) {
  const atanX = Math.atan(x)
  const derAtanX = 1 / (1 + Math.pow(x, 2)) // 1/(1+x^2)
  const fx =
    atanX + k1 * Math.pow(atanX, 3) + k2 * Math.pow(atanX, 5) + k3 * Math.pow(atanX, 7) + k4 * Math.pow(atanX, 9) - r_d

  const der_fx =
    derAtanX +
    3 * k1 * Math.pow(atanX, 2) * derAtanX +
    5 * k2 * Math.pow(atanX, 4) * derAtanX +
    7 * k3 * Math.pow(atanX, 6) * derAtanX +
    9 * k4 * Math.pow(atanX, 8) * derAtanX

  const nextX = x - fx / der_fx

  return nextX
}

去年一年有点忙,很长一段时间没有写博客,新的一年里我要立个flag给自己找点事情干干:
去年在公司的项目里我实现了一套基于2d canvas的图形框架库.自我感觉做的还可以,但是还是有不少设定上的遗憾在项目里已经不太好改,希望我能重头再写一套图形框架,对整体架构有更深的理解.拭目以待看看我能走到哪一步