code: def estimate_shape(
      x, shapeMU, shapePC, shapeEV, expression, s, R, t2d, lamb=3000
  ):
      """
      Args:
          x: (2, n). image points (to be fitted)
          shapeMU: (3n, 1)
          shapePC: (3n, n_sp)
          shapeEV: (n_sp, 1)
          expression: (3, n)
          s: scale
          R: (3, 3). rotation matrix
          t2d: (2,). 2d translation
          lambda: regulation coefficient
  
      Returns:
          shape_para: (n_sp, 1) shape parameters(coefficients)
      """
      x = x.copy()
      assert shapeMU.shape[0] == shapePC.shape[0]
      assert shapeMU.shape[0] == x.shape[1] * 3
  
      dof = shapePC.shape[1]
  
      n = x.shape[1]
      sigma = shapeEV
      t2d = np.array(t2d)
      P = np.array([[1, 0, 0], [0, 1, 0]], dtype=np.float32)
      A = s * P.dot(R)
  
      # --- calc pc
      pc_3d = np.resize(shapePC.T, [dof, n, 3])  # 199 x n x 3
      pc_3d = np.reshape(pc_3d, [dof * n, 3])
      pc_2d = pc_3d.dot(A.T.copy())  # 199 x n x 2
  
      pc = np.reshape(pc_2d, [dof, -1]).T  # 2n x 199
  
      # --- calc b
      # shapeMU
      mu_3d = np.resize(shapeMU, [n, 3]).T  # 3 x n
      # expression
      exp_3d = expression
      #
      b = A.dot(mu_3d + exp_3d) + np.tile(
          t2d[:, np.newaxis], [1, n]
      )  # 2 x n
      b = np.reshape(b.T, [-1, 1])  # 2n x 1
  
      # --- solve
      equation_left = np.dot(pc.T, pc) + lamb * np.diagflat(
          1 / sigma**2
      )
      x = np.reshape(x.T, [-1, 1])
      equation_right = np.dot(pc.T, x - b)
  
      shape_para = np.dot(np.linalg.inv(equation_left), equation_right)
  
      return shape_para
  
imports: import numpy
  from ..None import mesh
label_desc: Use image node_id 1 for calling a global function with example usage: estimate_shape(x, shapeMU, shapePC, shapeEV, expression, s, R, t2d, lamb) and returns: shape_para