1 /* GNUPLOT - voxelgrid.c */
2
3 /*[
4 * Copyright Ethan A Merritt 2019
5 *
6 * Gnuplot license:
7 *
8 * Permission to use, copy, and distribute this software and its
9 * documentation for any purpose with or without fee is hereby granted,
10 * provided that the above copyright notice appear in all copies and
11 * that both that copyright notice and this permission notice appear
12 * in supporting documentation.
13 *
14 * Permission to modify the software is granted, but not the right to
15 * distribute the complete modified source code. Modifications are to
16 * be distributed as patches to the released version. Permission to
17 * distribute binaries produced by compiling modified sources is granted,
18 * provided you
19 * 1. distribute the corresponding source modifications from the
20 * released version in the form of a patch file along with the binaries,
21 * 2. add special version identification to distinguish your version
22 * in addition to the base release version number,
23 * 3. provide your name and address as the primary contact for the
24 * support of your modified version, and
25 * 4. retain our contact information in regard to use of the base
26 * software.
27 * Permission to distribute the released version of the source code along
28 * with corresponding source modifications in the form of a patch file is
29 * granted with same provisions 2 through 4 for binary distributions.
30 *
31 * This software is provided "as is" without express or implied warranty
32 * to the extent permitted by applicable law.
33 *
34 * Alternative license:
35 *
36 * As an alternative to distributing code in this file under the gnuplot license,
37 * you may instead comply with the terms below. In this case, redistribution and
38 * use in source and binary forms, with or without modification, are permitted
39 * provided that the following conditions are met:
40 *
41 * Redistributions of source code must retain the above copyright notice, this
42 * list of conditions and the following disclaimer. Redistributions in binary
43 * form must reproduce the above copyright notice, this list of conditions and
44 * the following disclaimer in the documentation and/or other materials provided
45 * with the distribution.
46 *
47 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
48 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
49 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
50 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
51 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
52 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
53 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
54 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
55 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
56 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
57 * POSSIBILITY OF SUCH DAMAGE.
58 *
59 ]*/
60
61 /*
62 * This file implements a 3D grid of NxNxN evenly spaced voxels.
63 * The command
64 * set vgrid size N
65 * frees any previous grid, allocates a new one, and initializes it to zero.
66 * The grid range can be set explicitly using
67 * set vxrange [vxmin:vxmax]
68 * set vyrange [vymin:vymax]
69 * set vzrange [vzmin:vzmax]
70 * otherwise it inherits the xrange, yrange, and zrange limits present at the
71 * time of the first "vclear", "vfill", or "voxel() = " command.
72 *
73 * Grid voxels can be filled individually using the command
74 * voxel(x,y,z) = FOO
75 * or filled in the neighborhood of a data point using the command
76 * vfill $FILE using x:y:z:radius:(function)
77 * In this case the nearest grid point and all other grid points within
78 * a sphere centered about [x,y,z] are incremented by
79 * voxel(vx,vy,vz) += function( sqrt((x-vx)**2 + (y-vy)**2 + (z-vz)**2) )
80 * Once loaded the grid can either be referenced by splot commands, e.g.
81 * splot $FILE using x:y:z:(voxel(x,y,z)) with points lc palette
82 * or plotted using new splot options
83 * splot $GRID with {dots|points}
84 * splot $GRID with isosurface level <value>
85 *
86 * The grid can be cleared by
87 * vclear
88 * or deallocated by
89 * unset vgrid
90 */
91
92 #include "gp_types.h"
93 #include "alloc.h"
94 #include "axis.h"
95 #include "command.h"
96 #include "datablock.h"
97 #include "datafile.h"
98 #include "eval.h"
99 #include "graphics.h"
100 #include "graph3d.h"
101 #include "parse.h"
102 #include "util.h"
103 #include "variable.h"
104 #include "voxelgrid.h"
105
106 #ifdef VOXEL_GRID_SUPPORT
107
108 static vgrid *current_vgrid = NULL; /* active voxel grid */
109 static struct udvt_entry *udv_VoxelDistance = NULL; /* reserved user variable */
110 struct isosurface_opt isosurface_options;
111
112 /* Internal prototypes */
113 static void vfill( t_voxel *grid );
114 static void modify_voxels( t_voxel *grid, double x, double y, double z,
115 double radius, struct at_type *function );
116 static void vgrid_stats( vgrid *vgrid );
117
118 /* Purely local bookkeeping */
119 static int nvoxels_modified;
120 static struct at_type *density_function = NULL;
121
122 /*
123 * called on program entry and by "reset session"
124 */
125 void
init_voxelsupport()126 init_voxelsupport()
127 {
128 /* Make sure there is a user variable that can be used as a parameter
129 * to the function in the 5th spec of "vfill".
130 * Scripts can test if (exists("VoxelDistance")) to check for voxel support.
131 */
132 udv_VoxelDistance = add_udv_by_name("VoxelDistance");
133 udv_VoxelDistance->udv_value.type = CMPLX;
134
135 /* default state of other voxel-related structures */
136 isosurface_options.inside_offset = 1; /* inside color = outside + 1 */
137 isosurface_options.tessellation = 0; /* mixed triangles + quadrangles */
138 }
139
140 /*
141 * vx vy vz ranges must be established before the grid can be used
142 */
143 void
check_grid_ranges()144 check_grid_ranges()
145 {
146 if (current_vgrid == NULL)
147 int_error(NO_CARET, "vgrid must be set before use");
148
149 if (isnan(current_vgrid->vxmin) || isnan(current_vgrid->vxmax)) {
150 if ((axis_array[FIRST_X_AXIS].set_autoscale & AUTOSCALE_BOTH) == AUTOSCALE_NONE) {
151 current_vgrid->vxmin = axis_array[FIRST_X_AXIS].set_min;
152 current_vgrid->vxmax = axis_array[FIRST_X_AXIS].set_max;
153 } else
154 int_error(NO_CARET, "grid limits must be set before use");
155 }
156 if (isnan(current_vgrid->vymin) || isnan(current_vgrid->vymax)) {
157 if ((axis_array[FIRST_Y_AXIS].set_autoscale & AUTOSCALE_BOTH) == AUTOSCALE_NONE) {
158 current_vgrid->vymin = axis_array[FIRST_Y_AXIS].set_min;
159 current_vgrid->vymax = axis_array[FIRST_Y_AXIS].set_max;
160 } else
161 int_error(NO_CARET, "grid limits must be set before use");
162 }
163 if (isnan(current_vgrid->vzmin) || isnan(current_vgrid->vzmax)) {
164 if ((axis_array[FIRST_Z_AXIS].set_autoscale & AUTOSCALE_BOTH) == AUTOSCALE_NONE) {
165 current_vgrid->vzmin = axis_array[FIRST_Z_AXIS].set_min;
166 current_vgrid->vzmax = axis_array[FIRST_Z_AXIS].set_max;
167 } else
168 int_error(NO_CARET, "grid limits must be set before use");
169 }
170
171 current_vgrid->vxdelta = (current_vgrid->vxmax - current_vgrid->vxmin)
172 / (current_vgrid->size - 1);
173 current_vgrid->vydelta = (current_vgrid->vymax - current_vgrid->vymin)
174 / (current_vgrid->size - 1);
175 current_vgrid->vzdelta = (current_vgrid->vzmax - current_vgrid->vzmin)
176 / (current_vgrid->size - 1);
177 }
178
179 /*
180 * Initialize vgrid array
181 * set vgrid <name> {size <N>}
182 * - retrieve existing voxel grid or create a new one
183 * - initialize to N (defaults to N=100 for a new grid)
184 * - set current_vgrid to this grid.
185 */
186 void
set_vgrid()187 set_vgrid()
188 {
189 struct udvt_entry *grid = NULL;
190 int new_size = 100; /* default size */
191 char *name;
192
193 c_token++;
194 if (END_OF_COMMAND || !isletter(c_token+1))
195 int_error(c_token, "syntax: set vgrid $<gridname> {size N}");
196
197 /* Create or recycle a datablock with the requested name */
198 name = parse_datablock_name();
199 grid = add_udv_by_name(name);
200 if (grid->udv_value.type == VOXELGRID) {
201 /* Keep size of existing grid */
202 new_size = grid->udv_value.v.vgrid->size;
203 current_vgrid = grid->udv_value.v.vgrid;
204 } else {
205 /* Note: The only other variable type that starts with a $ is DATABLOCK */
206 free_value(&grid->udv_value);
207 current_vgrid = gp_alloc(sizeof(vgrid), "new vgrid");
208 memset(current_vgrid, 0, sizeof(vgrid));
209 current_vgrid->vxmin = not_a_number();
210 current_vgrid->vxmax = not_a_number();
211 current_vgrid->vymin = not_a_number();
212 current_vgrid->vymax = not_a_number();
213 current_vgrid->vzmin = not_a_number();
214 current_vgrid->vzmax = not_a_number();
215 grid->udv_value.v.vgrid = current_vgrid;
216 grid->udv_value.type = VOXELGRID;
217 }
218
219 if (equals(c_token, "size")) {
220 c_token++;
221 new_size = int_expression();
222 }
223
224 /* FIXME: arbitrary limit to reduce chance of memory exhaustion */
225 if (new_size < 10 || new_size > 256)
226 int_error(NO_CARET, "vgrid size must be an integer between 10 and 256");
227
228 /* Initialize storage for new voxel grid */
229 if (current_vgrid->size != new_size) {
230 current_vgrid->size = new_size;
231 current_vgrid->vdata = gp_realloc(current_vgrid->vdata,
232 new_size*new_size*new_size*sizeof(t_voxel), "voxel array");
233 memset(current_vgrid->vdata, 0, new_size*new_size*new_size*sizeof(t_voxel));
234 }
235
236 }
237
238 /*
239 * set vxrange [min:max]
240 */
241 void
set_vgrid_range()242 set_vgrid_range()
243 {
244 double gridmin, gridmax;
245 int save_token = c_token++;
246
247 if (!current_vgrid)
248 int_error(NO_CARET, "no voxel grid is active");
249
250 if (!equals(c_token,"["))
251 return;
252 c_token++;
253 gridmin = real_expression();
254 if (!equals(c_token,":"))
255 return;
256 c_token++;
257 gridmax = real_expression();
258 if (!equals(c_token,"]"))
259 return;
260 c_token++;
261
262 if (almost_equals(save_token, "vxr$ange")) {
263 current_vgrid->vxmin = gridmin;
264 current_vgrid->vxmax = gridmax;
265 }
266 if (almost_equals(save_token, "vyr$ange")) {
267 current_vgrid->vymin = gridmin;
268 current_vgrid->vymax = gridmax;
269 }
270 if (almost_equals(save_token, "vzr$ange")) {
271 current_vgrid->vzmin = gridmin;
272 current_vgrid->vzmax = gridmax;
273 }
274 }
275
276 /*
277 * show state of all voxel grids
278 */
279 void
show_vgrid()280 show_vgrid()
281 {
282 struct udvt_entry *udv;
283 vgrid *vgrid;
284
285 for (udv = first_udv; udv != NULL; udv = udv->next_udv) {
286 if (udv->udv_value.type == VOXELGRID) {
287 vgrid = udv->udv_value.v.vgrid;
288
289 fprintf(stderr, "\t%s:", udv->udv_name);
290 if (vgrid == current_vgrid)
291 fprintf(stderr, "\t(active)");
292 fprintf(stderr, "\n");
293 fprintf(stderr, "\t\tsize %d X %d X %d\n",
294 vgrid->size, vgrid->size, vgrid->size);
295 if (isnan(vgrid->vxmin) || isnan(vgrid->vxmax) || isnan(vgrid->vymin)
296 || isnan(vgrid->vymax) || isnan(vgrid->vzmin) || isnan(vgrid->vzmax)) {
297 fprintf(stderr, "\t\tgrid ranges not set\n");
298 continue;
299 }
300 fprintf(stderr, "\t\tvxrange [%g:%g] vyrange[%g:%g] vzrange[%g:%g]\n",
301 vgrid->vxmin, vgrid->vxmax, vgrid->vymin,
302 vgrid->vymax, vgrid->vzmin, vgrid->vzmax);
303 vgrid_stats(vgrid);
304 fprintf(stderr, "\t\tnon-zero voxel values: min %.2g max %.2g mean %.2g stddev %.2g\n",
305 vgrid->min_value, vgrid->max_value, vgrid->mean_value, vgrid->stddev);
306 fprintf(stderr, "\t\tnumber of zero voxels: %d (%.2f%%)\n",
307 vgrid->nzero,
308 100. * (double)vgrid->nzero / (vgrid->size*vgrid->size*vgrid->size));
309 }
310 }
311 }
312
313 /*
314 * run through the whole grid
315 * accumulate min/max, mean, and standard deviation of non-zero voxels
316 * TODO: median
317 */
318 static void
vgrid_stats(vgrid * vgrid)319 vgrid_stats(vgrid *vgrid)
320 {
321 double min = VERYLARGE;
322 double max = -VERYLARGE;
323 double sum = 0;
324 int nzero = 0;
325 t_voxel *voxel;
326 int N = vgrid->size;
327 int i;
328
329 /* bookkeeping for standard deviation */
330 double num = 0;
331 double delta = 0;
332 double delta2 = 0;
333 double mean = 0;
334 double mean2 = 0;
335
336 for (voxel = vgrid->vdata, i=0; i < N*N*N; voxel++, i++) {
337 if (*voxel == 0) {
338 nzero++;
339 continue;
340 }
341 sum += *voxel;
342 if (min > *voxel)
343 min = *voxel;
344 if (max < *voxel)
345 max = *voxel;
346
347 /* standard deviation */
348 num += 1.0;
349 delta = *voxel - mean;
350 mean += delta/num;
351 delta2 = *voxel - mean;
352 mean2 += delta * delta2;
353 }
354
355 vgrid->min_value = min;
356 vgrid->max_value = max;
357 vgrid->nzero = nzero;
358 if (num < 2) {
359 vgrid->mean_value = vgrid->stddev = not_a_number();
360 } else {
361 vgrid->mean_value = sum / (double)(N*N*N - nzero);
362 vgrid->stddev = sqrt(mean2 / (num - 1.0));
363 }
364
365 /* all zeros */
366 if (nzero == N*N*N) {
367 vgrid->min_value = 0;
368 vgrid->max_value = 0;
369 }
370 }
371
372 udvt_entry *
get_vgrid_by_name(char * name)373 get_vgrid_by_name(char *name)
374 {
375 struct udvt_entry *vgrid = get_udv_by_name(name);
376
377 if (!vgrid || vgrid->udv_value.type != VOXELGRID)
378 return NULL;
379 else
380 return vgrid;
381 }
382
383 /*
384 * initialize content of voxel grid to all zero
385 */
386 void
vclear_command()387 vclear_command()
388 {
389 vgrid *vgrid = current_vgrid;
390
391 c_token++;
392 if (!END_OF_COMMAND && equals(c_token, "$")) {
393 char *name = parse_datablock_name();
394 udvt_entry *vgrid_udv = get_vgrid_by_name(name);
395 if (!vgrid_udv)
396 int_error(c_token, "no such voxel grid");
397 vgrid = vgrid_udv->udv_value.v.vgrid;
398 }
399 if (vgrid && vgrid->size && vgrid->vdata) {
400 int size = vgrid->size;
401 memset(vgrid->vdata, 0, size*size*size * sizeof(t_voxel));
402 }
403 }
404
405 /*
406 * deallocate storage for a voxel grid
407 */
408 void
gpfree_vgrid(struct udvt_entry * grid)409 gpfree_vgrid(struct udvt_entry *grid)
410 {
411 if (grid->udv_value.type != VOXELGRID)
412 return;
413 free(grid->udv_value.v.vgrid->vdata);
414 free(grid->udv_value.v.vgrid);
415 if (grid->udv_value.v.vgrid == current_vgrid)
416 current_vgrid = NULL;
417 grid->udv_value.v.vgrid = NULL;
418 grid->udv_value.type = NOTDEFINED;
419 }
420
421 /*
422 * 'unset $vgrid' command
423 */
424 void
unset_vgrid()425 unset_vgrid()
426 {
427 struct udvt_entry *grid = NULL;
428 char *name;
429
430 if (END_OF_COMMAND || !equals(c_token,"$"))
431 int_error(c_token, "syntax: unset vgrid $<gridname>");
432
433 /* Look for a datablock with the requested name */
434 name = parse_datablock_name();
435 grid = get_vgrid_by_name(name);
436 if (!grid)
437 int_error(c_token, "no such vgrid");
438
439 gpfree_vgrid(grid);
440 }
441
442 /*
443 * "set isosurface {triangles|mixed}"
444 */
445 void
set_isosurface()446 set_isosurface()
447 {
448 while (!END_OF_COMMAND) {
449 c_token++;
450 if (almost_equals(c_token, "triang$les")) {
451 c_token++;
452 isosurface_options.tessellation = 1;
453 } else if (almost_equals(c_token, "mix$ed")) {
454 c_token++;
455 isosurface_options.tessellation = 0;
456 } else if (almost_equals(c_token, "inside$color")) {
457 c_token++;
458 if (END_OF_COMMAND)
459 isosurface_options.inside_offset = 1;
460 else
461 isosurface_options.inside_offset = int_expression();
462 } else if (almost_equals(c_token, "noin$sidecolor")) {
463 c_token++;
464 isosurface_options.inside_offset = 0;
465 } else
466 int_error(c_token,"unrecognized option");
467 }
468 }
469
470 void
show_isosurface()471 show_isosurface()
472 {
473 c_token++;
474 fprintf(stderr,"\tisosurfaces will use %s\n",
475 isosurface_options.tessellation != 0 ? "triangles only"
476 : "a mixture of triangles and quadrangles");
477 fprintf(stderr,"\tinside surface linetype offset by %d\n",
478 isosurface_options.inside_offset);
479 }
480
481 /*
482 * voxel(x,y,z) = expr()
483 */
484 void
voxel_command()485 voxel_command()
486 {
487 double vx, vy, vz;
488 int ivx, ivy, ivz;
489 int N;
490 t_voxel *voxel;
491
492 check_grid_ranges();
493
494 c_token++;
495 if (!equals(c_token++,"("))
496 int_error(c_token-1, "syntax: voxel(x,y,z) = newvalue");
497 vx = real_expression();
498 if (!equals(c_token++,","))
499 int_error(c_token-1, "syntax: voxel(x,y,z) = newvalue");
500 vy = real_expression();
501 if (!equals(c_token++,","))
502 int_error(c_token-1, "syntax: voxel(x,y,z) = newvalue");
503 vz = real_expression();
504 if (!equals(c_token++,")"))
505 int_error(c_token-1, "syntax: voxel(x,y,z) = newvalue");
506 if (!equals(c_token++,"="))
507 int_error(c_token-1, "syntax: voxel(x,y,z) = newvalue");
508
509 if (vx < current_vgrid->vxmin || vx > current_vgrid->vxmax
510 || vy < current_vgrid->vymin || vy > current_vgrid->vymax
511 || vz < current_vgrid->vzmin || vz > current_vgrid->vzmax) {
512 int_warn(NO_CARET, "voxel out of range");
513 (void)real_expression();
514 return;
515 }
516
517 ivx = ceil( (vx - current_vgrid->vxmin) / current_vgrid->vxdelta );
518 ivy = ceil( (vy - current_vgrid->vymin) / current_vgrid->vydelta );
519 ivz = ceil( (vz - current_vgrid->vzmin) / current_vgrid->vzdelta );
520
521 N = current_vgrid->size;
522 FPRINTF((stderr,"\tvgrid array index = %d\n", ivx + ivy * N + ivz * N*N));
523
524 voxel = ¤t_vgrid->vdata[ ivx + ivy * N + ivz * N*N ];
525 *voxel = real_expression();
526
527 return;
528 }
529
530 /*
531 * internal look-up function voxel(x,y,z)
532 */
533 t_voxel
voxel(double vx,double vy,double vz)534 voxel( double vx, double vy, double vz )
535 {
536 int ivx, ivy, ivz;
537 int N;
538
539 if (!current_vgrid)
540 return not_a_number();
541
542 if (vx < current_vgrid->vxmin || vx > current_vgrid->vxmax
543 || vy < current_vgrid->vymin || vy > current_vgrid->vymax
544 || vz < current_vgrid->vzmin || vz > current_vgrid->vzmax)
545 return not_a_number();
546
547 ivx = ceil( (vx - current_vgrid->vxmin) / current_vgrid->vxdelta );
548 ivy = ceil( (vy - current_vgrid->vymin) / current_vgrid->vydelta );
549 ivz = ceil( (vz - current_vgrid->vzmin) / current_vgrid->vzdelta );
550
551 N = current_vgrid->size;
552 return current_vgrid->vdata[ ivx + ivy * N + ivz * N*N ];
553 }
554
555 /*
556 * user-callable retrieval function voxel(x,y,z)
557 */
558 void
f_voxel(union argument * arg)559 f_voxel(union argument *arg)
560 {
561 struct value a;
562 double vx, vy, vz;
563
564 (void) arg; /* avoid -Wunused warning */
565 vz = real(pop(&a));
566 vy = real(pop(&a));
567 vx = real(pop(&a));
568
569 if (!current_vgrid)
570 int_error(NO_CARET, "no active voxel grid");
571
572 if (vx < current_vgrid->vxmin || vx > current_vgrid->vxmax
573 || vy < current_vgrid->vymin || vy > current_vgrid->vymax
574 || vz < current_vgrid->vzmin || vz > current_vgrid->vzmax) {
575 push(&(udv_NaN->udv_value));
576 return;
577 }
578
579 push( Gcomplex(&a, voxel(vx, vy, vz), 0.0) );
580 return;
581 }
582
583 /*
584 * "vfill" works very much like "plot" in that it reads from an input stream
585 * of data points specified by the same "using" "skip" "every" and other keyword
586 * syntax shared with the "plot" "splot" and "stats" commands.
587 * Basic example:
588 * vfill $FILE using x:y:z:radius:(function())
589 * However instead of creating a plot immediately, vfill modifies the content
590 * of a voxel grid by incrementing the value of each voxel within a radius.
591 * voxel(vx,vy,vz) += function( distance([vx,vy,vz], [x,y,z]) )
592 *
593 * vfill shares the routines df_open df_readline ... with the plot and splot
594 * commands.
595 */
596 void
vfill_command()597 vfill_command()
598 {
599 c_token++;
600 vfill(current_vgrid->vdata);
601 }
602
603 static void
vfill(t_voxel * grid)604 vfill(t_voxel *grid)
605 {
606 int plot_num = 0;
607 struct curve_points dummy_plot;
608 struct curve_points *this_plot = &dummy_plot;
609 t_value original_value_sample_var;
610 char *name_str;
611
612 check_grid_ranges();
613
614 /*
615 * This part is modelled on eval_plots()
616 */
617 plot_iterator = check_for_iteration();
618 while (TRUE) {
619 int sample_range_token;
620 int start_token = c_token;
621 int specs;
622
623 /* Forgive trailing comma on a multi-element plot command */
624 if (END_OF_COMMAND) {
625 if (plot_num == 0)
626 int_error(c_token, "data input source expected");
627 break;
628 }
629 plot_num++;
630
631 /* Check for a sampling range */
632 if (equals(c_token,"sample") && equals(c_token+1,"["))
633 c_token++;
634 sample_range_token = parse_range(SAMPLE_AXIS);
635 if (sample_range_token != 0)
636 axis_array[SAMPLE_AXIS].range_flags |= RANGE_SAMPLED;
637
638 /* FIXME: allow replacement of dummy variable? */
639
640 name_str = string_or_express(NULL);
641 if (!name_str)
642 int_error(c_token, "no input data source");
643 if (!strcmp(name_str, "@@"))
644 int_error(c_token-1, "filling from array not supported");
645 if (sample_range_token !=0 && *name_str != '+')
646 int_warn(sample_range_token, "Ignoring sample range in non-sampled data plot");
647
648 /* Dummy up a plot structure so that we can share code with plot command */
649 memset( this_plot, 0, sizeof(struct curve_points));
650 /* FIXME: what exactly do we need to fill in? */
651 this_plot->plot_type = DATA;
652 this_plot->noautoscale = TRUE;
653
654 /* Fixed number of input columns x:y:z:radius:(density_function) */
655 specs = df_open(name_str, 5, this_plot);
656
657 /* We will invoke density_function in modify_voxels rather than df_readline */
658 if (use_spec[4].at == NULL)
659 int_error(NO_CARET, "5th user spec to vfill must be an expression");
660 else {
661 free_at(density_function);
662 density_function = use_spec[4].at;
663 use_spec[4].at = NULL;
664 use_spec[4].column = 0;
665 }
666
667 /* Initialize user variable VoxelDistance used by density_function() */
668 Gcomplex(&udv_VoxelDistance->udv_value, 0.0, 0.0);
669
670 /* Store a pointer to the named variable used for sampling */
671 /* Save prior value of sample variables so we can restore them later */
672 this_plot->sample_var = add_udv_by_name(c_dummy_var[0]);
673 original_value_sample_var = this_plot->sample_var->udv_value;
674 this_plot->sample_var->udv_value.type = NOTDEFINED;
675
676 /* We don't support any further options */
677 if (almost_equals(c_token, "w$ith"))
678 int_error(c_token, "vfill does not support 'with' options");
679
680 /* This part is modelled on get_data().
681 * However we process each point as we come to it.
682 */
683 if (df_no_use_specs == 5) {
684 int j;
685 int ngood;
686 double v[MAXDATACOLS];
687 memset(v, 0, sizeof(v));
688
689 /* Initialize stats */
690 ngood = 0;
691 nvoxels_modified = 0;
692 fprintf(stderr,"vfill from %s :\n", name_str);
693
694 /* If the user has set an explicit locale for numeric input, apply it */
695 /* here so that it affects data fields read from the input file. */
696 set_numeric_locale();
697
698 /* Initial state */
699 df_warn_on_missing_columnheader = TRUE;
700
701 while ((j = df_readline(v, 5)) != DF_EOF) {
702 switch (j) {
703
704 case 0:
705 df_close();
706 int_error(this_plot->token, "Bad data on line %d of file %s",
707 df_line_number, df_filename ? df_filename : "");
708 continue;
709
710 case DF_UNDEFINED:
711 case DF_MISSING:
712 continue;
713
714 case DF_FIRST_BLANK:
715 case DF_SECOND_BLANK:
716 continue;
717
718 case DF_COLUMN_HEADERS:
719 case DF_FOUND_KEY_TITLE:
720 case DF_KEY_TITLE_MISSING:
721 continue;
722
723 default:
724 break; /* Not continue! */
725 }
726
727 /* Ignore out of range points */
728 /* FIXME: probably should be range + radius! */
729 if (v[0] < current_vgrid->vxmin || current_vgrid->vxmax < v[0])
730 continue;
731 if (v[1] < current_vgrid->vymin || current_vgrid->vymax < v[1])
732 continue;
733 if (v[2] < current_vgrid->vzmin || current_vgrid->vzmax < v[2])
734 continue;
735
736 /* Now we know for sure we will use the data point */
737 ngood++;
738
739 /* At this point get_data() would interpret the contents of v[]
740 * according to the current plot style and store it.
741 * vfill() has a fixed interpretation of v[] and will use it to
742 * modify the current voxel grid.
743 */
744 modify_voxels( grid, v[0], v[1], v[2], v[3], density_function );
745
746 } /* End of loop over input data points */
747
748 df_close();
749
750 /* We are finished reading user input; return to C locale for internal use */
751 reset_numeric_locale();
752
753 if (ngood == 0) {
754 if (!forever_iteration(plot_iterator))
755 int_warn(NO_CARET,"Skipping data file with no valid points");
756 this_plot->plot_type = NODATA;
757 }
758
759 /* print some basic stats */
760 fprintf(stderr, "\tnumber of points input: %8d\n", ngood);
761 fprintf(stderr, "\tnumber of voxels modified: %8d\n", nvoxels_modified);
762
763 } else if (specs < 0) {
764 this_plot->plot_type = NODATA;
765
766 } else {
767 int_error(NO_CARET, "vfill requires exactly 5 using specs x:y:z:radius:(func)");
768 }
769
770 /* restore original value of sample variables */
771 this_plot->sample_var->udv_value = original_value_sample_var;
772
773 /* Iterate-over-plot mechanism */
774 if (empty_iteration(plot_iterator) && this_plot)
775 this_plot->plot_type = NODATA;
776 if (forever_iteration(plot_iterator) && (this_plot->plot_type == NODATA)) {
777 /* nothing to do here */ ;
778 } else if (next_iteration(plot_iterator)) {
779 c_token = start_token;
780 continue;
781 }
782
783 plot_iterator = cleanup_iteration(plot_iterator);
784 if (equals(c_token, ",")) {
785 c_token++;
786 plot_iterator = check_for_iteration();
787 } else
788 break;
789
790 }
791
792 if (plot_num == 0)
793 int_error(c_token, "no data to plot");
794
795 plot_token = -1;
796
797 }
798
799 /* This is called by vfill for every data point.
800 * It modifies all voxels within radius of the point coordinates.
801 */
802 static void
modify_voxels(t_voxel * grid,double x,double y,double z,double radius,struct at_type * density_function)803 modify_voxels( t_voxel *grid, double x, double y, double z,
804 double radius, struct at_type *density_function )
805 {
806 struct value a;
807 int ix, iy, iz;
808 int ivx, ivy, ivz;
809 int nvx, nvy, nvz;
810 double vx, vy, vz, distance;
811 t_voxel *voxel;
812 int N;
813
814 if (x < current_vgrid->vxmin || x > current_vgrid->vxmax
815 || y < current_vgrid->vymin || y > current_vgrid->vymax
816 || z < current_vgrid->vzmin || z > current_vgrid->vzmax)
817 int_error(NO_CARET, "voxel out of range");
818
819 N = current_vgrid->size;
820
821 /* ivx, ivy, ivz are the indices of the nearest grid point */
822 ivx = ceil( (x - current_vgrid->vxmin) / current_vgrid->vxdelta );
823 ivy = ceil( (y - current_vgrid->vymin) / current_vgrid->vydelta );
824 ivz = ceil( (z - current_vgrid->vzmin) / current_vgrid->vzdelta );
825
826 /* nvx, nvy, nvz are the number of grid points within radius */
827 nvx = floor( radius / current_vgrid->vxdelta );
828 nvy = floor( radius / current_vgrid->vydelta );
829 nvz = floor( radius / current_vgrid->vzdelta );
830
831 /* Only print once */
832 if (nvoxels_modified == 0)
833 fprintf(stderr,
834 "\tradius %g gives a brick of %d voxels on x, %d voxels on y, %d voxels on z\n",
835 radius, 1+2*nvx, 1+2*nvy, 1+2*nvz);
836
837 /* The iteration covers a cube rather than a sphere */
838 evaluate_inside_using = TRUE;
839 for (ix = ivx - nvx; ix <= ivx + nvx; ix++) {
840 if (ix < 0 || ix >= N)
841 continue;
842 for (iy = ivy - nvy; iy <= ivy + nvy; iy++) {
843 if (iy < 0 || iy >= N)
844 continue;
845 for (iz = ivz - nvz; iz <= ivz + nvz; iz++) {
846 int index;
847
848 if (iz < 0 || iz >= N)
849 continue;
850 index = ix + iy * N + iz * N*N;
851 if (index < 0 || index > N*N*N)
852 continue;
853 voxel = ¤t_vgrid->vdata[index];
854
855 /* vx, vy, vz are the true coordinates of this voxel */
856 vx = current_vgrid->vxmin + ix * current_vgrid->vxdelta;
857 vy = current_vgrid->vymin + iy * current_vgrid->vydelta;
858 vz = current_vgrid->vzmin + iz * current_vgrid->vzdelta;
859 distance = sqrt( (vx-x)*(vx-x) + (vy-y)*(vy-y) + (vz-z)*(vz-z) );
860
861 /* Limit the summation to a sphere rather than a cube */
862 /* but always include the voxel nearest the target */
863 if (distance > radius && (ix != ivx || iy != ivy || iz != ivz))
864 continue;
865
866 /* Store in user variable VoxelDistance for use by density_function */
867 udv_VoxelDistance->udv_value.v.cmplx_val.real = distance;
868
869 evaluate_at(density_function, &a);
870 *voxel += real(&a);
871
872 /* Bookkeeping */
873 nvoxels_modified++;
874 }
875 }
876 }
877 evaluate_inside_using = FALSE;
878
879 return;
880 }
881
882 #endif /* VOXEL_GRID_SUPPORT */
883
884 #ifndef VOXEL_GRID_SUPPORT
885 # define NO_SUPPORT int_error(NO_CARET, "this gnuplot does not support voxel grids")
886
check_grid_ranges()887 void check_grid_ranges() { NO_SUPPORT; }
set_vgrid()888 void set_vgrid() { NO_SUPPORT; }
set_vgrid_range()889 void set_vgrid_range() { NO_SUPPORT; }
show_vgrid()890 void show_vgrid() { NO_SUPPORT; }
show_isosurface()891 void show_isosurface() { NO_SUPPORT; }
voxel_command()892 void voxel_command() { NO_SUPPORT; }
vfill_command()893 void vfill_command() { NO_SUPPORT; }
vclear_command()894 void vclear_command() {}
unset_vgrid()895 void unset_vgrid() {}
init_voxelsupport()896 void init_voxelsupport() {}
set_isosurface()897 void set_isosurface() {}
898
get_vgrid_by_name(char * c)899 udvt_entry *get_vgrid_by_name(char *c) { return NULL; }
gpfree_vgrid(struct udvt_entry * x)900 void gpfree_vgrid(struct udvt_entry *x) {}
901
902 #endif /* no VOXEL_GRID_SUPPORT */
903