3维曲面插值，matlab C++ 相关书籍

编写语言：C++；
实现效果：类似于matlab中的gridata()函数，或者拟合的曲面穿过采集点的曲面插值算法，可以提供相关书籍论文等资料优先，最好是使用C++语言编写的书籍；

https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/

你可以看看下面这个开源项目
Scattered data interpolation with multilevel B-Splines
https://github.com/ddemidov/mba

https://download.csdn.net/download/stea9527/10535956?spm=1005.2026.3001.5635&utm_medium=distribute.pc_relevant_ask_down.none-task-download-2~default~OPENSEARCH~Rate-4.pc_feed_download_top3ask&depth_1-utm_source=distribute.pc_relevant_ask_down.none-task-download-2~default~OPENSEARCH~Rate-4.pc_feed_download_top3ask

如果你想要matlab中的gridata()功能一模一样。
你先引用matlab中的 Octave 库和任何 Octave 支持库，
直接调用matlab内置函数，builtin-defun-decls.h有函数名称。

#include <iostream>
#include <octave/oct.h>
#include <octave/octave.h>
#include <octave/parse.h>
#include <octave/interpreter.h>
#include <octave/builtin-defun-decls.h>

或者使用点云库(http://www.pointclouds.org)和VoxelGrid

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;

//创建对象
pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_filtered(new pcl::PointCloud<pcl::PointXYZ>());

// 创建过滤对象和进程
pcl::VoxelGrid<pcl::PointXYZ> sor;
sor.setInputCloud (basic_cloud_ptr);
//在这里设置大小(dx, dy, dz)
sor.setLeafSize (0.1, 0.1, 1000);
sor.filter (*cloud_filtered);

下列示例不一定能实现你要的功能

示例1：
void Curve::readIn(GridData& gdata, Triple** data, unsigned int columns, unsigned int rows)
{
    gdata.setSize(columns,rows);
    
    ParallelEpiped range(Triple(DBL_MAX,DBL_MAX,DBL_MAX),Triple(-DBL_MAX,-DBL_MAX,-DBL_MAX));

    /* fill out the vertex array for the mesh. */
    for (unsigned i = 0; i != columns; ++i){
        for (unsigned j = 0; j != rows; ++j){
            gdata.vertices[i][j][0] = data[i][j].x; 
            gdata.vertices[i][j][1] = data[i][j].y;
            gdata.vertices[i][j][2] = data[i][j].z;

            if (data[i][j].x > range.maxVertex.x)
                range.maxVertex.x = data[i][j].x;
            if (data[i][j].y > range.maxVertex.y)
                range.maxVertex.y = data[i][j].y;
            if (data[i][j].z > range.maxVertex.z)
                range.maxVertex.z = data[i][j].z;
            if (data[i][j].x < range.minVertex.x)
                range.minVertex.x = data[i][j].x;
            if (data[i][j].y < range.minVertex.y)
                range.minVertex.y = data[i][j].y;
            if (data[i][j].z < range.minVertex.z)
                range.minVertex.z = data[i][j].z;
         }
    }
    gdata.setHull(range);
    emit readInFinished(title()->string());
}
示例2：
void Curve::sewPeriodic(GridData& gdata)
{
    // sewing
    Triple n;

    unsigned int columns = gdata.columns();
    unsigned int rows = gdata.rows();

    if (gdata.uperiodic()) {
        for (unsigned i = 0; i != columns; ++i) {
            n = Triple(gdata.normals[i][0][0] + gdata.normals[i][rows-1][0],
                gdata.normals[i][0][1] + gdata.normals[i][rows-1][1],
                gdata.normals[i][0][2] + gdata.normals[i][rows-1][2]);

            n.normalize();        
            gdata.normals[i][0][0] = gdata.normals[i][rows-1][0] = n.x;
            gdata.normals[i][0][1] = gdata.normals[i][rows-1][1] = n.y;
            gdata.normals[i][0][2] = gdata.normals[i][rows-1][2] = n.z;
        }
    }

    if (gdata.vperiodic()) {
        for (unsigned j = 0; j != rows; ++j) {
            n = Triple(gdata.normals[0][j][0] + gdata.normals[columns-1][j][0],
                gdata.normals[0][j][1] + gdata.normals[columns-1][j][1],
                gdata.normals[0][j][2] + gdata.normals[columns-1][j][2]);

            n.normalize();        
            gdata.normals[0][j][0] = gdata.normals[columns-1][j][0] = n.x;
            gdata.normals[0][j][1] = gdata.normals[columns-1][j][1] = n.y;
            gdata.normals[0][j][2] = gdata.normals[columns-1][j][2] = n.z;
        }
    }
}
示例3：
void GridPlot::sewPeriodic(GridData& gdata)
{
    // sewing

    Triple n;

    unsigned int columns = gdata.columns();
    unsigned int rows = gdata.rows();

    if (gdata.uperiodic())
    {
        for (unsigned i = 0; i != columns; ++i)
        {
            n = gdata.normals[i][0] + gdata.normals[i][rows-1];
            n.normalize();
            gdata.normals[i][0] = gdata.normals[i][rows-1] = n;
        }
    }
    if (gdata.vperiodic())
    {
        for (unsigned j = 0; j != rows; ++j)
        {
            n = gdata.normals[0][j] + gdata.normals[columns-1][j];
            n.normalize();
            gdata.normals[0][j] = gdata.normals[columns-1][j] = n;
        }
    }
}
示例4：
void writeVolume(GridData &s, AABB &bbox, int channels, Stream *stream) {
        stream->writeChar('V');
        stream->writeChar('O');
        stream->writeChar('L');
        stream->writeChar(3);

        int type = 1;
        stream->writeInt(type);

        stream->writeInt(s.size());
        stream->writeInt(s[0].size());
        stream->writeInt(s[0][0].size());

        stream->writeInt(channels);

        stream->writeSingle(bbox.min.x);
        stream->writeSingle(bbox.min.y);
        stream->writeSingle(bbox.min.z);
        stream->writeSingle(bbox.max.x);
        stream->writeSingle(bbox.max.y);
        stream->writeSingle(bbox.max.z);

        float *data = new float[s.size() * s[0].size() * s[0][0].size() * channels];
        for (int z = 0; z < s[0][0].size(); z++) {
            for (int y = 0; y < s[0].size(); y++) {
                for (int x = 0; x < s.size(); x++) {
                    for (int c = 0; c < channels; c++) {
                        data[((z * s[0].size() + y) * s.size() + x) * channels + c] = s[x][y][z][c];
                    }
                }
            }
        }
        stream->writeSingleArray(data, s.size() * s[0].size() * s[0][0].size() * channels);
        delete[] data;
    }
示例5：
void GridPlot::readIn(GridData& gdata, Triple** data, unsigned int columns, unsigned int rows)
{
    gdata.setSize(columns,rows);

    ParallelEpiped range(Triple(DBL_MAX,DBL_MAX,DBL_MAX),Triple(-DBL_MAX,-DBL_MAX,-DBL_MAX));

    /* fill out the vertex array for the mesh. */
    for (unsigned i = 0; i != columns; ++i)
    {
        for (unsigned j = 0; j != rows; ++j)
        {
            Triple& val = data[i][j];

            gdata.vertices[i][j] = val;

            if (val.x > range.maxVertex.x)
                range.maxVertex.x = val.x;
            if (val.y > range.maxVertex.y)
                range.maxVertex.y = val.y;
            if (val.z > range.maxVertex.z)
                range.maxVertex.z = val.z;
            if (val.x < range.minVertex.x)
                range.minVertex.x = val.x;
            if (val.y < range.minVertex.y)
                range.minVertex.y = val.y;
            if (val.z < range.minVertex.z)
                range.minVertex.z = val.z;
        }
    }
    gdata.setHull(range);
}
示例6：
void Curve::readIn(GridData& gdata, double** data, unsigned int columns, unsigned int rows,
                   double minx, double maxx, double miny, double maxy)
{
    gdata.setPeriodic(false,false);
    gdata.setSize(columns, rows);
    
    double dx = (maxx - minx) / (gdata.columns() - 1);
    double dy = (maxy - miny) / (gdata.rows() - 1);

    double tmin = DBL_MAX;
    double tmax = -DBL_MAX;

    /* fill out the vertex array for the mesh. */
    for (unsigned i = 0; i != columns; ++i){
        for (unsigned j = 0; j != rows; ++j){
            gdata.vertices[i][j][0] = minx + i*dx;
            gdata.vertices[i][j][1] = miny + j*dy;
            gdata.vertices[i][j][2] = data[i][j];

            double val = data[i][j];
            if (val > tmax)
                tmax = val;
            if (val < tmin)
                tmin = val;
         }
    }

    ParallelEpiped hull = ParallelEpiped(Triple(gdata.vertices[0][0][0], gdata.vertices[0][0][1], tmin), 
                                        Triple(gdata.vertices[gdata.columns() - 1][gdata.rows() - 1][0],
                                               gdata.vertices[gdata.columns() - 1][gdata.rows() - 1][1], tmax));

    gdata.setHull(hull);
    emit readInFinished(title()->string());
}
示例7：
void GridPlot::calcNormals(GridData& gdata)
{
    unsigned int rows = gdata.rows();
    unsigned int columns = gdata.columns();

    // normals
    
    Triple u, v, n;  // for cross product
    Triple vert_ip1j;
    Triple vert_ijp1;
    Triple vert_im1j;
    Triple vert_ijm1;

    for (unsigned i = 0; i < columns; ++i)
    {
        for (unsigned j = 0; j < rows; ++j)
        {
            Triple& n = gdata.normals[i][j];
            n = Triple(0,0,0);

            const Triple& vert_ij = gdata.vertices[i][j];

            if (i<columns-1 && j<rows-1)
            {
                vert_ip1j = gdata.vertices[i+1][j] - vert_ij;
                vert_ijp1 = gdata.vertices[i][j+1] - vert_ij;
                n += normalizedcross(vert_ip1j, vert_ijp1); // right hand system here !
            }

            if (i>0 && j>0)
            {
                vert_im1j = gdata.vertices[i-1][j] - vert_ij;
                vert_ijm1 = gdata.vertices[i][j-1] - vert_ij;
                n += normalizedcross(vert_im1j, vert_ijm1);
            }

            if (i>0 && j<rows-1)
            {
                n += normalizedcross(vert_ijp1, vert_im1j);
            }

            if (i<columns-1 && j>0)
            {
                n += normalizedcross(vert_ijm1, vert_ip1j);
            }

            n.normalize();
        }
    }
}
示例8：
void initS(GridData &s, const Vector3i &res) {
        s.resize(res.x);
        for (int i = 0; i < res.x; i++) {
            s[i].resize(res.y);
            for (int j = 0; j < res.y; j++) {
                s[i][j].resize(res.z);
            }
        }
    }
示例9：
void GridPlot::readIn(GridData& gdata, double** data, unsigned int columns, unsigned int rows
                      , double minx, double maxx, double miny, double maxy)
{
    gdata.setPeriodic(false,false);
    gdata.setSize(columns,rows);

    double dx = (maxx - minx) / (gdata.columns() - 1);
    double dy = (maxy - miny) / (gdata.rows() - 1);

    double tmin = DBL_MAX;
    double tmax = -DBL_MAX;

    /* fill out the vertex array for the mesh. */
    for (unsigned i = 0; i != columns; ++i)
    {
        for (unsigned j = 0; j != rows; ++j)
        {
            Triple& gdata_ij = gdata.vertices[i][j];
            double& val = data[i][j];

            gdata_ij.x = minx + i*dx;
            gdata_ij.y = miny + j*dy;
            gdata_ij.z = val;

            if (val > tmax)
                tmax = val;
            if (val < tmin)
                tmin = val;
        }
    }
    ParallelEpiped hull =
            ParallelEpiped(
                Triple(
                    gdata.vertices[0][0].x,
            gdata.vertices[0][0].y,
            tmin
            ),
            Triple(
                gdata.vertices[gdata.columns()-1][gdata.rows()-1].x,
            gdata.vertices[gdata.columns()-1][gdata.rows()-1].y,
            tmax
            )
            );

    gdata.setHull(hull);
}
示例10：
void Curve::calcNormals(GridData& gdata)
{
    unsigned int rows = gdata.rows();
    unsigned int columns = gdata.columns();
  
    // normals
    Triple u, v, n;  // for cross product

    for (unsigned i = 0; i != columns; ++i) {
        for (unsigned j = 0; j != rows; ++j) {
            n = Triple(0,0,0);

            if (i<columns-1 && j<rows-1) {
                /*    get two vectors to cross */
                u = Triple(gdata.vertices[i+1][j][0] - gdata.vertices[i][j][0],
                    gdata.vertices[i+1][j][1] - gdata.vertices[i][j][1],
                    gdata.vertices[i+1][j][2] - gdata.vertices[i][j][2]);

                v = Triple(
                gdata.vertices[i][j+1][0] - gdata.vertices[i][j][0],
                gdata.vertices[i][j+1][1] - gdata.vertices[i][j][1],
                gdata.vertices[i][j+1][2] - gdata.vertices[i][j][2]);

                /* get the normalized cross product */ 
                n += normalizedcross(u,v); // right hand system here !
            }

            if (i>0 && j<rows-1) {
                u = Triple(gdata.vertices[i][j+1][0] - gdata.vertices[i][j][0],
                    gdata.vertices[i][j+1][1] - gdata.vertices[i][j][1],
                    gdata.vertices[i][j+1][2] - gdata.vertices[i][j][2]);
              
                v = Triple(gdata.vertices[i-1][j][0] - gdata.vertices[i][j][0],
                    gdata.vertices[i-1][j][1] - gdata.vertices[i][j][1],
                    gdata.vertices[i-1][j][2] - gdata.vertices[i][j][2]);

                n += normalizedcross(u,v); 
            }

            if (i>0 && j>0) {
                u = Triple(gdata.vertices[i-1][j][0] - gdata.vertices[i][j][0],
                    gdata.vertices[i-1][j][1] - gdata.vertices[i][j][1],
                    gdata.vertices[i-1][j][2] - gdata.vertices[i][j][2]);

                v = Triple(gdata.vertices[i][j-1][0] - gdata.vertices[i][j][0],
                    gdata.vertices[i][j-1][1] - gdata.vertices[i][j][1],
                    gdata.vertices[i][j-1][2] - gdata.vertices[i][j][2]);

                n += normalizedcross(u,v);
            }

            if (i<columns-1 && j>0) {
                u = Triple(gdata.vertices[i][j-1][0] - gdata.vertices[i][j][0],
                    gdata.vertices[i][j-1][1] - gdata.vertices[i][j][1],
                    gdata.vertices[i][j-1][2] - gdata.vertices[i][j][2]);

                v = Triple(gdata.vertices[i+1][j][0] - gdata.vertices[i][j][0],
                    gdata.vertices[i+1][j][1] - gdata.vertices[i][j][1],
                    gdata.vertices[i+1][j][2] - gdata.vertices[i][j][2]);

                n += normalizedcross(u,v);
            }

            n.normalize();
            gdata.normals[i][j][0] = n.x;
            gdata.normals[i][j][1] = n.y;
            gdata.normals[i][j][2] = n.z;
        } 
    } 
}
示例11：
int main(int argc, char *argv[]) 
{
    char *f;
    int dpr = 10;
    GridNMC *gridbase;
    GridList *list;
    gridbase = new GridNMC();
    list = new GridList(gridbase);

    f = argv[1];

// First load the data
    fprintf (stderr, "Archivo %s  Datos %p %d\n", f, gridbase, gridbase->getN()); fflush(stderr);
    list->setList(f);

    if (list->size()<=0)
        return 1;

    fprintf (stderr, "Lista %d\n", list->size());
    

    printf ("latitude = ");
    for (int i=0; i < gridbase->getN(); i++) {
        printf("%g", gridbase->la(i));
        if (i < gridbase->getN()-1)
            printf(", ");
        else 
            printf(" ;\n");
        if (i > 0 && i % dpr==0)
            printf("\n");
    }

    printf ("longitude = ");
    for (int i=0; i < gridbase->getN(); i++) {
        printf("%g", gridbase->lo(i));
        if (i < gridbase->getN()-1)
            printf(", ");
        else 
            printf(" ;\n");
        if (i > 0 && i % dpr==0)
            printf("\n");
    }

    for (int j=0; j < list->size(); j++) 
    {
        GridData *d = list->getData(j);
        fprintf(stderr, "Data %d %g %g %p\n", j, d->getMin(), d->getMax(), d->getData());  
      
        printf ("temperature =\n");
        for (int i=0; i < gridbase->getN(); i++) {
            printf("%g", d->getData(i));
            if (i < gridbase->getN()-1)
                printf(", ");
            else 
                printf(" ;\n");
            if (i > 0 && i % dpr==0)
                printf("\n");
        }
    }    

    return 0;
}

```c++