将以下代码改为C++代码： import scipy.special as sp import numpy as np import numba from numba import njit,prange import math import trimesh as tri fileName="data/blub.obj" outName='./output/blub_rec.obj' # 参数 # 限制选取球谐基函数的带宽 bw=64 # 极坐标,经度0<=theta<2*pi,纬度0<=phi<pi; # (x,y,z)=r(sin(phi)cos(theta),sin(phi)sin(theta),cos(phi)) def get_angles(x,y,z): r=np.sqrt(x*x+y*y+z*z) x/=r y/=r z/=r phi=np.arccos(z) if phi==0: theta=0 theta=np.arccos(x/np.sin(phi)) if y/np.sin(phi)<0: theta+=math.pi return [theta,phi] if __name__=='__main__': # 载入网格 mesh=tri.load(fileName) # 获得网格顶点(x,y,z)对应的(theta,phi) numV=len(mesh.vertices) angles=np.zeros([numV,2]) for i in range(len(mesh.vertices)): v=mesh.vertices[i] [angles[i,0],angles[i,1]]=get_angles(v[0],v[1],v[2]) # 求解方程:x(theta,phi)=对m,l求和 a^m_lY^m_l(theta,phi) 解出系数a^m_l # 得到每个theta,phi对应的x X,Y,Z=np.zeros([numV,1]),np.zeros([numV,1]),np.zeros([numV,1]) for i in range(len(mesh.vertices)): X[i],Y[i],Z[i]=mesh.vertices[i,0],mesh.vertices[i,1],mesh.vertices[i,2] # 求出Y^m_l(theta,phi)作为矩阵系数 sph_harm_values=np.zeros([numV,(bw+1)*(bw+1)]) for i in range(numV): for l in range(bw): for m in range(-l,l+1): sph_harm_values[i,l*(l+1)+m]=sp.sph_harm(m,l,angles[i,0],angles[i,1]) print('系数矩阵维数:{}'.format(sph_harm_values.shape)) # 求解方程组,得到球谐分解系数 a_x=np.linalg.lstsq(sph_harm_values,X,rcond=None)[0] a_y=np.linalg.lstsq(sph_harm_values,Y,rcond=None)[0] a_z=np.linalg.lstsq(sph_harm_values,Z,rcond=None)[0] # 从系数恢复的x,y,z坐标,存为新的点云用于比较 x=np.matmul(sph_harm_values,a_x) y=np.matmul(sph_harm_values,a_y) z=np.matmul(sph_harm_values,a_z) with open(outName,'w') as output: for i in range(len(x)): output.write("v %f %f %f\n"%(x[i,0],y[i,0],z[i,0]))

#include <iostream>
#include <fstream>
#include <cmath>
#include <vector>
#include <Eigen/Dense>
#include <unsupported/Eigen/Splines>
#include <trimesh/trimesh.h>
using namespace Eigen;

std::string fileName = "data/blub.obj";
std::string outName = "./output/blub_rec.obj";
int bw = 64;

std::vector<double> get_angles(double x, double y, double z){
double r = std::sqrt(x*x + y*y + z*z);
double phi = std::acos(z/r);
double theta;
if(phi == 0){
theta = 0;
}else{
theta = std::acos(x/std::sin(phi));
if(y/std::sin(phi) < 0){
theta += M_PI;
}
}
return {theta, phi};
}

int main(){
// Load mesh
trimesh::TriMesh mesh;
mesh.read(fileName);

// Get (theta, phi) for each vertex
int numV = mesh.vertices.size();
std::vector<std::vector<double>> angles(numV, std::vector<double>(2));
for(int i=0; i<numV; ++i){
auto v = mesh.vertices[i];
angles[i] = get_angles(v[0], v[1], v[2]);
}

// Solve equation: x(theta, phi) = sum(a^m_l * Y^m_l(theta, phi)), and get coefficients a^m_l
// Get x, y, z for each (theta, phi)
std::vector<double> X(numV), Y(numV), Z(numV);
for(int i=0; i<numV; ++i){
X[i] = mesh.vertices[i][0];
Y[i] = mesh.vertices[i][1];
Z[i] = mesh.vertices[i][2];
}

// Get sph_harm_values as matrix coefficients
MatrixXd sph_harm_values(numV, (bw+1)*(bw+1));
for(int i=0; i<numV; ++i){
for(int l=0; l<bw; ++l){
for(int m=-l; m<=l; ++m){
sph_harm_values(i, l*(l+1)+m) = std::real(std::pow(-1, m) * boost::math::spherical_harmonic_r(l, m, angles[i][1], angles[i][0]));
}
}
}

// Solve the equation system to get the spherical harmonic decomposition coefficients
VectorXd a_x = sph_harm_values.colPivHouseholderQr().solve(X);
VectorXd a_y = sph_harm_values.colPivHouseholderQr().solve(Y);
VectorXd a_z = sph_harm_values.colPivHouseholderQr().solve(Z);

// Generate new point clouds from the coefficients
std::vector<double> new_X(numV), new_Y(numV), new_Z(numV);
for(int i=0; i<numV; ++i){
for(int l=0; l<bw; ++l){
for(int m=-l; m<=l; ++m){
new_X[i] += std::real(std::pow(-1, m) * a_x(l*(l+1)+m) * boost::math::spherical_harmonic_r(l, m, angles[i][1], angles[i][0]));
new_Y[i] += std::real(std::pow(-1, m) * a_y(l*(l+1)+m) * boost::math::spherical_harmonic_r(l, m, angles[i][1], angles[i][0]));
new_Z[i] += std::real(std::pow(-1, m) * a_z(l*(l+1)+m) * boost::math::spherical_harmonic_r(l, m, angles[i][1], angles[i][0]));
}
}
}

// Write new point clouds to obj file
std::ofstream output(outName);
for(int i=0; i<numV; ++i){
output << "v " << new_X[i] << " " << new_Y[i] << " " << new_Z[i] << std::endl;
}
output.close();

return 0;
}