Matlab griddata équivalent en C++

Question

Je cherche un C++ équivalent à la fonction Gridata de Matlab, ou n'importe quelle méthode d'interpolation globale 2D.

J'ai un code C++ qui utilise Eigen 3. Je vais avoir un vecteur propre qui contiendra les valeurs x, y et z, et deux matrices propres équivalentes à celles produites par Meshgrid dans Matlab. Je voudrais interpoler les valeurs z des vecteurs sur les points de la grille définis par les équivalents Meshgrid (qui dépasseront un peu l'extérieur des points d'origine, donc une extrapolation mineure est requise).

Je ne suis pas trop gêné par la précision - il n'a pas besoin d'être parfait. Cependant, je ne peux pas accepter NaN comme une solution - l'interpolation doit être calculée partout sur le maillage, indépendamment des lacunes de données. En d'autres termes, rester à l'intérieur de la coque convexe n'est pas une option.

Je préférerais ne pas écrire une interpolation à partir de zéro, mais si quelqu'un veut me montrer une recette assez bonne (et explicite), je vais essayer. Ce n'est pas la chose la plus haineuse à écrire (du moins dans un sens algorithmique), mais je ne veux pas réinventer la roue.

Effectivement ce que j'ai sont des emplacements de terrain dispersés, et je souhaite définir un maillage rectiligne qui suit nominalement une certaine distance sous la topographie pour une utilisation ultérieure. Une fois que j'ai les points de noeud, je serai bon.

Mes recherches jusqu'à présent:

La question posée ici: MATLAB functions in C++ a produit une réponse proche, mais malheureusement, la suggestion n'a pas été libre (SciMath).

J'ai essayé de comprendre la fonction d'interpolation utilisée dans Generic Mapping Tools, et j'ai été récompensé par un mal de tête. J'ai brièvement regardé dans la bibliothèque Algorithmes de la Grille (GRAL). Si quelqu'un a un commentaire, je l'apprécierais. Eigen a un paquet d'interpolation non supporté, mais il semble être juste pour des courbes (pas des surfaces).

Éditer: VTK a une fonctionnalité de matplotlib. Vraisemblablement, il doit y avoir une interpolation utilisée quelque part à des fins d'affichage. Est-ce que quelqu'un sait si c'est accessible et utilisable?

Merci.

Answer 1

Ceci est probablement un peu en retard, mais j'espère que cela aide quelqu'un.

Méthode 1.) Octave: Si vous venez de Matlab, l'une des façons consiste à incorporer le clone Matlab Octave directement dans le programme C++. Je n'ai pas beaucoup d'expérience, mais vous pouvez appeler les fonctions de la bibliothèque d'octave directement à partir d'un fichier cpp.

Voir ici, par exemple. http://www.gnu.org/software/octave/doc/interpreter/Standalone-Programs.html#Standalone-Programs

griddata est inclus dans le package de géométrie d'octave.

Méthode 2.) PCL: La façon dont je le fais est d'utiliser la bibliothèque de nuages ​​de points (http://www.pointclouds.org) et VoxelGrid. Vous pouvez définir les tailles de bacs x et y comme vous le souhaitez, puis définir une taille de bac très grande, qui vous donne une valeur z pour chaque bac x, y. Le hic est que les valeurs x, y et z sont le centroïde pour les points moyennés dans la poubelle, pas les centres de bin (ce qui est aussi pourquoi cela fonctionne pour cela).Vous devez donc masser les valeurs x, y lorsque vous avez terminé:

Ex: // lire dans une liste de valeurs séparées par des virgules (x, y, z) FICHIER * fp; fp = fopen ("points.xyz", "r");

//store them in PCL's point cloud format 
pcl::PointCloud<pcl::PointXYZ>::Ptr basic_cloud_ptr (new pcl::PointCloud<pcl::PointXYZ>); 

int numpts=0; 
double x,y,z; 
while(fscanf(fp, "%lg, %lg, %lg", &x, &y, &z)!=EOF) 
{ 
pcl::PointXYZ basic_point; 
basic_point.x = x; basic_point.y = y; basic_point.z = z; 
basic_cloud_ptr->points.push_back(basic_point); 
} 
fclose(fp); 
basic_cloud_ptr->width = (int) basic_cloud_ptr->points.size(); 
basic_cloud_ptr->height = 1; 

// create object for result 
pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_filtered(new pcl::PointCloud<pcl::PointXYZ>()); 

// create filtering object and process 
pcl::VoxelGrid<pcl::PointXYZ> sor; 
sor.setInputCloud (basic_cloud_ptr); 
//set the bin sizes here. (dx,dy,dz). for 2d results, make one of the bins larger 
//than the data set span in that axis 
sor.setLeafSize (0.1, 0.1, 1000); 
sor.filter (*cloud_filtered); 

Alors que cloud_filtered est maintenant un nuage de points qui contient un point pour chaque casier. Ensuite, je fais juste une matrice à 2 d et je passe par le nuage de points en assignant des points à leurs cases x, y si je veux une image, etc., comme cela serait produit par griddata. Cela fonctionne plutôt bien, et c'est beaucoup plus rapide que la grille de matlab pour les grands ensembles de données.