How to create a 3D X,Y,Z array from 2D faces so that contiguity between points is preserved

Question

I'm trying to implement an algorithm that from six 2-D matrices of X,Y,Z points store them in a 3D matrix of X,Y,Z points in such a manner that their connectivity is preserved. An example will clarify the problem.

I've to represent with three 3-D matrices

( M3D_X(:,:,:), M3D_Y(:,:,:),M_3D_Z(:,:,:) ) 

an hexahedral region of space. Of that hexahedral region will be known only its faces, stored in three 2-D matrices

( [face1_x(:,:,:),face1_y(:,:,:),face1_z(:,:,:)], [face2_x.....]....,],...,[...,face6_Z(:,:,:)] )

The only information known about connectivity between faces is that face_i have a side in common with face_j. It isn't given which side is in common. So,how to store the 2-D matrices in the 3-D matrix so that contiguity between faces is preserved ? Just to clarify, the interior of the hexahedral region will be created with 3-D transfinite interpolation.

Example ( with only two faces instead of six):

! []-->matrix ,--> column separator ;--> row separator 
! face_i_x= [p11_x, p12_x,..,p1n_x;
             ....................
             pm1_x,..........,pmn_x]

face1_x = [0.0,0.5,1.0;
           0.0,0.5,1.0;
           0.0,0.5,1.0]

face1_y = [0.0,0.0,0.0;
           0.5,0.5,0.5;
           1.0,1.0,1.0]

face1_z = [0.0,0.0,0.0;
           0.0,0.0,0.0;
           0.0,0.0,0.0]


face2_x = [0.0,0.5,1.0;
           0.0,0.5,1.0;
           0.0,0.5,1.0]

face2_y = [1.0,1.0,1.0;
           1.0,1.0,1.0;
           1.0,1.0,1.0]

face2_z = [0.0,0.0,0.0;
           0.5,0.5,0.5;
           1.0,1.0,1.0]

cube3D_x = (:,:,:)
cube3D_y = (:,:,:)
cube3D_z = (:,:,:)

face1 represent the face at z=0 of the unit cube. face2 represents the face at y=1.0 of the unit cube. To retain the simplicity of this example I will not consider face3,..,face6. Now I need to create a 3D matrix to representing the unit cube. I can start like this

cube3D_x(1,:,:)=face1_x
cube3D_y(1,:,:)=face1_y
cube3D_z(1,:,:)=face1_z

but then the second (and subsequent insertion must be handled with special care because face1 has a side in common with face2 (the edge at y=1.0, z=0.0).
How to insert the second face in the cube ?

For more clarity here is a python script for better visualization of the problem:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

face1x=np.array( [ [0.0,0.4,0.6,1.0] , [0.0,0.4,0.6,1.0] , [0.0,0.4,0.6,1.0],[0.0,0.4,0.6,1.0] ])

face1y=np.array( [ [0.0, 0.0, 0.0, 0.0], [0.3,0.3,0.3,0.3], [0.7,0.7,0.7,0.7],[1.0,1.0,1.0,1.0]])

face1z=np.array( [ [0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0],[0.0,0.0,0.0,0.0]] )



face2x=np.array( [ [0.0,0.4,0.6,1.0] , [0.0,0.4,0.6,1.0],[0.0,0.4,0.6,1.0] , [0.0,0.4,0.6,1.0],[0.0,0.4,0.6,1.0] ])

face2y=np.array( [ [0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0],[0.0,0.0,0.0,0.0],[0.0,0.0,0.0,0.0]] )

face2z=np.array( [ [0.0, 0.0, 0.0, 0.0], [0.3,0.3,0.3,0.3],[0.5,0.5,0.5,0.5] ,[0.7,0.7,0.7,0.7],[1.0,1.0,1.0,1.0]])

# because face2 is stored with arbitrary direction of indices, just simulate this randomness with some flipping and rotation

np.flip(face2x)
np.flip(face2y)
np.flip(face2z)

np.rot90(face2x)
np.rot90(face2y)
np.rot90(face2z)

# now plot the two faces
fig = plt.figure()
ax = fig.gca(projection='3d')

ax.plot_wireframe(face1x,face1y,face1z,color='b')

ax.plot_wireframe(face2x,face2y,face2z,color='r')
plt.show()

Also : https://en.wikipedia.org/wiki/Regular_grid I'm not interested in performance. Thank you very much.

----edit--- The same problem can be posed with lines instead of faces and faces instead of 3D volumes.

line1 = [1,2,3]
line2 = [5,4,3]
line3 = [5,6,7]
line4 = [9,8,7]
face= matrix(3,3)

One can start with

face(1,:)= line1 

and he obtains

face = [1,2,3 ;
        *,*,* ;
        *,*,*  ]

now if the line2 is inserted like this

face(:,3)=line2 

the result will be

face = [1,2,5 ;
        *,*,4 ;
        *,*,3 ]

when the correct solution is

face(:,3)= reverse(line2) 

(because point 3 is in common with line1 and line2 ) to obtain

face = [1,2,3 ;
        *,*,4 ;
        *,*,5 ]

Hope this makes the problem clearer.

Answer 1

So If I get it right you got 6x 2D matrices of 3D points representing 6 surfaces of some shape (that connects together) in any flip/mirror/order state between each-other and want to construct single 3D matrix with empty interior (for example point (0,0,0) and its surface will contain the 6 surfaces reordered reoriented so their edges match.

As you always have 6 3D surfaces with 2D topology your 3D topology can be fixed which simplifies things alot.

1. define 3D topology

   It can be any I decided to use this one:

                           -------------
                          /           /
                         /     5     /
                      y /           /------
                       -------------      |
                       x      |           |
                              |     3     |
           /|                 |           |          /|
          / |                 |           |         / |
         /  |               y -------------        /  |
        |   |                 x                   |   |
        | 0 |                                     | 1 |
        |  /                                      |  /
        | /    -------------                      | /
       y|/x    |           |                     y|/x
               |           |
               |     2     |
               |           |------------
               |           |          /
             y -------------   4     /
               x     y /            /
                       -------------
                       x
   

   The numbers 0..5 represent your face (surfaces/faces f0,f1,...f5) index and x,y index origin and direction within each 2D matrix.

2. convert your 6 surfaces to match topology

   You will need 2 arrays of 6 values. One telling if input matrix is already used in topology (us[6]) and the other holding index of input surface mapped to topology f0...f5 (id[6]). At start set all to -1. The us can hold also the reverse of id if needed.

   The first face f0 can be hardcoded as first input surface as is so

   id[0]=0;
   us[0]=0;
   

   Now its just a matter of going through all edges of f0 and find matching edge in all unused yet surfaces. So each surface has 4 edges but each edge have 2 possible directions. so that is 8 edge tests per face.

   Once matching edge is found we need to orient the matching face so it matches our selected topology. So we need to mirror x, mirror y and or swap x/y.

   This will obtain f2,f3,f4,f5.

   Now just take any of the new found faces and find f0 by matching its shared edge. Again orient f0 so it matches topology.

3. copy surfaces into 3D matrix

   so allocate 3D matrix size (obtained from resolutions of f0(x,y),f2(x) and now just copy all the f0..f5 to their location in matrix... You can fill up the interior points with (0,0,0) or interpolate between the surfaces...

Here small C++/VCL/OpenGL example:

arrays.h 1D/2D/3D containers

//---------------------------------------------------------------------------
#ifndef _arrays_h
#define _arrays_h
/*---------------------------------------------------------------------------
// inline safety constructors/destructors add this to any class/struct used as T
T()     {}
T(T& a) { *this=a; }
~T()    {}
T* operator = (const T *a) { *this=*a; return this; }
//T* operator = (const T &a) { ...copy... return this; }
//-------------------------------------------------------------------------*/
template <class T> class array1D
    {
public:
    T       *dat;   // items array
    int     nx;     // allocated size
    array1D()   { dat=NULL; free(); }
    array1D(array1D& a) { *this=a; }
    ~array1D()  { free(); }
    array1D* operator = (const array1D *a) { *this=*a; return this; }
    array1D* operator = (const array1D &a)
        {
        int x;
        allocate(a.nx);
        for (x=0;x<nx;x++) dat[x]=a.dat[x];
        return this;
        }
    T& operator[] (int x){ return dat[x]; }
    void free()
        {
        if (dat!=NULL) delete[] dat;
        nx=0; dat=NULL;
        }
    int  allocate(int _nx)
        {
        free();
        dat=new T[_nx];
        nx=_nx;
        }
    };
//---------------------------------------------------------------------------
template <class T> class array2D
    {
public:
    T       **dat;  // items array
    int     nx,ny;  // allocated size
    array2D()   { dat=NULL; free(); }
    array2D(array2D& a) { *this=a; }
    ~array2D()  { free(); }
    array2D* operator = (const array2D *a) { *this=*a; return this; }
    array2D* operator = (const array2D &a)
        {
        int x,y;
        allocate(a.nx,a.ny);
        for (x=0;x<nx;x++)
         for (y=0;y<ny;y++)
          dat[x][y]=a.dat[x][y];
        return this;
        }
    T* operator[] (int x){ return dat[x]; }
    void free()
        {
        if (dat!=NULL)
            {
            if (dat[0]!=NULL) delete[] dat[0];
            delete[] dat;
            }
        nx=0; ny=0; dat=NULL;
        }
    int  allocate(int _nx,int _ny)
        {
        free();
        dat=new T*[_nx];
        dat[0]=new T[_nx*_ny];
        nx=_nx; ny=_ny;
        for (int x=1;x<nx;x++) dat[x]=dat[x-1]+ny;
        }
    };
//---------------------------------------------------------------------------
template <class T> class array3D
    {
public:
    T       ***dat;     // items array
    int     nx,ny,nz;   // allocated size
    array3D()   { dat=NULL; free(); }
    array3D(array3D& a) { *this=a; }
    ~array3D()  { free(); }
    array3D* operator = (const array3D *a) { *this=*a; return this; }
    array3D* operator = (const array3D &a)
        {
        int x,y,z;
        allocate(a.nx,a.ny,a.nz);
        for (x=0;x<nx;x++)
         for (y=0;y<ny;y++)
          for (z=0;z<nz;z++)
           dat[x][y][z]=a.dat[x][y][z];
        return this;
        }
    T* operator[] (int x){ return dat[x]; }
    void free()
        {
        if (dat!=NULL)
            {
            if (dat[0]!=NULL)
                {
                if (dat[0][0]!=NULL) delete[] dat[0][0];
                delete[] dat[0];
                }
            delete[] dat;
            }
        nx=0; ny=0; nz=0; dat=NULL;
        }
    int  allocate(int _nx,int _ny,int _nz)
        {
        int x,y;
        free();
        dat=new T**[_nx];
        dat[0]=new T*[_nx*_ny];
        dat[0][0]=new T[_nx*_ny*_nz];
        nx=_nx; ny=_ny; nz=_nz;
        for (x=1;x<nx;x++) dat[x]=dat[x-1]+ny;
        for (x=0;x<nx;x++)
         for (y=0;y<ny;y++)
          dat[x][y]=dat[0][0]+(x*ny*nz)+(y*nz);
        }
    };
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------

gl_simple.h just simple OpenGL

parametric_grid.h this is the main source you need to digest

//---------------------------------------------------------------------------
#ifndef _parametric_grid_h
#define _parametric_grid_h
//---------------------------------------------------------------------------
#include "arrays.h"
//---------------------------------------------------------------------------
struct point
    {
    float x,y,z;
    point()     {}
    point(point& a) { *this=a; }
    ~point()    {}
    point* operator = (const point *a) { *this=*a; return this; }
    //point* operator = (const point &a) { ...copy... return this; }
    bool point::operator==(const point &a)
        {
        float d;
        d =(x-a.x)*(x-a.x);
        d+=(y-a.y)*(y-a.y);
        d+=(z-a.z)*(z-a.z);
        return (d<=0.001);
        }
    bool point::operator!=(const point &a)
        {
        float d;
        d =(x-a.x)*(x-a.x);
        d+=(y-a.y)*(y-a.y);
        d+=(z-a.z)*(z-a.z);
        return (d>0.001);
        }
    };
typedef array1D<float>  param[3];   // parametrers <0,1>
typedef array2D<point>  surface[6]; // surface
typedef array3D<point>  volume;     // volume
//---------------------------------------------------------------------------
void swap_xy(array2D<point> &f)
    {
    int ix,iy;
    array2D<point> f0;
    f0=f;
    f.allocate(f0.ny,f0.nx);
    for (ix=0;ix<f.nx;ix++)
     for (iy=0;iy<f.ny;iy++)
      f.dat[ix][iy]=f0.dat[iy][ix];
    }
//---------------------------------------------------------------------------
void mirror_x(array2D<point> &f)
    {
    int ix0,ix1,iy;
    point p;
    for (ix0=0,ix1=f.nx-1;ix0<ix1;ix0++,ix1--)
     for (iy=0;iy<f.ny;iy++)
      { p=f.dat[ix0][iy]; f.dat[ix0][iy]=f.dat[ix1][iy]; f.dat[ix1][iy]=p; }
    }
//---------------------------------------------------------------------------
void mirror_y(array2D<point> &f)
    {
    int ix,iy0,iy1;
    point p;
    for (ix=0;ix<f.nx;ix++)
     for (iy0=0,iy1=f.ny-1;iy0<iy1;iy0++,iy1--)
      { p=f.dat[ix][iy0]; f.dat[ix][iy0]=f.dat[ix][iy1]; f.dat[ix][iy1]=p; }
    }
//---------------------------------------------------------------------------
void surface2volume(volume &v,surface &s)   // normalize surface and convert it to volume
    {
    point p;
    float *pp=(float*)&p;
    int ix,iy,iz,nx,ny,nz,mx,my,i,j,k,id[6],us[6];

    // find which surface belongs to which face f0..f5 and store it to id[]
    // also mirror/rotate to match local x,y orientations
    id[0]=0; us[0]=0;                   // first face hard coded as 0
    for (i=1;i<6;i++) id[i]=-1;         // all the rest not known yet
    for (i=1;i<6;i++) us[i]=-1;         // unused
    nx=s[0].nx;
    ny=s[0].ny;

    for (iy=0,j=id[0],k=4,i=1;i<6;i++)  // find f4 so it share f0(iy=0) edge
     if (us[i]<0)                       // only in unknonw yet surfaces
        {
        mx=s[i].nx-1;
        my=s[i].ny-1;
        if ((s[j].nx==s[i].nx)&&(s[j].ny==s[i].ny)) // try x,y matched
            {
            for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[   ix][   iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k;                                 swap_xy(s[i]); break; }
            for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[   ix][my-iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k;                 mirror_y(s[i]); swap_xy(s[i]); break; }
            for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][   iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]);                 swap_xy(s[i]); break; }
            for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][my-iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); swap_xy(s[i]); break; }
            }
        if ((s[j].nx==s[i].ny)&&(s[j].ny==s[i].nx)) // try x,y swapped
            {
            for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[   iy][   ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k;                                 break; }
            for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[   iy][my-ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k;                 mirror_y(s[i]); break; }
            for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][   ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]);                 break; }
            for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][my-ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); break; }
            }
        } if (id[k]<0) return;          // fail
    nz=s[id[k]].nx;                     // 3th resolution
    v.allocate(nx,ny,nz);               // allocate volume container

    for (iy=ny-1,j=id[0],k=5,i=1;i<6;i++)// find f5 so it share f0(iy=ny-1) edge
     if (us[i]<0)                       // only in unknonw yet surfaces
        {
        mx=s[i].nx-1;
        my=s[i].ny-1;
        if ((s[j].nx==s[i].nx)&&(s[j].ny==s[i].ny)) // try x,y matched
            {
            for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[   ix][   iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k;                                 swap_xy(s[i]); break; }
            for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[   ix][my-iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k;                 mirror_y(s[i]); swap_xy(s[i]); break; }
            for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][   iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]);                 swap_xy(s[i]); break; }
            for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][my-iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); swap_xy(s[i]); break; }
            }
        if ((s[j].nx==s[i].ny)&&(s[j].ny==s[i].nx)) // try x,y swapped
            {
            for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[   iy][   ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]);                 break; }
            for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[   iy][my-ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); break; }
            for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][   ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k;                                 break; }
            for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][my-ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k;                 mirror_y(s[i]); break; }
            }
        } if (id[k]<0) return;          // fail

    for (ix=0,j=id[0],k=2,i=1;i<6;i++)  // find f2 so it share f0(ix=0) edge
     if (us[i]<0)                       // only in unknonw yet surfaces
        {
        mx=s[i].nx-1;
        my=s[i].ny-1;
        if ((s[j].nx==s[i].nx)&&(s[j].ny==s[i].ny)) // try x,y matched
            {
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[   ix][   iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k;                                 break; }
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[   ix][my-iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k;                 mirror_y(s[i]); break; }
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][   iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]);                 break; }
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][my-iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); break; }
            }
        if ((s[j].nx==s[i].ny)&&(s[j].ny==s[i].nx)) // try x,y swapped
            {
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[   iy][   ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k;                                 swap_xy(s[i]); break; }
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[   iy][my-ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k;                 mirror_y(s[i]); swap_xy(s[i]); break; }
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][   ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]);                 swap_xy(s[i]); break; }
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][my-ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); swap_xy(s[i]); break; }
            }
        } if (id[k]<0) return;          // fail

    for (ix=nx-1,j=id[0],k=3,i=1;i<6;i++)// find f3 so it share f0(ix=nx-1) edge
     if (us[i]<0)                       // only in unknonw yet surfaces
        {
        mx=s[i].nx-1;
        my=s[i].ny-1;
        if ((s[j].nx==s[i].nx)&&(s[j].ny==s[i].ny)) // try x,y matched
            {
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[   ix][   iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]);                 break; }
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[   ix][my-iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); break; }
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][   iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k;                                 break; }
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][my-iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k;                 mirror_y(s[i]); break; }
            }
        if ((s[j].nx==s[i].ny)&&(s[j].ny==s[i].nx)) // try x,y swapped
            {
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[   iy][   ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k;                 mirror_y(s[i]); swap_xy(s[i]); break; }
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[   iy][my-ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k;                                 swap_xy(s[i]); break; }
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][   ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); swap_xy(s[i]); break; }
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][my-ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]);                 swap_xy(s[i]); break; }
            }
        } if (id[k]<0) return;          // fail

    for (ix=nz-1,j=id[2],k=1,i=1;i<6;i++)// find f1 so it share f2(ix=nz-1) edge
     if (us[i]<0)                       // only in unknonw yet surfaces
        {
        mx=s[i].nx-1;
        my=s[i].ny-1;
        if ((s[j].nx==s[i].nx)&&(s[j].ny==s[i].ny)) // try x,y matched
            {
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[   ix][   iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]);                 break; }
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[   ix][my-iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); break; }
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][   iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k;                                 break; }
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][my-iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k;                 mirror_y(s[i]); break; }
            }
        if ((s[j].nx==s[i].ny)&&(s[j].ny==s[i].nx)) // try x,y swapped
            {
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[   iy][   ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k;                 mirror_y(s[i]); swap_xy(s[i]); break; }
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[   iy][my-ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k;                                 swap_xy(s[i]); break; }
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][   ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); swap_xy(s[i]); break; }
            for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][my-ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]);                 swap_xy(s[i]); break; }
            }
        } if (id[k]<0) return;          // fail
    i=0;
    // copy surfaces into volume matrix
    for (ix=0;ix<nx;ix++)
     for (iy=0;iy<ny;iy++)
      for (iz=0;iz<nz;iz++)
        {
        // interior points
        p.x=0.0;
        p.y=0.0;
        p.z=0.0;
        // surface points
        if (iz==   0) p=s[id[0]].dat[ix][iy];
        if (iz==nz-1) p=s[id[1]].dat[ix][iy];
        if (ix==   0) p=s[id[2]].dat[iz][iy];
        if (ix==nx-1) p=s[id[3]].dat[iz][iy];
        if (iy==   0) p=s[id[4]].dat[iz][ix];
        if (iy==ny-1) p=s[id[5]].dat[iz][ix];
        v.dat[ix][iy][iz]=p;
        }
    }
//---------------------------------------------------------------------------
void draw(const surface &s) // render surface
    {
    int ix,iy,i;
    DWORD col[6]=   // colors per face
        {
        0x00202060,
        0x00202030,
        0x00206020,
        0x00203020,
        0x00602020,
        0x00302020
        };
    glBegin(GL_LINES);
    for (i=0;i<6;i++)
        {
        glColor4ubv((BYTE*)(&col[i]));
        for (ix=0;ix<s[i].nx;ix++)
         for (iy=1;iy<s[i].ny;iy++)
            {
            glVertex3fv((float*)&(s[i].dat[ix][iy  ]));
            glVertex3fv((float*)&(s[i].dat[ix][iy-1]));
            }
        for (ix=1;ix<s[i].nx;ix++)
         for (iy=0;iy<s[i].ny;iy++)
            {
            glVertex3fv((float*)&(s[i].dat[ix  ][iy]));
            glVertex3fv((float*)&(s[i].dat[ix-1][iy]));
            }
        }
    glEnd();
    }
//---------------------------------------------------------------------------
void draw(const volume &v)
    {
    int ix,iy,iz,i;
    // all points
    glBegin(GL_POINTS);
    for (ix=0;ix<v.nx;ix++)
     for (iy=0;iy<v.ny;iy++)
      for (iz=0;iz<v.nz;iz++)
       glVertex3fv((float*)&(v.dat[ix][iy][iz]));
    glEnd();
    // surface edges for visual check if points are ordered as should
    glBegin(GL_LINE_STRIP); for (ix=     0,iy=     0,iz=     0;ix<v.nx;ix++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
    glBegin(GL_LINE_STRIP); for (ix=     0,iy=     0,iz=v.nz-1;ix<v.nx;ix++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
    glBegin(GL_LINE_STRIP); for (ix=     0,iy=v.ny-1,iz=     0;ix<v.nx;ix++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
    glBegin(GL_LINE_STRIP); for (ix=     0,iy=v.ny-1,iz=v.nz-1;ix<v.nx;ix++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();

    glBegin(GL_LINE_STRIP); for (ix=     0,iy=     0,iz=     0;iy<v.ny;iy++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
    glBegin(GL_LINE_STRIP); for (ix=     0,iy=     0,iz=v.nz-1;iy<v.ny;iy++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
    glBegin(GL_LINE_STRIP); for (ix=v.nx-1,iy=     0,iz=     0;iy<v.ny;iy++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
    glBegin(GL_LINE_STRIP); for (ix=v.nx-1,iy=     0,iz=v.nz-1;iy<v.ny;iy++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();

    glBegin(GL_LINE_STRIP); for (ix=     0,iy=     0,iz=     0;iz<v.nz;iz++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
    glBegin(GL_LINE_STRIP); for (ix=     0,iy=v.nz-1,iz=     0;iz<v.nz;iz++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
    glBegin(GL_LINE_STRIP); for (ix=v.nx-1,iy=     0,iz=     0;iz<v.nz;iz++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
    glBegin(GL_LINE_STRIP); for (ix=v.nx-1,iy=v.nz-1,iz=     0;iz<v.nz;iz++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
    }
//---------------------------------------------------------------------------
void cube(surface &s,const param &t)    // cube -> surface
    {
    point p;
    float *pp=(float*)&p;
    int iu,iv,nx,ny,nz,i,u,v;
    // count of points per axis
    nx=t[0].nx;
    ny=t[1].nx;
    nz=t[2].nx;
    if ((!nx)||(!ny)||(!nz)) return;
    // create 6 faces grid
    for (i=0;i<6;i++)
        {
        if (i==0){ p.z=t[2].dat[   0]; u=0; v=1; } // z = first coordinate
        if (i==1){ p.z=t[2].dat[nz-1]; u=0; v=1; } // z = last coordinate
        if (i==2){ p.y=t[0].dat[   0]; u=0; v=2; } // y = first coordinate
        if (i==3){ p.y=t[0].dat[ny-1]; u=0; v=2; } // y = last coordinate
        if (i==4){ p.x=t[1].dat[   0]; u=2; v=1; } // x = first coordinate
        if (i==5){ p.x=t[2].dat[nx-1]; u=2; v=1; } // x = last coordinate
        s[i].allocate(t[u].nx,t[v].nx);         // allocate face resolution
        for (iu=0;iu<s[i].nx;iu++)              // process all cells
         for (iv=0;iv<s[i].ny;iv++)
            {
            pp[u]=t[u].dat[iu];
            pp[v]=t[v].dat[iv];
            s[i].dat[iu][iv]=p;
            }
        }        
    }
//---------------------------------------------------------------------------
void cubic(surface &s,const param &t)   // hardcoded cubic curved object -> surface
    {
    point p;
    float *pp=(float*)&p;
    int ix,iy,iz,nx,ny,nz,i,j;

    // 3 hard coded cubic control points
    float cp[3][4][3]=  // [parameter][point][coordinate]
        {
        {{ 0.00,0.00,0.00 },
         { 0.50,0.50,0.00 },
         { 1.00,0.50,0.00 },
         { 1.50,0.00,0.00 }},
        {{ 0.00,0.00,0.00 },
         { 0.10,0.15,0.00 },
         { 0.15,0.30,0.00 },
         { 0.16,0.45,0.00 }},
        {{ 0.00,0.00,0.00 },
         { 0.00,0.10,0.15 },
         { 0.00,0.10,0.30 },
         { 0.00,0.00,0.45 }}
        };
    float ce[3][4][3];  // cubic coeficients
    float tx[4],ty[4],tz[4],d1,d2;
    for (ix=0;ix<3;ix++)
     for (iy=0;iy<3;iy++)
        {
        d1=0.5*(cp[ix][2][iy]-cp[ix][0][iy]);
        d2=0.5*(cp[ix][3][iy]-cp[ix][1][iy]);
        ce[ix][0][iy]=cp[ix][1][iy];
        ce[ix][1][iy]=d1;
        ce[ix][2][iy]=(3.0*(cp[ix][2][iy]-cp[ix][1][iy]))-(2.0*d1)-d2;
        ce[ix][3][iy]=d1+d2+(2.0*(-cp[ix][2][iy]+cp[ix][1][iy]));
        }

    // count of points per axis
    nx=t[0].nx;
    ny=t[1].nx;
    nz=t[2].nx;
    if ((!nx)||(!ny)||(!nz)) return;
    // allocate faces
    s[0].allocate(nx,ny);
    s[1].allocate(nx,ny);
    s[2].allocate(ny,nz);
    s[3].allocate(ny,nz);
    s[4].allocate(nz,nx);
    s[5].allocate(nz,nx);
    // process all points
    for (ix=0;ix<nx;ix++)
     for (iy=0;iy<ny;iy++)
      for (iz=0;iz<nz;iz++)
        {
        // compute t,t^2,t^3 per each parameter
        for (tx[0]=1.0,i=1;i<4;i++) tx[i]=tx[i-1]*t[0].dat[ix];
        for (ty[0]=1.0,i=1;i<4;i++) ty[i]=ty[i-1]*t[1].dat[iy];
        for (tz[0]=1.0,i=1;i<4;i++) tz[i]=tz[i-1]*t[2].dat[iz];
        // compute position as superposition of 3 cubics
        for (i=0;i<3;i++)
            {
            pp[i]=0.0;
            for (j=0;j<4;j++) pp[i]+=ce[0][j][i]*tx[j];
            for (j=0;j<4;j++) pp[i]+=ce[1][j][i]*ty[j];
            for (j=0;j<4;j++) pp[i]+=ce[2][j][i]*tz[j];
            }
        // copy posiiton to corresponding surfaces
        if (iz==   0) s[0].dat[ix][iy]=p;
        if (iz==nz-1) s[1].dat[ix][iy]=p;
        if (ix==   0) s[2].dat[iy][iz]=p;
        if (ix==nx-1) s[3].dat[iy][iz]=p;
        if (iy==   0) s[4].dat[iz][ix]=p;
        if (iy==ny-1) s[5].dat[iz][ix]=p;
        }
    }
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------

usage: simple single form VCL app without any components

//---------------------------------------------------------------------------
#include <vcl.h>            // VCL stuff (ignore)
#pragma hdrstop             // VCL stuff (ignore)
#include "Unit1.h"          // VCL stuff (header of this window)
#include "gl_simple.h"      // my GL init (source included)
#include "parametric_grid.h"
//---------------------------------------------------------------------------
#pragma package(smart_init) // VCL stuff (ignore)
#pragma resource "*.dfm"    // VCL stuff (ignore)
TForm1 *Form1;              // VCL stuff (this window)
surface sur;                // 6 2D matrices (surfaces in any order/orientation before normalization)
volume vol;                 // 3D matrix (volume)
//---------------------------------------------------------------------------
void gl_draw()
    {
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    glDisable(GL_CULL_FACE);
    glEnable(GL_DEPTH_TEST);

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    glTranslatef(-0.9,-0.75,-2.0);
    glRotatef(45.0,0.5,0.5,0.0);

    glPointSize(4.0);
    glColor3f(0.1,0.5,0.7); draw(vol);
    glPointSize(1.0);
//  draw(sur);

    glFlush();
    SwapBuffers(hdc);
    }
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner)
    {
    // this is called on window startup
    gl_init(Handle);                // init OpenGL 1.0

    int i;
    param t;                        // source coordinates
/*
    t[0].allocate(5); i=0;
    t[0][i]=0.000; i++;
    t[0][i]=0.125; i++;
    t[0][i]=0.250; i++;
    t[0][i]=0.500; i++;
    t[0][i]=1.000; i++;

    t[1].allocate(5); i=0;
    t[1][i]=0.000; i++;
    t[1][i]=0.250; i++;
    t[1][i]=0.500; i++;
    t[1][i]=0.750; i++;
    t[1][i]=1.000; i++;

    t[2].allocate(6); i=0;
    t[2][i]=0.0; i++;
    t[2][i]=0.2; i++;
    t[2][i]=0.4; i++;
    t[2][i]=0.6; i++;
    t[2][i]=0.8; i++;
    t[2][i]=1.0; i++;
*/
  
    t[0].allocate(6); i=0;
    t[0][i]=0.0; i++;
    t[0][i]=0.2; i++;
    t[0][i]=0.4; i++;
    t[0][i]=0.6; i++;
    t[0][i]=0.8; i++;
    t[0][i]=1.0; i++;
    t[1].allocate(6); i=0;
    t[1][i]=0.0; i++;
    t[1][i]=0.2; i++;
    t[1][i]=0.4; i++;
    t[1][i]=0.6; i++;
    t[1][i]=0.8; i++;
    t[1][i]=1.0; i++;
    t[2].allocate(6); i=0;
    t[2][i]=0.0; i++;
    t[2][i]=0.2; i++;
    t[2][i]=0.4; i++;
    t[2][i]=0.6; i++;
    t[2][i]=0.8; i++;
    t[2][i]=1.0; i++;


//  cube(sur,t);
    cubic(sur,t);
    surface2volume(vol,sur);
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormDestroy(TObject *Sender)
    {
    // this is called before window exits
    gl_exit();                      // exit OpenGL
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormResize(TObject *Sender)
    {
    // this is called on each window resize (and also after startup)
    gl_resize(ClientWidth,ClientHeight);
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormPaint(TObject *Sender)
    {
    // this is called whnewer app needs repaint
    gl_draw();
    }
//---------------------------------------------------------------------------

This is preview for code above:

I debug by trying all combinations of mirror/swap for each surfac (first does not count as its used as is) and rendering the shared edges if they overlap exactly and tweaked the reordering of found faces until they did so the code should be working for any input (that connects)