c - 将三维数组分割成更小的 "cubes"

Question

我目前正在研究这个:我生成一个包含多个 .vtr 文件的 Paraview .vtm 文件。每个 .vtr 文件都包含值和坐标，如下所示，假设我正在处理 8 的维度:

<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，带有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 */
};

“轴”存储沿其中一个轴的值​​。立方体本身并不担心与轴相关的代码，只是将其委托(delegate)给它的三个成员轴。 “立方体”是您的数据对象。 (在下面的实现中，数据表示为 x-major，这意味着 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;
}