之前我写了一篇关于针孔相机和鱼眼相机从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的图形框架库.自我感觉做的还可以,但是还是有不少设定上的遗憾在项目里已经不太好改,希望我能重头再写一套图形框架,对整体架构有更深的理解.拭目以待看看我能走到哪一步

这段时间在优化3d场景坐标(x,y,z)映射到相机的2d图片(u,v)的代码.在之前的实现里,我这边前端实现了无畸变的坐标转换,但是在一些畸变比较大的图片里,坐标点显示偏差较大,功能无法使用.

之前我们标注平台的3d立方体映射到相机图片上是通过后端计算,对需求的实时性上面有很大的影响.于是近期我在尝试前端去实现相机畸变的矫正研究.

https://zhuanlan.zhihu.com/p/94445955
首先我网上搜索相机畸变方面文章,主要参考了知乎,然后咨询了一下作者本人,

大概理解了坐标转换的原理和顺序:
已知世界坐标(x,y,z)
1.根据相机矩阵转换为相机坐标(x相,y相,z相) z轴正方向为相机朝向,x轴向右,y轴向下
2.坐标归一化得到 (x归,y归,1) = (x相/z相,y相/z相,1)
3.得到距离 r = sqrt(x归^2+y归^2)
4.这时候进行相机畸变调整 得到x',y'

x' = x归*(1+k1*r^2+k2*r^4+k3*r^6)+2*p1*x归*y归+p2*(r^2+2*x归^2)
y' = y归*(1+k1*r^2+k2*r^4+k3*r^6)+2*p2*x归*y归+p1*(r^2+2*y归^2)

5.(x',y',1)乘以相机内参得到最终(u,v)值

开发完成后,发现有一部分照片的畸变调整还是不行,反而更离谱
仔细调试代码+借阅后台代码进行查看

发现后台是根据不同相机做不同处理的,于是再次询问了上面知乎的作者Mver,确定了鱼眼相机的处理逻辑和普通针孔相机处理逻辑不一致,又翻阅了关于鱼眼相机畸变处理的资料
https://zhuanlan.zhihu.com/p/577186117
https://blog.csdn.net/qq_39642978/article/details/112758040
得到结论

已知世界坐标(x,y,z)

转换为相机的相对坐标
相机坐标归一化(x相/z相,y相/z相,1) => (x归,y归,1)
计算 r = sqrt(x归^2+y归^2)
计算 theta = atan(r)
计算矫正 r_d = theta(1+k1theta^2+k2theta^4+k3theta^6+k4*theta^8)
计算 scale = r_d / r
计算 x' = x归 scale ; y' = y归 scale
(x',y',1) 乘以相机内参得到图片像素(u,v)

其中做到第5布计算矫正r_d的时候,卡住了我一晚上没成功,最终通过另外一个博客得到正解(上面博客介绍应该也是对的,只不过参数值和opencv不同)
https://blog.csdn.net/qq_16137569/article/details/112398976

到这步总算所有畸变都处理成功,最终测试成功时,那一刻感到从心底升起的喜悦,只感觉神经气爽,海阔天空
我想这就是高级多巴胺吧,是让我能在求知的路上一路前进的动力来源.O(∩_∩)O
在这里感谢Mver的帮助,虽然他肯定不知道