将一个三维数组分割成较小的“立方体”

Question

我目前正在做这个工作:我生成一个Paraview .vtm文件，其中包含几个.vtr文件。假设我处理的维度为8，每个.vtr文件都包含值和坐标，如下所示：

<PointData Scalars="U">
  <DataArray type="Float32" Name="U" format="ascii">
    <!-- 8*8*8 values -->
  </DataArray>
</PointData>
<Coordinates>
  <DataArray type="Float32" Name="x" format="ascii">
    <!-- 8 x values -->
  </DataArray>
  <DataArray type="Float32" Name="y" format="ascii">
    <!-- 8 y values -->
  </DataArray>
  <DataArray type="Float32" Name="z" format="ascii">
    <!-- 8 z values -->
  </DataArray>
</Coordinates>

我使用一个四维数组来存储我的值：float ****tab，with tab[s][x][y][z]，其中：

• s是当前的拆分步骤。每次我开始处理下一个.vtr文件时，它都会增加。
• x, y, z值。

现在给我带来麻烦的是:我必须放置这些点的坐标可以是任何东西。它可以是常量(遵循一个步骤，比如0、0.1、0.2等等)，也可以不是。

我将坐标存储在三个数组中：x[], y[], z[]。我的目标是将一组值切割成更小的立方体。假设我将我的值分成8个文件(2^3个文件)，我必须为8个小立方体检索正确的坐标。我也想不出办法来。

我很确定我的数据结构选择很糟糕，有人能帮我一下吗？

编辑:

下面是生成我的四星数组的函数：

float**** fill_array_random4d(int split, int size)
{
  float**** ret;
  ret = malloc(sizeof(float***) * split);
  for (int i = 0; i < split; i++)
  {
    ret[i] = malloc(sizeof (float**) * size);
    for (int j = 0; j < size; j++)
    {
      ret[i][j] = malloc(sizeof (float*) * size);
      for (int k = 0; k < size; k++)
      {
        ret[i][j][k] = malloc(sizeof (float) * size);
        for (int l = 0; l < size; l++)
          ret[i][j][k][l] = rand() % 100;
      }
    }
  }

  return ret;
}

这是个很基本的东西。现在我用的是随机值。下面是我如何创建和填充我的x, y, z数组：

float *x, *y, *z;
  x = malloc(sizeof (float) * size);
  y = malloc(sizeof (float) * size);
  z = malloc(sizeof (float) * size);
  for (int i = 0; i < size * split; i++)
    x[i] = step * i;
  for (int i = 0; i < size * split; i++)
    y[i] = step * i;
  for (int i = 0; i < size * split; i++)
    z[i] = step * i;

它仍然是非常基本的，最后是输出文件中的坐标的函数，遵循vtk遗留格式：

void print_Coordinates(FILE *file, float *x, float *y, float *z, int size, int split)
{
  fprintf(file, "      <Coordinates>\n");
  for (int i = 0; i < 3; i++)
  {
    const char *text1 = "        <DataArray type=\"Float32\" Name=\"";
    const char *text2 = "\" format=\"ascii\">\n";
    fprintf(file, "%s%c%s", text1, 'x' + i, text2);
    for (int j = 0; j < size; j++)
    {
      if (i == 0)
        fprintf(file, "          %f\n", x[j]);
      else if (i == 1)
        fprintf(file, "          %f\n", y[j]);
      else
        fprintf(file, "          %f\n", z[j]);
    }
    fprintf(file, "        </DataArray>\n");
  }
  fprintf(file, "      </Coordinates>\n");
}

所以，是的，它根本做不到我想做的事。下面是结果的截图：

​

所有的立方体都在一起。在我前面使用的代码中，我有几个多维数据集(每个文件一个)，但是它们是对角对齐的(这也不是很好)。

Answer 1

正如您所承认的，您的数据结构存在一些问题：

• 一维s似乎不一致:数据结构是否应该包括原始的和较小的立方体？这不容易做到，因为较小的立方体有其他的尺寸。
• 您有许多独立的数据：(随机)数据、坐标和数组尺寸。为了表示多维数据集，您需要跟踪所有这些。我建议创建一个结构来将相关数据放在一起。
• 用三元指针表示三维数组的方法本身并没有什么问题，但是设计导致了许多零碎的分配。具有恒定维数的多维数组可能更好地表示为一个“平面”内存块。

我建议两种结构：

typedef struct Cube Cube;
typedef struct Axis Axis;

struct Axis {
    int n;              /* number of values */
    float *data;        /* graduation values */
};

struct Cube {
    Axis *x, *y, *z;    /* Axes of the cube */
    float *data;        /* x-major data */
};

“轴”沿轴之一存储值。多维数据集本身并不担心与轴相关的代码，而只是将其委托给它的三个成员轴。“多维数据集”是您的数据对象。(在下面的实现中，数据表示为x- meaning，这意味着x循环是最外层的，z循环是最内部的。你可以通过交换回路来达到这个目的。

如果您有一个填充的多维数据集对象，则可以通过创建一个较小维度的多维数据集并从轴和多维数据集数据复制相关数据范围来提取子多维数据集。如果你想要覆盖整个立方体，你可以在你走的时候提取并写这些立方体，或者把它们存储在一个立方体数组中，例如Cube *small[8]，用于分割成两个方向。(这就像您原来的s索引，只是每个多维数据集可能有它自己的维度。)

下面是使用(附加简单的)测试main实现这种行为的方法：

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

typedef struct Cube Cube;
typedef struct Axis Axis;

struct Axis {
    int n;              /* number of values */
    float *data;        /* graduation values */
};

struct Cube {
    Axis *x, *y, *z;    /* Axes of the cube */
    float *data;        /* x-major data */
};

/*
 *      Create a new axis with a constant step.
 */
Axis *axis_new(int n, float start, float step)
{
    Axis *axis = malloc(sizeof(*axis));
    float *p;

    axis->n = n;
    axis->data = malloc(n * sizeof(*axis->data));

    p = axis->data;
    while (n--) {
        *p = start;
        start += step;
        p++;
    }

    return axis;
}

/*
 *      Destroy and clean up axis
 */
void axis_delete(Axis *axis)
{
    if (axis) {
        free(axis->data);
        free(axis);
    }
}

/*
 *      Write axis in XML format to given file
 */
void axis_write(const Axis *axis, FILE *f, const char *name)
{
    float *p = axis->data;
    int n = axis->n;

    fprintf(f, "  <DataArray type=\"Float32\" "
        "Name=\"%s\" format=\"ascii\">\n", name);

    fprintf(f, "    ");
    while (n--) {
        fprintf(f, " %g", *p++);
    }
    fprintf(f, "\n");
    fprintf(f, "  </DataArray>\n");
}

/*
 *      Create a new axis that is a sub-axis of orig.
 */
Axis *axis_slice(const Axis *orig, int start, int len)
{
    Axis *axis = axis_new(len, 0, 0);
    memcpy(axis->data, orig->data + start, len * sizeof(*axis->data));

    return axis;
}



/*
 *      Create a cube of zero values for the given axes
 */
Cube *cube_new(Axis *x, Axis *y, Axis *z)
{
    Cube *cube = malloc(sizeof(*cube));
    int dim = x->n * y->n * z->n;

    cube->x = x;
    cube->y = y;
    cube->z = z;

    cube->data = malloc(dim * sizeof(*cube->data));

    return cube;
}

/*
 *      Destroy and clean up cube
 */
void cube_delete(Cube *cube)
{
    if (cube) {
        axis_delete(cube->x);
        axis_delete(cube->y);
        axis_delete(cube->z);

        free(cube->data);
        free(cube);
    }
}

float *cube_at(const Cube *cube, int x, int y, int z)
{
    int pos = (x * cube->y->n + y) * cube->z->n + z;
    return cube->data + pos;
}

/*
 *      Populate all x, y, z values according to the function func.
 */
void cube_populate(Cube *cube, float (*func)(float x, float y, float z))
{
    int i, j, k;
    float *p = cube->data;

    for (i = 0; i < cube->x->n; i++) {
        float x = cube->x->data[i];

        for (j = 0; j < cube->y->n; j++) {
            float y = cube->y->data[j];

            for (k = 0; k < cube->z->n; k++) {
                float z = cube->z->data[k];

                *p++ = func(x, y, z);
            }
        }
    }    
}

/*
 *      Write cube to given file.
 */
void cube_write(const Cube *cube, FILE *f)
{
    float *p = cube->data;
    int n = cube->x->n * cube->y->n * cube->z->n;

    fprintf(f, "<PointData Scalars=\"U\">\n");
    fprintf(f, "  <DataArray type=\"Float32\" Name=\"U\" format=\"ascii\">\n");

    while (n--) {
        fprintf(f, " %g", *p++);
    }
    fprintf(f, "\n");

    fprintf(f, "  </DataArray>\n");
    fprintf(f, "</PointData>\n");

    fprintf(f, "<Coordinates>\n");
    axis_write(cube->x, f, "x");
    axis_write(cube->y, f, "y");
    axis_write(cube->z, f, "z");
    fprintf(f, "</Coordinates>\n");
}

/*
 *      Create a new cube that is a sub-cube of orig.
 */
Cube *cube_slice(const Cube *orig,
    int x, int dx, int y, int dy, int z, int dz)
{
    Cube *cube;
    float *p;
    int i, j, k;

    if (x + dx > orig->x->n) return NULL;
    if (y + dy > orig->y->n) return NULL;
    if (z + dz > orig->z->n) return NULL;

    cube = cube_new(
        axis_slice(orig->x, x, dx),
        axis_slice(orig->y, y, dy),
        axis_slice(orig->z, z, dz));

    p = cube->data;

    for (i = 0; i < dx; i++) {
        for (j = 0; j < dy; j++) {
            for (k = 0; k < dz; k++) {
                *p++ = *cube_at(orig, x + i, y + j, z + k);
            }
        }
    }    

    return cube;
}

/*
 *      Example appliaction
 */
float dist2(float x, float y, float z)
{
    return x*x + y*y + z*z;
}

int main()
{
    Cube *cube = cube_new(
        axis_new(4, 0, 0.1),
        axis_new(4, 0, 0.1),
        axis_new(4, 0, 0.1));
    int i, j, k;

    cube_populate(cube, dist2);

    for (i = 0; i < 2; i++) {
        for (j = 0; j < 2; j++) {
            for (k = 0; k < 2; k++) {
                Cube *sub = cube_slice(cube, 2*i, 2, 2*j, 2, 2*k, 2);

                cube_write(sub, stdout);
                printf("--\n");
                cube_delete(sub);
            }
        }
    }

    cube_delete(cube);

    return 0;
}