API tools faq
paste
Advertisement
Guest User

grid.h

a guest
Dec 18th, 2014
342
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 26.86 KB | None | 0 0
raw download clone embed print report
  1. #ifndef   GRID_H
  2. #define   GRID_H
  3. #include <stdlib.h>
  4. #include <string.h>
  5. #include <errno.h>
  6.  
  7. /*
  8.  * This file is in public domain.
  9. */
  10.  
  11. typedef struct {
  12.     double          minimum;
  13.     double          maximum;
  14. } span;
  15.  
  16. struct grid {
  17.     /* Grid size */
  18.     int             xsize;
  19.     int             ysize;
  20.     int             zsize;
  21.  
  22.     /* Total number of grid points */
  23.     size_t          size;   /* xsize * ysize * zsize */
  24.  
  25.     /* Strides, allowing any major */
  26.     size_t          xstride;  /* 1 */
  27.     size_t          ystride;  /* xsize */
  28.     size_t          zstride;  /* xsize * ysize */
  29.  
  30.     /* Cell containing minimum/maximum samples */
  31.     size_t          minimum_cell;
  32.     size_t          maximum_cell;
  33.     double          minimum_sample;
  34.     double          maximum_sample;
  35.  
  36.     /* Location of the (0,0,0) grid point */
  37.     double          xorigin;
  38.     double          yorigin;
  39.     double          zorigin;
  40.  
  41.     /* Grid unit distances along each axis */
  42.     double          xunit;
  43.     double          yunit;
  44.     double          zunit;
  45.  
  46.     /* Scalar field sampling function */
  47.     void           *field_data;
  48.     double        (*field)(void *field_data, double x, double y, double z);
  49.  
  50.     /* Location stack */
  51.     size_t          stack_size;
  52.     size_t          stack_used;
  53.     size_t         *stack;
  54.  
  55.     /* Grid cells */
  56.     unsigned char  *cell;
  57.  
  58.     /* Sample cache */
  59.     double          sample[];
  60. };
  61.  
  62. #define  CELL_UNKNOWN  0U
  63. #define  CELL_SAMPLED  1U
  64. #define  CELL_PUSHED   2U
  65. #define  CELL_WALKED   4U
  66.  
  67. #define  GRID_X(grid, index) (((index) / (grid)->xstride) % (grid)->xsize)
  68. #define  GRID_Y(grid, index) (((index) / (grid)->ystride) % (grid)->ysize)
  69. #define  GRID_Z(grid, index) (((index) / (grid)->zstride) % (grid)->zsize)
  70. #define  GRID_I(grid, x, y, z) ((size_t)(x) * (grid)->xstride + (size_t)(y) * (grid)->ystride + (size_t)(z) * (grid)->zstride)
  71.  
  72.  
  73. static inline void grid_clear(struct grid *const g)
  74. {
  75.     unsigned char       *p = g->cell;
  76.     unsigned char *const q = g->cell + g->size;
  77.  
  78.     /* Retain only CELL_SAMPLED flag. */
  79.     while (p < q)
  80.         *(p++) &= CELL_SAMPLED;
  81.  
  82.     /* Clear the stack, too. */
  83.     g->stack_used = 0;
  84. }
  85.  
  86.  
  87. static inline size_t grid_pop(struct grid *const g)
  88. {
  89.     /* Skip already walked indices. */
  90.     while (g->stack_used > 0 && (g->cell[g->stack[g->stack_used - 1]] & CELL_WALKED) != 0)
  91.         g->stack_used--;
  92.  
  93.     if (g->stack_used > 0) {
  94.         const size_t i = g->stack[--g->stack_used];
  95.         /* Mark popped index walked, and return it. */
  96.         g->cell[i] |= CELL_WALKED;
  97.         return i;
  98.     } else
  99.         return g->size; /* Stack empty. */
  100. }
  101.  
  102.  
  103. static inline int grid_push(struct grid *const g, const size_t i)
  104. {
  105.     /* Do not push already pushed or walked cells. */
  106.     if (g->cell[i] & (CELL_PUSHED | CELL_WALKED))
  107.         return 0;
  108.  
  109.     /* Grow stack if necessary. */
  110.     if (g->stack_used >= g->stack_size) {
  111.         const size_t  new_size = (g->stack_used | 1023) + 1025;
  112.         size_t       *new_stack;
  113.  
  114.         new_stack = realloc(g->stack, new_size * sizeof g->stack[0]);
  115.         if (new_stack == NULL)
  116.             return errno = ENOMEM;
  117.  
  118.         g->stack = new_stack;
  119.         g->stack_size = new_size;
  120.     }
  121.  
  122.     /* Mark cell pushed, */
  123.     g->cell[i] |= CELL_PUSHED;
  124.  
  125.     /* and append it to the stack. */
  126.     g->stack[g->stack_used++] = i;
  127.  
  128.     return 0;
  129. }
  130.  
  131.  
  132. static inline double grid_sample(struct grid *const g, const size_t x, const size_t y, const size_t z)
  133. {
  134.     const size_t i = x * g->xstride + y * g->ystride + z * g->zstride;
  135.  
  136.     if (!(g->cell[i] & CELL_SAMPLED)) {
  137.         g->sample[i] = g->field(g->field_data, g->xorigin + (double)x * g->xunit,
  138.                                                g->yorigin + (double)y * g->yunit,
  139.                                                g->zorigin + (double)z * g->zunit);
  140.         g->cell[i] |= CELL_SAMPLED;
  141.     }
  142.  
  143.     return g->sample[i];
  144. }
  145.  
  146. static span grid_cell(struct grid *const g, const size_t x, const size_t y, const size_t z, double corner[2][2][2])
  147. {
  148.     const size_t xstride = g->xstride,
  149.                  ystride = g->ystride,
  150.                  zstride = g->zstride;
  151.     const double field_x = g->xorigin + (double)x * g->xunit,
  152.                  field_y = g->yorigin + (double)y * g->yunit,
  153.                  field_z = g->zorigin + (double)z * g->zunit,
  154.                  unit_x = g->xunit,
  155.                  unit_y = g->yunit,
  156.                  unit_z = g->zunit;
  157.     unsigned char *const cell = g->cell + x * xstride + y * ystride + z * zstride;
  158.     double *const sample = g->sample + x * xstride + y * ystride + z * zstride;
  159.     span result;
  160.  
  161.     if (!(cell[0] & CELL_SAMPLED)) {
  162.         sample[0] = g->field(g->field_data, field_x, field_y, field_z);
  163.         cell[0] |= CELL_SAMPLED;
  164.     }
  165.     corner[0][0][0] = sample[0];
  166.  
  167.     {   const size_t i = xstride;
  168.         if (!(cell[i] & CELL_SAMPLED)) {
  169.             sample[i] = g->field(g->field_data, field_x + unit_x, field_y, field_z);
  170.             cell[i] |= CELL_SAMPLED;
  171.         }
  172.         corner[1][0][0] = sample[i];
  173.     }
  174.  
  175.     {   const size_t i = ystride;
  176.         if (!(cell[i] & CELL_SAMPLED)) {
  177.             sample[i] = g->field(g->field_data, field_x, field_y + unit_y, field_z);
  178.             cell[i] |= CELL_SAMPLED;
  179.         }
  180.         corner[0][1][0] = sample[i];
  181.     }
  182.  
  183.     {   const size_t i = xstride + ystride;
  184.         if (!(cell[i] & CELL_SAMPLED)) {
  185.             sample[i] = g->field(g->field_data, field_x + unit_x, field_y + unit_y, field_z);
  186.             cell[i] |= CELL_SAMPLED;
  187.         }
  188.         corner[1][1][0] = sample[i];
  189.     }
  190.  
  191.     {   const size_t i = zstride;
  192.         if (!(cell[i] & CELL_SAMPLED)) {
  193.             sample[i] = g->field(g->field_data, field_x, field_y, field_z + unit_z);
  194.             cell[i] |= CELL_SAMPLED;
  195.         }
  196.         corner[0][0][1] = sample[i];
  197.     }
  198.  
  199.     {   const size_t i = xstride + zstride;
  200.         if (!(cell[i] & CELL_SAMPLED)) {
  201.             sample[i] = g->field(g->field_data, field_x + unit_x, field_y, field_z + unit_z);
  202.             cell[i] |= CELL_SAMPLED;
  203.         }
  204.         corner[1][0][1] = sample[i];
  205.     }
  206.  
  207.     {   const size_t i = ystride + zstride;
  208.         if (!(cell[i] & CELL_SAMPLED)) {
  209.             sample[i] = g->field(g->field_data, field_x, field_y + unit_y, field_z + unit_z);
  210.             cell[i] |= CELL_SAMPLED;
  211.         }
  212.         corner[0][1][1] = sample[i];
  213.     }
  214.  
  215.     {   const size_t i = xstride + ystride + zstride;
  216.         if (!(cell[i] & CELL_SAMPLED)) {
  217.             sample[i] = g->field(g->field_data, field_x + unit_x, field_y + unit_y, field_z + unit_z);
  218.             cell[i] |= CELL_SAMPLED;
  219.         }
  220.         corner[1][1][1] = sample[i];
  221.     }
  222.  
  223.     result.minimum = corner[0][0][0];
  224.     if (result.minimum > corner[0][0][1])
  225.         result.minimum = corner[0][0][1];
  226.     if (result.minimum > corner[0][1][0])
  227.         result.minimum = corner[0][1][0];
  228.     if (result.minimum > corner[0][1][1])
  229.         result.minimum = corner[0][1][1];
  230.     if (result.minimum > corner[1][0][0])
  231.         result.minimum = corner[1][0][0];
  232.     if (result.minimum > corner[1][0][1])
  233.         result.minimum = corner[1][0][1];
  234.     if (result.minimum > corner[1][1][0])
  235.         result.minimum = corner[1][1][0];
  236.     if (result.minimum > corner[1][1][1])
  237.         result.minimum = corner[1][1][1];
  238.  
  239.     result.maximum = corner[0][0][0];
  240.     if (result.maximum < corner[0][0][1])
  241.         result.maximum = corner[0][0][1];
  242.     if (result.maximum < corner[0][1][0])
  243.         result.maximum = corner[0][1][0];
  244.     if (result.maximum < corner[0][1][1])
  245.         result.maximum = corner[0][1][1];
  246.     if (result.maximum < corner[1][0][0])
  247.         result.maximum = corner[1][0][0];
  248.     if (result.maximum < corner[1][0][1])
  249.         result.maximum = corner[1][0][1];
  250.     if (result.maximum < corner[1][1][0])
  251.         result.maximum = corner[1][1][0];
  252.     if (result.maximum < corner[1][1][1])
  253.         result.maximum = corner[1][1][1];
  254.  
  255.     return result;
  256. }
  257.  
  258.  
  259. static inline span grid_cell_span(struct grid *const g, const size_t x, const size_t y, const size_t z)
  260. {
  261.     const size_t xstride = g->xstride,
  262.                  ystride = g->ystride,
  263.                  zstride = g->zstride;
  264.     const double field_x = g->xorigin + (double)x * g->xunit,
  265.                  field_y = g->yorigin + (double)y * g->yunit,
  266.                  field_z = g->zorigin + (double)z * g->zunit,
  267.                  unit_x = g->xunit,
  268.                  unit_y = g->yunit,
  269.                  unit_z = g->zunit;
  270.     unsigned char *const cell = g->cell + x * xstride + y * ystride + z * zstride;
  271.     double *const sample = g->sample + x * xstride + y * ystride + z * zstride;
  272.     span result;
  273.  
  274.     if (!(cell[0] & CELL_SAMPLED)) {
  275.         sample[0] = g->field(g->field_data, field_x, field_y, field_z);
  276.         cell[0] |= CELL_SAMPLED;
  277.     }
  278.     result.minimum = sample[0];
  279.     result.maximum = sample[0];
  280.  
  281.     {   const size_t i = xstride;
  282.         if (!(cell[i] & CELL_SAMPLED)) {
  283.             sample[i] = g->field(g->field_data, field_x + unit_x, field_y, field_z);
  284.             cell[i] |= CELL_SAMPLED;
  285.         }
  286.         {   const double s = sample[i];
  287.             if (result.minimum > s)
  288.                 result.minimum = s;
  289.             if (result.maximum < s)
  290.                 result.maximum = s;
  291.         }
  292.     }
  293.  
  294.     {   const size_t i = ystride;
  295.         if (!(cell[i] & CELL_SAMPLED)) {
  296.             sample[i] = g->field(g->field_data, field_x, field_y + unit_y, field_z);
  297.             cell[i] |= CELL_SAMPLED;
  298.         }
  299.         {   const double s = sample[i];
  300.             if (result.minimum > s)
  301.                 result.minimum = s;
  302.             if (result.maximum < s)
  303.                 result.maximum = s;
  304.         }
  305.     }
  306.  
  307.     {   const size_t i = xstride + ystride;
  308.         if (!(cell[i] & CELL_SAMPLED)) {
  309.             sample[i] = g->field(g->field_data, field_x + unit_x, field_y + unit_y, field_z);
  310.             cell[i] |= CELL_SAMPLED;
  311.         }
  312.         {   const double s = sample[i];
  313.             if (result.minimum > s)
  314.                 result.minimum = s;
  315.             if (result.maximum < s)
  316.                 result.maximum = s;
  317.         }
  318.     }
  319.  
  320.     {   const size_t i = zstride;
  321.         if (!(cell[i] & CELL_SAMPLED)) {
  322.             sample[i] = g->field(g->field_data, field_x, field_y, field_z + unit_z);
  323.             cell[i] |= CELL_SAMPLED;
  324.         }
  325.         {   const double s = sample[i];
  326.             if (result.minimum > s)
  327.                 result.minimum = s;
  328.             if (result.maximum < s)
  329.                 result.maximum = s;
  330.         }
  331.     }
  332.  
  333.     {   const size_t i = xstride + zstride;
  334.         if (!(cell[i] & CELL_SAMPLED)) {
  335.             sample[i] = g->field(g->field_data, field_x + unit_x, field_y, field_z + unit_z);
  336.             cell[i] |= CELL_SAMPLED;
  337.         }
  338.         {   const double s = sample[i];
  339.             if (result.minimum > s)
  340.                 result.minimum = s;
  341.             if (result.maximum < s)
  342.                 result.maximum = s;
  343.         }
  344.     }
  345.  
  346.     {   const size_t i = ystride + zstride;
  347.         if (!(cell[i] & CELL_SAMPLED)) {
  348.             sample[i] = g->field(g->field_data, field_x, field_y + unit_y, field_z + unit_z);
  349.             cell[i] |= CELL_SAMPLED;
  350.         }
  351.         {   const double s = sample[i];
  352.             if (result.minimum > s)
  353.                 result.minimum = s;
  354.             if (result.maximum < s)
  355.                 result.maximum = s;
  356.         }
  357.     }
  358.  
  359.     {   const size_t i = xstride + ystride + zstride;
  360.         if (!(cell[i] & CELL_SAMPLED)) {
  361.             sample[i] = g->field(g->field_data, field_x + unit_x, field_y + unit_y, field_z + unit_z);
  362.             cell[i] |= CELL_SAMPLED;
  363.         }
  364.         {   const double s = sample[i];
  365.             if (result.minimum > s)
  366.                 result.minimum = s;
  367.             if (result.maximum < s)
  368.                 result.maximum = s;
  369.         }
  370.     }
  371.  
  372.     return result;
  373. }
  374.  
  375.  
  376. static int grid_minimum_cell(struct grid *const g)
  377. {
  378.     unsigned char *const cell = g->cell;
  379.     const size_t xstride = g->xstride,
  380.                  ystride = g->ystride,
  381.                  zstride = g->zstride,
  382.                  xsize = g->xsize,
  383.                  ysize = g->ysize,
  384.                  zsize = g->zsize,
  385.                  end = g->size;
  386.     size_t    x, y, z;
  387.     size_t    i, best_i;
  388.     double    s, best_s;
  389.  
  390.     grid_clear(g);
  391.  
  392.     best_s = grid_sample(g, 0, 0, 0);
  393.     best_i = 0;
  394.  
  395.     s = grid_sample(g, xsize-1, 0, 0);
  396.     if (s < best_s) {
  397.         best_s = s;
  398.         best_i = (xsize-1) * xstride;
  399.     }
  400.  
  401.     s = grid_sample(g, 0, ysize-1, 0);
  402.     if (s < best_s) {
  403.         best_s = s;
  404.         best_i = (ysize-1) * ystride;
  405.     }
  406.  
  407.     s = grid_sample(g, xsize-1, ysize-1, 0);
  408.     if (s < best_s) {
  409.         best_s = s;
  410.         best_i = (xsize-1) * xstride + (ysize-1) * ystride;
  411.     }
  412.  
  413.     s = grid_sample(g, 0, 0, zsize-1);
  414.     if (s < best_s) {
  415.         best_s = s;
  416.         best_i = (zsize-1) * zstride;
  417.     }
  418.  
  419.     s = grid_sample(g, xsize-1, 0, zsize-1);
  420.     if (s < best_s) {
  421.         best_s = s;
  422.         best_i = (xsize-1) * xstride + (zsize-1) * zstride;
  423.     }
  424.  
  425.     s = grid_sample(g, 0, ysize-1, zsize - 1);
  426.     if (s < best_s) {
  427.         best_s = s;
  428.         best_i = (ysize-1) * ystride + (zsize-1) * zstride;
  429.     }
  430.  
  431.     s = grid_sample(g, xsize-1, ysize-1, zsize-1);
  432.     if (s < best_s) {
  433.         best_s = s;
  434.         best_i = (xsize-1) * xstride + (ysize-1) * ystride + (zsize-1) * zstride;
  435.     }
  436.  
  437.     if (grid_push(g, best_i))
  438.         return errno;
  439.  
  440.     while (1) {
  441.         i = grid_pop(g);
  442.         if (i >= end)
  443.             break;
  444.  
  445.         x = (i / xstride) % xsize;
  446.         y = (i / ystride) % ysize;
  447.         z = (i / zstride) % zsize;
  448.         s = grid_sample(g, x, y, z);
  449.  
  450.         if (s < best_s || (s == best_s && i < best_i)) {
  451.             best_i = i;
  452.             best_s = s;
  453.         }
  454.  
  455.         if (x > 0 && grid_sample(g, x-1, y, z) <= best_s)
  456.             if (grid_push(g, i - xstride))
  457.                 return errno;
  458.  
  459.         if (y > 0 && grid_sample(g, x, y-1, z) <= best_s)
  460.             if (grid_push(g, i - ystride))
  461.                 return errno;
  462.  
  463.         if (z > 0 && grid_sample(g, x, y, z-1) <= best_s)
  464.             if (grid_push(g, i - zstride))
  465.                 return errno;
  466.  
  467.         if (x+1 < xsize && grid_sample(g, x+1, y, z) <= best_s)
  468.             if (grid_push(g, i + xstride))
  469.                 return errno;
  470.  
  471.         if (y+1 < ysize && grid_sample(g, x, y+1, z) <= best_s)
  472.             if (grid_push(g, i + ystride))
  473.                 return errno;
  474.  
  475.         if (z+1 < zsize && grid_sample(g, x, y, z+1) <= best_s)
  476.             if (grid_push(g, i + zstride))
  477.                 return errno;
  478.     }
  479.  
  480.     x = (best_i / xstride) % xsize;
  481.     if (x+1 >= xsize)
  482.         x = xsize - 2;
  483.  
  484.     y = (best_i / ystride) % ysize;
  485.     if (y+1 >= ysize)
  486.         y = ysize - 2;
  487.  
  488.     z = (best_i / zstride) % zsize;
  489.     if (z+1 >= zsize)
  490.         z = zsize - 2;
  491.  
  492.     g->minimum_cell = x*xstride + y*ystride + z*zstride;
  493.     g->minimum_sample = best_s;
  494.  
  495.     return 0;
  496. }
  497.  
  498.  
  499. static int grid_maximum_cell(struct grid *const g)
  500. {
  501.     unsigned char *const cell = g->cell;
  502.     const size_t xstride = g->xstride,
  503.                  ystride = g->ystride,
  504.                  zstride = g->zstride,
  505.                  xsize = g->xsize,
  506.                  ysize = g->ysize,
  507.                  zsize = g->zsize,
  508.                  end = g->size;
  509.     size_t    x, y, z;
  510.     size_t    i, best_i;
  511.     double    s, best_s;
  512.  
  513.     grid_clear(g);
  514.  
  515.     best_s = grid_sample(g, 0, 0, 0);
  516.     best_i = 0;
  517.  
  518.     s = grid_sample(g, xsize-1, 0, 0);
  519.     if (s >= best_s) {
  520.         best_s = s;
  521.         best_i = (xsize-1) * xstride;
  522.     }
  523.  
  524.     s = grid_sample(g, 0, ysize-1, 0);
  525.     if (s >= best_s) {
  526.         best_s = s;
  527.         best_i = (ysize-1) * ystride;
  528.     }
  529.  
  530.     s = grid_sample(g, xsize-1, ysize-1, 0);
  531.     if (s >= best_s) {
  532.         best_s = s;
  533.         best_i = (xsize-1) * xstride + (ysize-1) * ystride;
  534.     }
  535.  
  536.     s = grid_sample(g, 0, 0, zsize-1);
  537.     if (s >= best_s) {
  538.         best_s = s;
  539.         best_i = (zsize-1) * zstride;
  540.     }
  541.  
  542.     s = grid_sample(g, xsize-1, 0, zsize-1);
  543.     if (s >= best_s) {
  544.         best_s = s;
  545.         best_i = (xsize-1) * xstride + (zsize-1) * zstride;
  546.     }
  547.  
  548.     s = grid_sample(g, 0, ysize-1, zsize - 1);
  549.     if (s >= best_s) {
  550.         best_s = s;
  551.         best_i = (ysize-1) * ystride + (zsize-1) * zstride;
  552.     }
  553.  
  554.     s = grid_sample(g, xsize-1, ysize-1, zsize-1);
  555.     if (s >= best_s) {
  556.         best_s = s;
  557.         best_i = (xsize-1) * xstride + (ysize-1) * ystride + (zsize-1) * zstride;
  558.     }
  559.  
  560.     if (grid_push(g, best_i))
  561.         return errno;
  562.  
  563.     while (1) {
  564.         i = grid_pop(g);
  565.         if (i >= end)
  566.             break;
  567.  
  568.         x = (i / xstride) % xsize;
  569.         y = (i / ystride) % ysize;
  570.         z = (i / zstride) % zsize;
  571.         s = grid_sample(g, x, y, z);
  572.  
  573.         if (s > best_s || (s == best_s && i > best_i)) {
  574.             best_i = i;
  575.             best_s = s;
  576.         }
  577.  
  578.         if (x > 0 && grid_sample(g, x-1, y, z) >= best_s)
  579.             if (grid_push(g, i - xstride))
  580.                 return errno;
  581.  
  582.         if (y > 0 && grid_sample(g, x, y-1, z) >= best_s)
  583.             if (grid_push(g, i - ystride))
  584.                 return errno;
  585.  
  586.         if (z > 0 && grid_sample(g, x, y, z-1) >= best_s)
  587.             if (grid_push(g, i - zstride))
  588.                 return errno;
  589.  
  590.         if (x+1 < xsize && grid_sample(g, x+1, y, z) >= best_s)
  591.             if (grid_push(g, i + xstride))
  592.                 return errno;
  593.  
  594.         if (y+1 < ysize && grid_sample(g, x, y+1, z) >= best_s)
  595.             if (grid_push(g, i + ystride))
  596.                 return errno;
  597.  
  598.         if (z+1 < zsize && grid_sample(g, x, y, z+1) >= best_s)
  599.             if (grid_push(g, i + zstride))
  600.                 return errno;
  601.     }
  602.  
  603.     x = (best_i / xstride) % xsize;
  604.     if (x > 0)
  605.         x--;
  606.  
  607.     y = (best_i / ystride) % ysize;
  608.     if (y > 0)
  609.         y--;
  610.  
  611.     z = (best_i / zstride) % zsize;
  612.     if (z > 0)
  613.         z--;
  614.  
  615.     g->maximum_cell = x*xstride + y*ystride + z*zstride;
  616.     g->maximum_sample = best_s;
  617.  
  618.     return 0;
  619. }
  620.  
  621.  
  622. static size_t grid_find_cell(struct grid *const g, const double c)
  623. {
  624.     const size_t xsize = g->xsize,
  625.                  ysize = g->ysize,
  626.                  zsize = g->zsize,
  627.                  xstride = g->xstride,
  628.                  ystride = g->ystride,
  629.                  zstride = g->zstride;
  630.     size_t  xmin = (g->minimum_cell / xstride) % xsize,
  631.             ymin = (g->minimum_cell / ystride) % ysize,
  632.             zmin = (g->minimum_cell / zstride) % zsize;
  633.     size_t  xmax = (g->maximum_cell / xstride) % xsize,
  634.             ymax = (g->maximum_cell / ystride) % ysize,
  635.             zmax = (g->maximum_cell / zstride) % zsize;
  636.     span    s, smin, smax;
  637.  
  638.     smin = grid_cell_span(g, xmin, ymin, zmin);
  639.     smax = grid_cell_span(g, xmax, ymax, zmax);
  640.  
  641.     /* c outside the sample range? */
  642.     if (c < smin.minimum || c > smax.maximum)
  643.         return g->size; /* Not found */
  644.  
  645.     /* c within the minimum cell? */
  646.     if (smin.minimum <= c && c <= smin.maximum)
  647.         return g->minimum_cell;
  648.  
  649.     /* c within the maximum cell? */
  650.     if (smax.maximum <= c && c <= smax.maximum)
  651.         return g->maximum_cell;
  652.  
  653.     /* Do a binary search to locate a cell spanning c. */
  654.     while (1) {
  655.         const size_t x = (xmin + xmax) / 2;
  656.         const size_t y = (ymin + ymax) / 2;
  657.         const size_t z = (zmin + zmax) / 2;
  658.  
  659.         s = grid_cell_span(g, x, y, z);
  660.         if (s.minimum <= c && c <= s.maximum)
  661.             return x * xstride + y * ystride + z * zstride;
  662.  
  663.         if (x == xmin && y == ymin && z == zmin)
  664.             return xmax * xstride + ymax * ystride + zmax * zstride;
  665.  
  666.         if (s.minimum > c) {
  667.             smax = s;
  668.             xmax = x;
  669.             ymax = y;
  670.             zmax = z;
  671.             continue;
  672.         }
  673.  
  674.         if (s.maximum < c) {
  675.             smin = s;
  676.             xmin = x;
  677.             ymin = y;
  678.             zmin = z;
  679.             continue;
  680.         }
  681.  
  682.         if (s.maximum > c) {
  683.             smax = s;
  684.             xmax = x;
  685.             ymax = y;
  686.             zmax = z;
  687.             continue;
  688.         }
  689.  
  690.         if (s.minimum < c) {
  691.             smin = s;
  692.             xmin = x;
  693.             ymin = y;
  694.             zmin = z;
  695.             continue;
  696.         }
  697.     }
  698.  
  699.     /* Only possibility! */
  700.     return xmax * xstride + ymax * ystride + zmax * zstride;
  701. }
  702.  
  703. static int grid_isosurface(struct grid *const g, const double value, void *custom,
  704.                            int (*callback)(void *custom, struct grid *g, int x, int y, int z, double value, double sample[2][2][2]))
  705. {
  706.     double corner[2][2][2];
  707.     size_t i, x, y, z;
  708.     span s;
  709.     int err;
  710.  
  711.     if (g == NULL || callback == NULL)
  712.         return errno = EINVAL;
  713.  
  714.     /* Use a binary search from minimum to maximum cells,
  715.      * to find a suitable starting cell. */    
  716.     i = grid_find_cell(g, value);
  717.     if (i >= g->size)
  718.         return errno = ENOENT;
  719.  
  720.     grid_clear(g);
  721.  
  722.     err = grid_push(g, i);
  723.     if (err)
  724.         return errno = err;
  725.  
  726.     while (1) {
  727.  
  728.         i = grid_pop(g);
  729.         if (i >= g->size)
  730.             break;
  731.  
  732.         x = (i / g->xstride) % g->xsize;
  733.         y = (i / g->ystride) % g->ysize;
  734.         z = (i / g->zstride) % g->zsize;
  735.         s = grid_cell(g, x, y, z, corner);
  736.  
  737.         if (s.minimum <= value && s.maximum >= value) {
  738.             err = callback(custom, g, x, y, z, value, corner);
  739.             if (err) {
  740.                 errno = 0;
  741.                 return err;
  742.             }
  743.         }
  744.  
  745.         /* Neighboring cells sharing a face */
  746.  
  747.         if (x > 0 && !(g->cell[i - g->xstride] & CELL_PUSHED)) {
  748.             s = grid_cell_span(g, x-1, y, z);
  749.             if (s.minimum <= value && s.maximum >= value)
  750.                 if (grid_push(g, i - g->xstride))
  751.                     return errno;
  752.         }
  753.  
  754.         if (y > 0 && !(g->cell[i - g->ystride] & CELL_PUSHED)) {
  755.             s = grid_cell_span(g, x, y-1, z);
  756.             if (s.minimum <= value && s.maximum >= value)
  757.                 if (grid_push(g, i - g->ystride))
  758.                     return errno;
  759.         }
  760.  
  761.         if (z > 0 && !(g->cell[i - g->zstride] & CELL_PUSHED)) {
  762.             s = grid_cell_span(g, x, y, z-1);
  763.             if (s.minimum <= value && s.maximum >= value)
  764.                 if (grid_push(g, i - g->zstride))
  765.                     return errno;
  766.         }
  767.  
  768.         if (x < g->xsize-2 && !(g->cell[i + g->xstride] & CELL_PUSHED)) {
  769.             s = grid_cell_span(g, x+1, y, z);
  770.             if (s.minimum <= value && s.maximum >= value)
  771.                 if (grid_push(g, i + g->xstride))
  772.                     return errno;
  773.         }
  774.  
  775.         if (y < g->ysize-2 && !(g->cell[i + g->ystride] & CELL_PUSHED)) {
  776.             s = grid_cell_span(g, x, y+1, z);
  777.             if (s.minimum <= value && s.maximum >= value)
  778.                 if (grid_push(g, i + g->ystride))
  779.                     return errno;
  780.         }
  781.  
  782.         if (z < g->zsize-2 && !(g->cell[i + g->zstride] & CELL_PUSHED)) {
  783.             s = grid_cell_span(g, x, y, z+1);
  784.             if (s.minimum <= value && s.maximum >= value)
  785.                 if (grid_push(g, i + g->zstride))
  786.                     return errno;
  787.         }
  788.     }
  789.  
  790.     /* Success. */
  791.     return errno = 0;
  792. }
  793.  
  794.  
  795. static struct grid *grid_free(struct grid *const g)
  796. {
  797.     if (g != NULL) {
  798.         free(g->stack);
  799.  
  800.         g->xsize = 0;
  801.         g->ysize = 0;
  802.         g->zsize = 0;
  803.  
  804.         g->size = 0;
  805.  
  806.         g->xstride = 0;
  807.         g->ystride = 0;
  808.         g->zstride = 0;
  809.  
  810.         g->minimum_cell = 0;
  811.         g->maximum_cell = 0;
  812.         g->minimum_sample = 0.0;
  813.         g->maximum_sample = 0.0;
  814.  
  815.         g->xorigin = 0.0;
  816.         g->yorigin = 0.0;
  817.         g->zorigin = 0.0;
  818.  
  819.         g->xunit = 0.0;
  820.         g->yunit = 0.0;
  821.         g->zunit = 0.0;
  822.  
  823.         g->field_data = NULL;
  824.         g->field = NULL;
  825.  
  826.         g->stack_size = 0;
  827.         g->stack_used = 0;
  828.         g->stack = NULL;
  829.  
  830.         g->cell = NULL;
  831.  
  832.         free(g);
  833.     }
  834.     return NULL;
  835. }
  836.  
  837.  
  838. static struct grid *grid_new(const int xsize, const int ysize, const int zsize,
  839.                              const double xmin, const double ymin, const double zmin,
  840.                              const double xmax, const double ymax, const double zmax,
  841.                              void *field_data,
  842.                              double (*field)(void *data, double x, double y, double z))
  843. {
  844.     const size_t size = (size_t)xsize * (size_t)ysize * (size_t)zsize;
  845.     const size_t base = sizeof (struct grid) + size * sizeof (double);
  846.     const size_t bytes = base + size;
  847.     struct grid *g;
  848.     size_t       i;
  849.     int          err;
  850.  
  851.     /* Make sure parameters make sense. */
  852.     if (xsize < 2 || ysize < 2 || zsize < 2 || field == NULL) {
  853.         errno = EINVAL;
  854.         return NULL;
  855.     }
  856.  
  857.     /* Try to detect size overflow. */
  858.     if (bytes <= base || bytes <= size ||
  859.         (size_t)((base - sizeof (struct grid)) / sizeof (double)) != size ||
  860.         (size_t)(size / (size_t)xsize) != (size_t)ysize * (size_t)zsize ||
  861.         (size_t)(size / (size_t)ysize / (size_t)zsize) != (size_t)xsize) {
  862.         errno = ENOMEM;
  863.         return NULL;
  864.     }
  865.  
  866.     /* Allocate; both cell and sample arrays are part of the structure allocation. */
  867.     g = malloc(bytes);
  868.     if (g == NULL) {
  869.         errno = ENOMEM;
  870.         return NULL;
  871.     }
  872.  
  873.     /* Initialize the fields. */
  874.     g->xsize = xsize;
  875.     g->ysize = ysize;
  876.     g->zsize = zsize;
  877.  
  878.     g->size = size;
  879.  
  880.     g->xstride = 1;
  881.     g->ystride = (size_t)xsize;
  882.     g->zstride = (size_t)xsize * (size_t)ysize;
  883.  
  884.     g->xorigin = xmin;
  885.     g->yorigin = ymin;
  886.     g->zorigin = zmin;
  887.  
  888.     g->xunit = (xmax - xmin) / (double)(xsize - 1);
  889.     g->yunit = (ymax - ymin) / (double)(ysize - 1);
  890.     g->zunit = (zmax - zmin) / (double)(zsize - 1);
  891.  
  892.     g->field_data = field_data;
  893.     g->field = field;
  894.  
  895.     g->stack_size = 0;
  896.     g->stack_used = 0;
  897.     g->stack = NULL;
  898.  
  899.     /* Cells are allocated as part of the structure, following the last cached sample. */
  900.     g->cell = (unsigned char *)&(g->sample[size]);
  901.  
  902.     /* Clear sample cache. This is actually not really necessary, just "nice". */
  903.     g->sample[0] = 0.0;
  904.     for (i = 1; i < size; i++)
  905.         memcpy(g->sample + i, g->sample, sizeof g->sample[0]);
  906.  
  907.     /* Clear cells to initial state. This is actually necessary. */
  908.     memset(g->cell, CELL_UNKNOWN, size);
  909.  
  910.     /* Locate the cell containing the minimum sample. */
  911.     err = grid_minimum_cell(g);
  912.     if (err) {
  913.         free(g);
  914.         errno = err;
  915.         return NULL;
  916.     }
  917.  
  918.     /* Locate the cell containing the maximum sample. */
  919.     err = grid_maximum_cell(g);
  920.     if (err) {
  921.         free(g);
  922.         errno = err;
  923.         return NULL;
  924.     }
  925.  
  926.     /* All done successfully. */
  927.     errno = 0;
  928.     return g;
  929. }
  930.  
  931. #endif /* GRID_H */
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!