PCL - MLS代碼研讀（四）- 偏微分計算函數

前言

本篇延續PCL - MLS代碼研讀（二），繼續介紹偏微分計算函數。

getPolynomialPartialDerivative

還記得在PCL - MLS代碼研讀（二）提到過的PolynomialPartialDerivative結構體，它被用來儲存MLS多項式的偏微分。

pcl::MLSResult::PolynomialPartialDerivative
pcl::MLSResult::getPolynomialPartialDerivative (const double u, const double v) const
{
  // Compute the displacement along the normal using the fitted polynomial
  // and compute the partial derivatives needed for estimating the normal
  PolynomialPartialDerivative d{};
  Eigen::VectorXd u_pow (order + 2), v_pow (order + 2);
  int j = 0;

  //如果階數為2
  //d.z = c_vec[0] + v*c_vec[1] + v**2*c_vec[2] + u*c_vec[3] + u*v*c_vec[4] + u**2*c_vec[5]
  //d.z_u = c_vec[3] + v*c_vec[4] + 2*u*c_vec[5]
  //d.z_uv = c_vec[4]
  d.z = d.z_u = d.z_v = d.z_uu = d.z_vv = d.z_uv = 0;
  //u_pow (ui): pow(u, ui)
  //v_pow (vi): pow(v, vi)
  u_pow (0) = v_pow (0) = 1;
  for (int ui = 0; ui <= order; ++ui)
  {
    for (int vi = 0; vi <= order - ui; ++vi)
    {
      //j: (order-0) + (order-1) + ... + (order-(ui-1)) + vi
      // Compute displacement along normal
      d.z += u_pow (ui) * v_pow (vi) * c_vec[j];

      // Compute partial derivatives
      //與d.z相比:u_pow(ui)對u微分變成ui * u_pow (ui - 1)
      if (ui >= 1)
        d.z_u += c_vec[j] * ui * u_pow (ui - 1) * v_pow (vi);

      //與d.z相比:v_pow(vi)對v微分變成vi * v_pow(vi-1)
      if (vi >= 1)
        d.z_v += c_vec[j] * vi * u_pow (ui) * v_pow (vi - 1);

      //與d.z相比:u_pow (ui) * v_pow (vi)變成ui * u_pow (ui - 1) * vi * v_pow (vi - 1)
      if (ui >= 1 && vi >= 1)
        d.z_uv += c_vec[j] * ui * u_pow (ui - 1) * vi * v_pow (vi - 1);

      //與d.z相比:u_pow(ui)對u微分兩次變成ui * (ui - 1) * u_pow (ui - 2)
      if (ui >= 2)
        d.z_uu += c_vec[j] * ui * (ui - 1) * u_pow (ui - 2) * v_pow (vi);

      //與d.z相比:v_pow(vi)對v微分變成vi * (vi - 1) * v_pow (vi - 2)
      if (vi >= 2)
        d.z_vv += c_vec[j] * vi * (vi - 1) * u_pow (ui) * v_pow (vi - 2);

      //填滿v的次方的快取v_pow
      /**
       * 因為在外圈第一次迭代時(ui==0),內圈的上限最大(order-ui => order)
       * 在外圈>1次迭代,內圈的上限會逐漸減小
       * 所以這裡只需要在外圈第一次迭代時(ui == 0)把v_pow算好就好
       **/
      if (ui == 0)
        v_pow (vi + 1) = v_pow (vi) * v;

      ++j;
    }
    //填滿u的次方的快取u_pow
    u_pow (ui + 1) = u_pow (ui) * u;
  }

  return (d);
}

首先注意到第一層for循環的ui是由0到order，而第二層for循環的vi是由0到order-ui，這是因為ui+vi的上限為order。

另外有兩個長度為order+2的Eigen::VectorXfd：u_pow及v_pow，用於緩存u的次方及v的次方。在for循環之前，將第零個元素初始化為1：

u_pow (0) = v_pow (0) = 1;

在內層for循環結束之前，準備好下一次迭代要用的v的次方項：

    for (int vi = 0; vi <= order - ui; ++vi)
    {
      //...

      //填滿v的次方的快取v_pow
      /**
       * 因為在外圈第一次迭代時(ui==0),內圈的上限最大(order-ui => order)
       * 在外圈>1次迭代,內圈的上限會逐漸減小
       * 所以這裡只需要在外圈第一次迭代時(ui == 0)把v_pow算好就好
       **/
      if (ui == 0)
        v_pow (vi + 1) = v_pow (vi) * v;

      ++j;
    }

注意到這裡只有在ui==0時才會更新v_pow，這是因為在第一次迭代時，內圈循環的上限最大（order-ui = order -0 = order），所以只需要在第一次迭代時把v_pow算過一次即可。

在外層for循環結束之前，也準備好下一次迭代要用的u的次方項：

  for (int ui = 0; ui <= order; ++ui)
  {
    for (int vi = 0; vi <= order - ui; ++vi)
    {
      //...
    }
    //填滿u的次方的快取u_pow
    u_pow (ui + 1) = u_pow (ui) * u;
  }

在for循環裡面，主要是計算$z,\frac{\partial z}{\partial u},\frac{\partial z}{\partial v},\frac{\partial z}{\partial u\partial u},\frac{\partial z}{\partial v\partial v}$和$\frac{\partial z}{\partial u\partial v}$。

在曲面多項式階數為2的情況下，他們的數學式子如下：

$z=a+b\ast v+c\ast v^2+d\ast u+e\ast u\ast v+f\ast u^2$

$\frac{\partial z}{\partial u}=d+e\ast v+2\ast f\ast u$

$\frac{\partial z}{\partial v}=b+2\ast c\ast v+e\ast u$

$\frac{\partial z}{\partial v\partial v}=2\ast c$

$\frac{\partial z}{\partial u\partial u}=2\ast f$

$\frac{\partial z}{\partial u\partial v}=e$

來看看代碼中$z$的計算：

      //j: (order-0) + (order-1) + ... + (order-(ui-1)) + vi
      // Compute displacement along normal
      d.z += u_pow (ui) * v_pow (vi) * c_vec[j];

每一個項次就是簡單地將係數c_vec[j]乘上$u^{ui}$和$v^{vi}$，而$z$就是各項的總和。

再來看看$\frac{\partial z}{\partial u\partial v}$的計算：

      //與d.z相比:u_pow (ui) * v_pow (vi)變成ui * u_pow (ui - 1) * vi * v_pow (vi - 1)
      if (ui >= 1 && vi >= 1)
        d.z_uv += c_vec[j] * ui * u_pow (ui - 1) * vi * v_pow (vi - 1);

對照$z$的代碼，發現就是將u_pow(ui)替換成ui * u_pow (ui - 1)，將v_pow(vi)替換成vi * v_pow(vi - 1)，僅此而已。

其餘的部分，將數學式與代碼對照，應該可以很容易地看懂。