grid.c

#include <stdio.h>
#include <string.h>
#include <errno.h>
#include "grid.h"

/*
 * This file is in public domain.
*/

static unsigned long  field_evaluations = 0UL;

static double field(void *unused, double x, double y, double z)
{
    field_evaluations++;
    return x*x*x + y*y*y + z*z*z + x + y - 0.125 * (x*y + x*z + y*z + x*y*z);
}

static void check_edge(FILE *out, const double s,
                       const double x0, const double y0, const double z0, const double s0,
                       const double x1, const double y1, const double z1, const double s1)
{
    if (s0 < s && s < s1) {
        const double p = (s - s0) / (s1 - s0);
        fprintf(out, "%.5f %.5f %.5f\n", x0 + p*(x1-x0), y0 + p*(y1-y0), z0 + p*(z1-z0));
    } else
    if (s0 > s && s > s1) {
        const double p = (s - s1) / (s0 - s1);
        fprintf(out, "%.5f %.5f %.5f\n", x1 + p*(x0-x1), y1 + p*(y0-y1), z1 + p*(z0-z1));
    } else
    if (s0 == s && s1 != s) {
        fprintf(out, "%.5f %.5f %.5f\n", x0, y0, z0);
    } else
    if (s0 != s && s1 == s) {
        fprintf(out, "%.5f %.5f %.5f\n", x1, y1, z1);
    }
}

static unsigned long  cells_checked = 0UL;

static int check_cell(void *out, struct grid *g, int x, int y, int z, double c, double sample[2][2][2])
{
    const double x0 = g->xorigin + (double)x * g->xunit,
                 y0 = g->yorigin + (double)y * g->yunit,
                 z0 = g->zorigin + (double)z * g->zunit;
    const double x1 = x0 + g->xunit,
                 y1 = y0 + g->yunit,
                 z1 = z0 + g->zunit;

    cells_checked++;

#ifdef CELL_EDGE_INTERSECTIONS

    /* These are the 12 possible intersections at the cell edges,
     * in case you intend to use them for e.g. polygons.
    */

    check_edge(out, c, x0, y0, z0, sample[0][0][0], x1, y0, z0, sample[1][0][0]);
    check_edge(out, c, x0, y0, z0, sample[0][0][0], x0, y1, z0, sample[0][1][0]);
    check_edge(out, c, x0, y0, z0, sample[0][0][0], x0, y0, z1, sample[0][0][1]);

    check_edge(out, c, x1, y0, z0, sample[1][0][0], x1, y1, z0, sample[1][1][0]);
    check_edge(out, c, x1, y0, z0, sample[1][0][0], x1, y0, z1, sample[1][0][1]);

    check_edge(out, c, x0, y1, z0, sample[0][1][0], x1, y1, z0, sample[1][1][0]);
    check_edge(out, c, x0, y1, z0, sample[0][1][0], x0, y1, z1, sample[0][1][1]);

    check_edge(out, c, x0, y0, z1, sample[0][0][1], x1, y0, z1, sample[1][0][1]);
    check_edge(out, c, x0, y0, z1, sample[0][0][1], x0, y1, z1, sample[0][1][1]);

    check_edge(out, c, x1, y1, z0, sample[1][1][0], x1, y1, z1, sample[1][1][1]);
    check_edge(out, c, x1, y0, z1, sample[1][0][1], x1, y1, z1, sample[1][1][1]);
    check_edge(out, c, x0, y1, z1, sample[0][1][1], x1, y1, z1, sample[1][1][1]);

#else

    /* These are not (all) valid edges, but if you are just visualizing the
     * approximate isosurface as a point cloud, these give a slightly better
     * distributed point set. */

    check_edge(out, c, x0, y0, z0, sample[0][0][0], x1, y0, z0, sample[1][0][0]);
    check_edge(out, c, x0, y0, z0, sample[0][0][0], x0, y1, z0, sample[0][1][0]);
    check_edge(out, c, x0, y0, z0, sample[0][0][0], x0, y0, z1, sample[0][0][1]);
    check_edge(out, c, x0, y0, z0, sample[0][0][0], x0, y1, z1, sample[0][1][1]);
    check_edge(out, c, x0, y0, z0, sample[0][0][0], x1, y0, z1, sample[1][0][1]);
    check_edge(out, c, x0, y0, z0, sample[0][0][0], x1, y1, z0, sample[1][1][0]);
    check_edge(out, c, x0, y0, z0, sample[0][0][0], x1, y1, z1, sample[1][1][1]);

#endif

    return 0;
}



int main(int argc, char *argv[])
{
    struct grid *g;
    double min, max, level;
    int n, result;
    char dummy;

    if (argc != 5 || !strcmp(argv[1], "-h") || !strcmp(argv[1], "--help")) {
        fprintf(stderr, "\n");
        fprintf(stderr, "Usage: %s [ -h | --help ]\n", argv[0]);
        fprintf(stderr, "       %s N MIN MAX LEVEL\n", argv[0]);
        fprintf(stderr, "\n");
        return EXIT_FAILURE;
    }

    if (sscanf(argv[1], " %d %c", &n, &dummy) != 1 || n < 1) {
        fprintf(stderr, "%s: Invalid N.\n", argv[1]);
        return EXIT_FAILURE;
    }
    if (sscanf(argv[2], " %lf %c", &min, &dummy) != 1) {
        fprintf(stderr, "%s: Invalid minimum.\n", argv[2]);
        return EXIT_FAILURE;
    }
    if (sscanf(argv[3], " %lf %c", &max, &dummy) != 1 || max <= min) {
        fprintf(stderr, "%s: Invalid maximum.\n", argv[3]);
        return EXIT_FAILURE;
    }
    if (sscanf(argv[4], " %lf %c", &level, &dummy) != 1) {
        fprintf(stderr, "%s: Invalid level.\n", argv[4]);
        return EXIT_FAILURE;
    }

    g = grid_new(n+1, n+1, n+1, min, min, min, max, max, max, NULL, field);
    if (g == NULL) {
        fprintf(stderr, "Cannot construct grid: %s.\n", strerror(errno));
        return EXIT_FAILURE;
    }
    if (level < g->minimum_sample || level > g->maximum_sample) {
        fprintf(stderr, "%s: Level is outside the field (%g to %g).\n", argv[4], g->minimum_sample, g->maximum_sample);
        return EXIT_FAILURE;
    }

    fprintf(stderr, "Grid initialized using %lu samples,\n", field_evaluations);
    fprintf(stderr, "field (%.3f, %.3f, %.3f) - (%.3f, %.3f, %.3f).\n",
                    g->xorigin, g->yorigin, g->zorigin,
                    g->xorigin + (double)(g->xsize - 1) * g->xunit,
                    g->yorigin + (double)(g->ysize - 1) * g->yunit,
                    g->zorigin + (double)(g->zsize - 1) * g->zunit);

    result = grid_isosurface(g, level, stdout, check_cell);
    if (result != 0) {
        fprintf(stderr, "Failed to locate isosurface [%d]: %s.\n", result, strerror(errno));
        return EXIT_FAILURE;
    }

    fprintf(stderr, "Isosurface passed through %lu cells (%.5f%%).\n", cells_checked, 100.0 * (double)cells_checked / (double)(g->xsize - 1) / (double)(g->ysize - 1) / (double)(g->zsize - 1));
    fprintf(stderr, "Initialization and isosurface required %lu samples (%.5f%%) total.\n", field_evaluations, 100.0 * (double)field_evaluations / (double)g->size);

    g = grid_free(g);

    return EXIT_SUCCESS;
}