HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
UT_VoxelArray.C
Go to the documentation of this file.
1 /*
2  * PROPRIETARY INFORMATION. This software is proprietary to
3  * Side Effects Software Inc., and is not to be reproduced,
4  * transmitted, or disclosed in any way without written permission.
5  *
6  * NAME: UT_VoxelArray.C ( UT Library, C++)
7  *
8  * COMMENTS:
9  * Tiled Voxel Array Implementation.
10  */
11 
12 #include "UT_VoxelArray.h"
13 #include "UT_BoundingBox.h"
14 #include "UT_Filter.h"
15 #include "UT_IStream.h"
16 #include "UT_JSONParser.h"
17 #include "UT_JSONWriter.h"
18 #include "UT_SharedMemoryManager.h"
19 #include "UT_StackBuffer.h"
20 #include "UT_VoxelArrayJSON.h"
21 #include <SYS/SYS_Align.h>
22 #include <SYS/SYS_Math.h>
23 
24 ///
25 /// fpreal16 conversion functions
26 ///
27 inline fpreal16
29 inline fpreal16
31 inline fpreal16
33 inline fpreal16
34 UTvoxelConvertFP16(int8 a) { return a; }
35 inline fpreal16
37 inline fpreal16
39 inline fpreal16
41 inline fpreal16
42 UTvoxelConvertFP16(UT_Vector2 a) { return a.x(); }
43 inline fpreal16
44 UTvoxelConvertFP16(UT_Vector3 a) { return a.x(); }
45 inline fpreal16
46 UTvoxelConvertFP16(UT_Vector4 a) { return a.x(); }
47 
48 template <typename T>
50 {
51  return T(v);
52 }
53 
54 template <>
56 {
57  return UT_Vector3F(v, v, v);
58 }
59 template <>
61 {
62  return UT_Vector3D(v, v, v);
63 }
64 template <>
66 {
67  return UT_Vector4F(v, v, v, v);
68 }
69 template <>
71 {
72  return UT_Vector4D(v, v, v, v);
73 }
74 
75 ///
76 /// VoxelTileCompress definitions
77 ///
78 template <typename T>
79 void
81  T &min, T &max) const
82 {
83  int x, y, z;
84 
85  min = getValue(tile, 0, 0, 0);
86  max = min;
87 
88  // Generic approach is to just use the getValue() ability.
89  for (z = 0; z < tile.zres(); z++)
90  {
91  for (y = 0; y < tile.yres(); y++)
92  {
93  for (x = 0; x < tile.xres(); x++)
94  {
95  tile.expandMinMax(getValue(tile, x, y, z), min, max);
96  }
97  }
98  }
99  return;
100 }
101 
102 //
103 // VoxelTile definitions.
104 //
105 
106 template <typename T>
108 {
109  myRes[0] = 0;
110  myRes[1] = 0;
111  myRes[2] = 0;
112 
113  myCompressionType = COMPRESS_CONSTANT;
114  myForeignData = false;
115 
116  if (sizeof(T) <= sizeof(T*))
117  {
118  myData = 0;
119  }
120  else
121  {
122  myData = UT_VOXEL_ALLOC(sizeof(T));
123 
124  // It is not accidental that this is not a static_cast!
125  // There isn't a UT_Vector3(fpreal) constructor (as we have a foolish
126  // UT_Vector3(fpreal *) constructor that would cause massive problems)
127  // but there is a UT_Vector3 operator=(fpreal)!
128  ((T *)myData)[0] = 0;
129  }
130 }
131 
132 template <typename T>
134 {
135  freeData();
136 }
137 
138 template <typename T>
140 {
141  myData = 0;
142 
143  // Use assignment operator.
144  *this = src;
145 }
146 
147 template <typename T>
148 const UT_VoxelTile<T> &
150 {
151  if (&src == this)
152  return *this;
153 
154  freeData();
155 
156  myRes[0] = src.myRes[0];
157  myRes[1] = src.myRes[1];
158  myRes[2] = src.myRes[2];
159 
160  myCompressionType = src.myCompressionType;
161  switch (myCompressionType)
162  {
163  case COMPRESS_RAW:
164  myData = UT_VOXEL_ALLOC(
165  sizeof(T) * myRes[0] * myRes[1] * myRes[2]);
166  memcpy(myData, src.myData, sizeof(T) * myRes[0] * myRes[1] * myRes[2]);
167  break;
168  case COMPRESS_CONSTANT:
169  if (inlineConstant())
170  {
171  myData = src.myData;
172  }
173  else
174  {
175  myData = UT_VOXEL_ALLOC(sizeof(T));
176  memcpy(myData, src.myData, sizeof(T));
177  }
178  break;
179  case COMPRESS_RAWFULL:
180  myData = UT_VOXEL_ALLOC(
181  sizeof(T) * TILESIZE * TILESIZE * myRes[2]);
182  memcpy(myData, src.myData,
183  sizeof(T) * TILESIZE * TILESIZE * myRes[2]);
184  break;
185  case COMPRESS_FPREAL16:
186  myData = UT_VOXEL_ALLOC(
187  sizeof(fpreal16) * myRes[0] * myRes[1] * myRes[2]);
188  memcpy(myData, src.myData,
189  sizeof(fpreal16) * myRes[0] * myRes[1] * myRes[2]);
190  break;
191  default:
192  {
193  UT_VoxelTileCompress<T> *engine;
194 
195  engine = getCompressionEngine(myCompressionType);
196  myData = UT_VOXEL_ALLOC(engine->getDataLength(src));
197  memcpy(myData, src.myData, engine->getDataLength(src));
198  break;
199  }
200  }
201 
202  return *this;
203 }
204 
205 template <typename T>
206 bool
208 {
209  switch (myCompressionType)
210  {
211  case COMPRESS_RAW:
212  // Do the assignment.
213  ((T *)myData)[ ((z * myRes[1]) + y) * myRes[0] + x ] = t;
214  return true;
215 
216  case COMPRESS_CONSTANT:
217  if (rawConstVal() == t)
218  {
219  // Trivially true.
220  return true;
221  }
222  return false;
223 
224  case COMPRESS_RAWFULL:
225  ((T *)myData)[ ((z * TILESIZE) + y) * TILESIZE + x ] = t;
226  return true;
227 
228  case COMPRESS_FPREAL16:
229  return false;
230  }
231 
232  // Use the compression engine.
233  UT_VoxelTileCompress<T> *engine;
234 
235  engine = getCompressionEngine(myCompressionType);
236  return engine->writeThrough(*this, x, y, z, t);
237 }
238 
239 template <typename T>
240 T
241 UT_VoxelTile<T>::lerp(int x, int y, int z, fpreal32 fx, fpreal32 fy, fpreal32 fz) const
242 {
243  T vx, vx1, vy, vy1, vz;
244 
245  switch (myCompressionType)
246  {
247  case COMPRESS_RAW:
248  {
249  T *data = (T *) myData;
250  int offset = (z * myRes[1] + y) * myRes[0] + x;
251  int yinc = myRes[0];
252  int zinc = myRes[0] * myRes[1];
253 
254  // Lerp x:x+1, y, z
255  vx = lerpValues(data[offset], data[offset+1], fx);
256  // Lerp x:x+1, y+1, z
257  vx1 = lerpValues(data[offset+yinc], data[offset+yinc+1], fx);
258 
259  // Lerp x:x+1, y:y+1, z
260  vy = lerpValues(vx, vx1, fy);
261 
262  // Lerp x:x+1, y, z+1
263  vx = lerpValues(data[offset+zinc], data[offset+zinc+1], fx);
264  // Lerp x:x+1, y+1, z+1
265  vx1 = lerpValues(data[offset+zinc+yinc], data[offset+zinc+yinc+1], fx);
266 
267  // Lerp x:x+1, y:y+1, z+1
268  vy1 = lerpValues(vx, vx1, fy);
269 
270  // Lerp x:x+1, y:y+1, z:z+1
271  vz = lerpValues(vy, vy1, fz);
272  break;
273  }
274  case COMPRESS_RAWFULL:
275  {
276  T *data = (T *) myData;
277  int offset = (z * TILESIZE + y) * TILESIZE + x;
278  int yinc = TILESIZE;
279  int zinc = TILESIZE * TILESIZE;
280 
281  // Lerp x:x+1, y, z
282  vx = lerpValues(data[offset], data[offset+1], fx);
283  // Lerp x:x+1, y+1, z
284  vx1 = lerpValues(data[offset+yinc], data[offset+yinc+1], fx);
285 
286  // Lerp x:x+1, y:y+1, z
287  vy = lerpValues(vx, vx1, fy);
288 
289  // Lerp x:x+1, y, z+1
290  vx = lerpValues(data[offset+zinc], data[offset+zinc+1], fx);
291  // Lerp x:x+1, y+1, z+1
292  vx1 = lerpValues(data[offset+zinc+yinc], data[offset+zinc+yinc+1], fx);
293 
294  // Lerp x:x+1, y:y+1, z+1
295  vy1 = lerpValues(vx, vx1, fy);
296 
297  // Lerp x:x+1, y:y+1, z:z+1
298  vz = lerpValues(vy, vy1, fz);
299  break;
300  }
301  case COMPRESS_FPREAL16:
302  {
303  fpreal16 *data = (fpreal16 *) myData;
304  int offset = (z * myRes[1] + y) * myRes[0] + x;
305  int yinc = myRes[0];
306  int zinc = myRes[0] * myRes[1];
307  fpreal16 vx, vx1, vy, vy1, vz;
308 
309  // Lerp x:x+1, y, z
310  vx = SYSlerp(data[offset], data[offset+1], fx);
311  // Lerp x:x+1, y+1, z
312  vx1 = SYSlerp(data[offset+yinc], data[offset+yinc+1], fx);
313 
314  // Lerp x:x+1, y:y+1, z
315  vy = SYSlerp(vx, vx1, fy);
316 
317  // Lerp x:x+1, y, z+1
318  vx = SYSlerp(data[offset+zinc], data[offset+zinc+1], fx);
319  // Lerp x:x+1, y+1, z+1
320  vx1 = SYSlerp(data[offset+zinc+yinc], data[offset+zinc+yinc+1], fx);
321 
322  // Lerp x:x+1, y:y+1, z+1
323  vy1 = SYSlerp(vx, vx1, fy);
324 
325  // Lerp x:x+1, y:y+1, z:z+1
326  vz = SYSlerp(vy, vy1, fz);
327 
328  return UTvoxelConvertFromFP16<T>(vz);
329  }
330  case COMPRESS_CONSTANT:
331  {
332  // This is quite trivial to do a trilerp on.
333  vz = rawConstVal();
334  break;
335  }
336 
337  default:
338  {
339  UT_VoxelTileCompress<T> *engine;
340 
341  engine = getCompressionEngine(myCompressionType);
342  // Lerp x:x+1, y, z
343  vx = lerpValues(engine->getValue(*this, x, y, z),
344  engine->getValue(*this, x+1, y, z),
345  fx);
346  // Lerp x:x+1, y+1, z
347  vx1 = lerpValues(engine->getValue(*this, x, y+1, z),
348  engine->getValue(*this, x+1, y+1, z),
349  fx);
350 
351  // Lerp x:x+1, y:y+1, z
352  vy = lerpValues(vx, vx1, fy);
353 
354  // Lerp x:x+1, y, z+1
355  vx = lerpValues(engine->getValue(*this, x, y, z+1),
356  engine->getValue(*this, x+1, y, z+1),
357  fx);
358  // Lerp x:x+1, y+1, z+1
359  vx1 = lerpValues(engine->getValue(*this, x, y+1, z+1),
360  engine->getValue(*this, x+1, y+1, z+1),
361  fx);
362 
363  // Lerp x:x+1, y:y+1, z+1
364  vy1 = lerpValues(vx, vx1, fy);
365 
366  // Lerp x:x+1, y:y+1, z:z+1
367  vz = lerpValues(vy, vy1, fz);
368  break;
369  }
370  }
371 
372  return vz;
373 }
374 
375 template <typename T>
376 template <int AXIS2D>
377 T
378 UT_VoxelTile<T>::lerpAxis(int x, int y, int z, fpreal32 fx, fpreal32 fy, fpreal32 fz) const
379 {
380  T vx, vx1, vy, vy1, vz;
381 
382  switch (myCompressionType)
383  {
384  case COMPRESS_RAW:
385  {
386  T *data = (T *) myData;
387  int offset = (z * myRes[1] + y) * myRes[0] + x;
388  int yinc = myRes[0];
389  int zinc = myRes[0] * myRes[1];
390 
391  // Lerp x:x+1, y, z
392  if (AXIS2D != 0)
393  vx = lerpValues(data[offset],
394  data[offset+1],
395  fx);
396  else
397  vx = data[offset];
398 
399  if (AXIS2D != 1)
400  {
401  // Lerp x:x+1, y+1, z
402  if (AXIS2D != 0)
403  vx1= lerpValues(data[offset+yinc],
404  data[offset+yinc+1],
405  fx);
406  else
407  vx1 = data[offset+yinc];
408  // Lerp x:x+1, y:y+1, z
409  vy = lerpValues(vx, vx1, fy);
410  }
411  else
412  vy = vx;
413 
414  if (AXIS2D != 2)
415  {
416  // Lerp x:x+1, y, z+1
417  if (AXIS2D != 0)
418  vx = lerpValues(data[offset+zinc],
419  data[offset+zinc+1],
420  fx);
421  else
422  vx = data[offset+zinc];
423 
424  if (AXIS2D != 1)
425  {
426  // Lerp x:x+1, y+1, z+1
427  if (AXIS2D != 0)
428  vx1= lerpValues(data[offset+zinc+yinc],
429  data[offset+zinc+yinc+1],
430  fx);
431  else
432  vx1 = data[offset+zinc+yinc];
433  // Lerp x:x+1, y:y+1, z+1
434  vy1 = lerpValues(vx, vx1, fy);
435  }
436  else
437  vy1 = vx;
438 
439  // Lerp x:x+1, y:y+1, z:z+1
440  vz = lerpValues(vy, vy1, fz);
441  }
442  else
443  vz = vy;
444  break;
445  }
446  case COMPRESS_RAWFULL:
447  {
448  T *data = (T *) myData;
449  int offset = (z * TILESIZE + y) * TILESIZE + x;
450  int yinc = TILESIZE;
451  int zinc = TILESIZE * TILESIZE;
452 
453  // Lerp x:x+1, y, z
454  if (AXIS2D != 0)
455  vx = lerpValues(data[offset],
456  data[offset+1],
457  fx);
458  else
459  vx = data[offset];
460 
461  if (AXIS2D != 1)
462  {
463  // Lerp x:x+1, y+1, z
464  if (AXIS2D != 0)
465  vx1= lerpValues(data[offset+yinc],
466  data[offset+yinc+1],
467  fx);
468  else
469  vx1 = data[offset+yinc];
470  // Lerp x:x+1, y:y+1, z
471  vy = lerpValues(vx, vx1, fy);
472  }
473  else
474  vy = vx;
475 
476  if (AXIS2D != 2)
477  {
478  // Lerp x:x+1, y, z+1
479  if (AXIS2D != 0)
480  vx = lerpValues(data[offset+zinc],
481  data[offset+zinc+1],
482  fx);
483  else
484  vx = data[offset+zinc];
485 
486  if (AXIS2D != 1)
487  {
488  // Lerp x:x+1, y+1, z+1
489  if (AXIS2D != 0)
490  vx1= lerpValues(data[offset+zinc+yinc],
491  data[offset+zinc+yinc+1],
492  fx);
493  else
494  vx1 = data[offset+zinc+yinc];
495  // Lerp x:x+1, y:y+1, z+1
496  vy1 = lerpValues(vx, vx1, fy);
497  }
498  else
499  vy1 = vx;
500 
501  // Lerp x:x+1, y:y+1, z:z+1
502  vz = lerpValues(vy, vy1, fz);
503  }
504  else
505  vz = vy;
506  break;
507  }
508  case COMPRESS_FPREAL16:
509  {
510  fpreal16 *data = (fpreal16 *) myData;
511  int offset = (z * myRes[1] + y) * myRes[0] + x;
512  int yinc = myRes[0];
513  int zinc = myRes[0] * myRes[1];
514  fpreal16 vx, vx1, vy, vy1, vz;
515 
516  // Lerp x:x+1, y, z
517  if (AXIS2D != 0)
518  vx = SYSlerp(data[offset],
519  data[offset+1],
520  fx);
521  else
522  vx = data[offset];
523 
524  if (AXIS2D != 1)
525  {
526  // Lerp x:x+1, y+1, z
527  if (AXIS2D != 0)
528  vx1= SYSlerp(data[offset+yinc],
529  data[offset+yinc+1],
530  fx);
531  else
532  vx1 = data[offset+yinc];
533  // Lerp x:x+1, y:y+1, z
534  vy = SYSlerp(vx, vx1, fy);
535  }
536  else
537  vy = vx;
538 
539  if (AXIS2D != 2)
540  {
541  // Lerp x:x+1, y, z+1
542  if (AXIS2D != 0)
543  vx = SYSlerp(data[offset+zinc],
544  data[offset+zinc+1],
545  fx);
546  else
547  vx = data[offset+zinc];
548 
549  if (AXIS2D != 1)
550  {
551  // Lerp x:x+1, y+1, z+1
552  if (AXIS2D != 0)
553  vx1= SYSlerp(data[offset+zinc+yinc],
554  data[offset+zinc+yinc+1],
555  fx);
556  else
557  vx1 = data[offset+zinc+yinc];
558  // Lerp x:x+1, y:y+1, z+1
559  vy1 = SYSlerp(vx, vx1, fy);
560  }
561  else
562  vy1 = vx;
563 
564  // Lerp x:x+1, y:y+1, z:z+1
565  vz = SYSlerp(vy, vy1, fz);
566  }
567  else
568  vz = vy;
569 
570  return UTvoxelConvertFromFP16<T>(vz);
571  }
572  case COMPRESS_CONSTANT:
573  {
574  // This is quite trivial to do a bilerp on.
575  vz = rawConstVal();
576  break;
577  }
578 
579  default:
580  {
581  UT_VoxelTileCompress<T> *engine;
582 
583  engine = getCompressionEngine(myCompressionType);
584  // Lerp x:x+1, y, z
585  if (AXIS2D != 0)
586  vx = lerpValues(engine->getValue(*this, x, y, z),
587  engine->getValue(*this, x+1, y, z),
588  fx);
589  else
590  vx = engine->getValue(*this, x, y, z);
591 
592  if (AXIS2D != 1)
593  {
594  // Lerp x:x+1, y+1, z
595  if (AXIS2D != 0)
596  vx1= lerpValues(engine->getValue(*this, x, y+1, z),
597  engine->getValue(*this, x+1, y+1, z),
598  fx);
599  else
600  vx1 = engine->getValue(*this, x, y+1, z);
601  // Lerp x:x+1, y:y+1, z
602  vy = lerpValues(vx, vx1, fy);
603  }
604  else
605  vy = vx;
606 
607  if (AXIS2D != 2)
608  {
609  // Lerp x:x+1, y, z+1
610  if (AXIS2D != 0)
611  vx = lerpValues(engine->getValue(*this, x, y, z+1),
612  engine->getValue(*this, x+1, y, z+1),
613  fx);
614  else
615  vx = engine->getValue(*this, x, y, z+1);
616 
617  if (AXIS2D != 1)
618  {
619  // Lerp x:x+1, y+1, z+1
620  if (AXIS2D != 0)
621  vx1= lerpValues(engine->getValue(*this, x, y+1, z+1),
622  engine->getValue(*this, x+1, y+1, z+1),
623  fx);
624  else
625  vx1 = engine->getValue(*this, x, y+1, z+1);
626  // Lerp x:x+1, y:y+1, z+1
627  vy1 = lerpValues(vx, vx1, fy);
628  }
629  else
630  vy1 = vx;
631 
632  // Lerp x:x+1, y:y+1, z:z+1
633  vz = lerpValues(vy, vy1, fz);
634  }
635  else
636  vz = vy;
637  break;
638  }
639  }
640 
641  return vz;
642 }
643 
644 template <typename T>
645 bool
646 UT_VoxelTile<T>::extractSample(int x, int y, int z, T *sample) const
647 {
648  switch (myCompressionType)
649  {
650  case COMPRESS_RAW:
651  {
652  T *data = (T *) myData;
653  int offset = (z * myRes[1] + y) * myRes[0] + x;
654  int yinc = myRes[0];
655  int zinc = myRes[0] * myRes[1];
656 
657  sample[0] = data[offset];
658  sample[1] = data[offset+1];
659  sample[2+0] = data[offset+yinc];
660  sample[2+1] = data[offset+yinc+1];
661  sample[4+0] = data[offset+zinc];
662  sample[4+1] = data[offset+zinc+1];
663  sample[4+2+0] = data[offset+zinc+yinc];
664  sample[4+2+1] = data[offset+zinc+yinc+1];
665  break;
666  }
667  case COMPRESS_RAWFULL:
668  {
669  T *data = (T *) myData;
670  int offset = (z * TILESIZE + y) * TILESIZE + x;
671  int yinc = TILESIZE;
672  int zinc = TILESIZE * TILESIZE;
673 
674  sample[0] = data[offset];
675  sample[1] = data[offset+1];
676  sample[2+0] = data[offset+yinc];
677  sample[2+1] = data[offset+yinc+1];
678  sample[4+0] = data[offset+zinc];
679  sample[4+1] = data[offset+zinc+1];
680  sample[4+2+0] = data[offset+zinc+yinc];
681  sample[4+2+1] = data[offset+zinc+yinc+1];
682  break;
683  }
684  case COMPRESS_FPREAL16:
685  {
686  fpreal16 *data = (fpreal16 *) myData;
687  int offset = (z * myRes[1] + y) * myRes[0] + x;
688  int yinc = myRes[0];
689  int zinc = myRes[0] * myRes[1];
690 
691  sample[0] = data[offset];
692  sample[1] = data[offset+1];
693  sample[2+0] = data[offset+yinc];
694  sample[2+1] = data[offset+yinc+1];
695  sample[4+0] = data[offset+zinc];
696  sample[4+1] = data[offset+zinc+1];
697  sample[4+2+0] = data[offset+zinc+yinc];
698  sample[4+2+1] = data[offset+zinc+yinc+1];
699  break;
700  }
701  case COMPRESS_CONSTANT:
702  {
703  sample[0] = rawConstVal();
704  return true;
705  }
706 
707  default:
708  {
709  UT_VoxelTileCompress<T> *engine;
710 
711  engine = getCompressionEngine(myCompressionType);
712  // Lerp x:x+1, y, z
713  sample[0] = engine->getValue(*this, x, y, z);
714  sample[1] = engine->getValue(*this, x+1, y, z);
715  sample[2+0] = engine->getValue(*this, x, y+1, z);
716  sample[2+1] = engine->getValue(*this, x+1, y+1, z);
717  sample[4+0] = engine->getValue(*this, x, y, z+1);
718  sample[4+1] = engine->getValue(*this, x+1, y, z+1);
719  sample[4+2+0] = engine->getValue(*this, x, y+1, z+1);
720  sample[4+2+1] = engine->getValue(*this, x+1, y+1, z+1);
721  break;
722  }
723  }
724  return false;
725 }
726 
727 template <typename T>
728 bool
729 UT_VoxelTile<T>::extractSamplePlus(int x, int y, int z, T *sample) const
730 {
731  switch (myCompressionType)
732  {
733  case COMPRESS_RAW:
734  {
735  T *data = (T *) myData;
736  int offset = (z * myRes[1] + y) * myRes[0] + x;
737  int yinc = myRes[0];
738  int zinc = myRes[0] * myRes[1];
739 
740  sample[0] = data[offset-1];
741  sample[1] = data[offset+1];
742  sample[2+0] = data[offset-yinc];
743  sample[2+1] = data[offset+yinc];
744  sample[4+0] = data[offset-zinc];
745  sample[4+1] = data[offset+zinc];
746  sample[6] = data[offset];
747  break;
748  }
749  case COMPRESS_RAWFULL:
750  {
751  T *data = (T *) myData;
752  int offset = (z * TILESIZE + y) * TILESIZE + x;
753  int yinc = TILESIZE;
754  int zinc = TILESIZE * TILESIZE;
755 
756  sample[0] = data[offset-1];
757  sample[1] = data[offset+1];
758  sample[2+0] = data[offset-yinc];
759  sample[2+1] = data[offset+yinc];
760  sample[4+0] = data[offset-zinc];
761  sample[4+1] = data[offset+zinc];
762  sample[6] = data[offset];
763  break;
764  }
765  case COMPRESS_FPREAL16:
766  {
767  fpreal16 *data = (fpreal16 *) myData;
768  int offset = (z * myRes[1] + y) * myRes[0] + x;
769  int yinc = myRes[0];
770  int zinc = myRes[0] * myRes[1];
771 
772  sample[0] = data[offset-1];
773  sample[1] = data[offset+1];
774  sample[2+0] = data[offset-yinc];
775  sample[2+1] = data[offset+yinc];
776  sample[4+0] = data[offset-zinc];
777  sample[4+1] = data[offset+zinc];
778  sample[6] = data[offset];
779  break;
780  }
781  case COMPRESS_CONSTANT:
782  {
783  sample[0] = rawConstVal();
784  return true;
785  }
786 
787  default:
788  {
789  UT_VoxelTileCompress<T> *engine;
790 
791  engine = getCompressionEngine(myCompressionType);
792  // Lerp x:x+1, y, z
793  sample[0] = engine->getValue(*this, x-1, y, z);
794  sample[1] = engine->getValue(*this, x+1, y, z);
795  sample[2+0] = engine->getValue(*this, x, y-1, z);
796  sample[2+1] = engine->getValue(*this, x, y+1, z);
797  sample[4+0] = engine->getValue(*this, x, y, z+1);
798  sample[4+1] = engine->getValue(*this, x, y, z-1);
799  sample[6] = engine->getValue(*this, x, y, z);
800  break;
801  }
802  }
803  return false;
804 }
805 
806 template <typename T>
807 bool
808 UT_VoxelTile<T>::extractSampleCube(int x, int y, int z, T *sample) const
809 {
810  switch (myCompressionType)
811  {
812  case COMPRESS_RAW:
813  {
814  T *data = (T *) myData;
815  int offset = (z * myRes[1] + y) * myRes[0] + x;
816  int yinc = myRes[0];
817  int zinc = myRes[0] * myRes[1];
818  int sampidx = 0;
819 
820  offset -= zinc;
821  for (int dz = -1; dz <= 1; dz++)
822  {
823  offset -= yinc;
824  for (int dy = -1; dy <= 1; dy++)
825  {
826  sample[sampidx] = data[offset-1];
827  sample[sampidx+1] = data[offset];
828  sample[sampidx+2] = data[offset+1];
829  sampidx += 3;
830  offset += yinc;
831  }
832  offset -= yinc * 3;
833  offset += zinc;
834  }
835  break;
836  }
837  case COMPRESS_RAWFULL:
838  {
839  T *data = (T *) myData;
840  int offset = (z * TILESIZE + y) * TILESIZE + x;
841  int yinc = TILESIZE;
842  int zinc = TILESIZE * TILESIZE;
843  int sampidx = 0;
844 
845  offset -= zinc;
846  for (int dz = -1; dz <= 1; dz++)
847  {
848  offset -= yinc;
849  for (int dy = -1; dy <= 1; dy++)
850  {
851  sample[sampidx] = data[offset-1];
852  sample[sampidx+1] = data[offset];
853  sample[sampidx+2] = data[offset+1];
854  sampidx += 3;
855  offset += yinc;
856  }
857  offset -= yinc * 3;
858  offset += zinc;
859  }
860  break;
861  }
862  case COMPRESS_FPREAL16:
863  {
864  fpreal16 *data = (fpreal16 *) myData;
865  int offset = (z * myRes[1] + y) * myRes[0] + x;
866  int yinc = myRes[0];
867  int zinc = myRes[0] * myRes[1];
868  int sampidx = 0;
869 
870  offset -= zinc;
871  for (int dz = -1; dz <= 1; dz++)
872  {
873  offset -= yinc;
874  for (int dy = -1; dy <= 1; dy++)
875  {
876  sample[sampidx] = data[offset-1];
877  sample[sampidx+1] = data[offset];
878  sample[sampidx+2] = data[offset+1];
879  sampidx += 3;
880  offset += yinc;
881  }
882  offset -= yinc * 3;
883  offset += zinc;
884  }
885  break;
886  }
887  case COMPRESS_CONSTANT:
888  {
889  sample[0] = rawConstVal();
890  return true;
891  }
892 
893  default:
894  {
895  UT_VoxelTileCompress<T> *engine;
896 
897  engine = getCompressionEngine(myCompressionType);
898  int sampidx = 0;
899  // Lerp x:x+1, y, z
900  for (int dz = -1; dz <= 1; dz++)
901  {
902  for (int dy = -1; dy <= 1; dy++)
903  {
904  for (int dx = -1; dx <= 1; dx++)
905  {
906  sample[sampidx++] = engine->getValue(*this, x+dx, y+dy, z+dz);
907  }
908  }
909  }
910  break;
911  }
912  }
913  return false;
914 }
915 
916 template <typename T>
917 template <int AXIS2D>
918 bool
919 UT_VoxelTile<T>::extractSampleAxis(int x, int y, int z, T *sample) const
920 {
921  switch (myCompressionType)
922  {
923  case COMPRESS_RAW:
924  {
925  T *data = (T *) myData;
926  int offset = (z * myRes[1] + y) * myRes[0] + x;
927  int yinc = myRes[0];
928  int zinc = myRes[0] * myRes[1];
929 
930  sample[0] = data[offset];
931  if (AXIS2D != 0)
932  sample[1] = data[offset+1];
933  if (AXIS2D != 1)
934  {
935  sample[2+0] = data[offset+yinc];
936  if (AXIS2D != 0)
937  sample[2+1] = data[offset+yinc+1];
938  }
939  if (AXIS2D != 2)
940  {
941  sample[4+0] = data[offset+zinc];
942  if (AXIS2D != 0)
943  sample[4+1] = data[offset+zinc+1];
944  if (AXIS2D != 1)
945  {
946  sample[4+2+0] = data[offset+zinc+yinc];
947  if (AXIS2D != 0)
948  sample[4+2+1] = data[offset+zinc+yinc+1];
949  }
950  }
951  break;
952  }
953  case COMPRESS_RAWFULL:
954  {
955  T *data = (T *) myData;
956  int offset = (z * TILESIZE + y) * TILESIZE + x;
957  int yinc = TILESIZE;
958  int zinc = TILESIZE * TILESIZE;
959 
960  sample[0] = data[offset];
961  if (AXIS2D != 0)
962  sample[1] = data[offset+1];
963  if (AXIS2D != 1)
964  {
965  sample[2+0] = data[offset+yinc];
966  if (AXIS2D != 0)
967  sample[2+1] = data[offset+yinc+1];
968  }
969  if (AXIS2D != 2)
970  {
971  sample[4+0] = data[offset+zinc];
972  if (AXIS2D != 0)
973  sample[4+1] = data[offset+zinc+1];
974  if (AXIS2D != 1)
975  {
976  sample[4+2+0] = data[offset+zinc+yinc];
977  if (AXIS2D != 0)
978  sample[4+2+1] = data[offset+zinc+yinc+1];
979  }
980  }
981  break;
982  }
983  case COMPRESS_FPREAL16:
984  {
985  fpreal16 *data = (fpreal16 *) myData;
986  int offset = (z * myRes[1] + y) * myRes[0] + x;
987  int yinc = myRes[0];
988  int zinc = myRes[0] * myRes[1];
989 
990  sample[0] = data[offset];
991  if (AXIS2D != 0)
992  sample[1] = data[offset+1];
993  if (AXIS2D != 1)
994  {
995  sample[2+0] = data[offset+yinc];
996  if (AXIS2D != 0)
997  sample[2+1] = data[offset+yinc+1];
998  }
999  if (AXIS2D != 2)
1000  {
1001  sample[4+0] = data[offset+zinc];
1002  if (AXIS2D != 0)
1003  sample[4+1] = data[offset+zinc+1];
1004  if (AXIS2D != 1)
1005  {
1006  sample[4+2+0] = data[offset+zinc+yinc];
1007  if (AXIS2D != 0)
1008  sample[4+2+1] = data[offset+zinc+yinc+1];
1009  }
1010  }
1011  break;
1012  }
1013  case COMPRESS_CONSTANT:
1014  {
1015  sample[0] = rawConstVal();
1016  return true;
1017  }
1018 
1019  default:
1020  {
1021  UT_VoxelTileCompress<T> *engine;
1022 
1023  engine = getCompressionEngine(myCompressionType);
1024  // Lerp x:x+1, y, z
1025  sample[0] = engine->getValue(*this, x, y, z);
1026  if (AXIS2D != 0)
1027  sample[1] = engine->getValue(*this, x+1, y, z);
1028  if (AXIS2D != 1)
1029  {
1030  sample[2+0] = engine->getValue(*this, x, y+1, z);
1031  if (AXIS2D != 0)
1032  sample[2+1] = engine->getValue(*this, x+1, y+1, z);
1033  }
1034  if (AXIS2D != 2)
1035  {
1036  sample[4+0] = engine->getValue(*this, x, y, z+1);
1037  if (AXIS2D != 0)
1038  sample[4+1] = engine->getValue(*this, x+1, y, z+1);
1039  if (AXIS2D != 1)
1040  {
1041  sample[4+2+0] = engine->getValue(*this, x, y+1, z+1);
1042  if (AXIS2D != 0)
1043  sample[4+2+1] = engine->getValue(*this, x+1, y+1, z+1);
1044  }
1045  }
1046  break;
1047  }
1048  }
1049  return false;
1050 }
1051 
1052 #if 0
1053 template <typename T>
1054 T
1055 UT_VoxelTile<T>::lerp(v4uf frac, int x, int y, int z) const
1056 {
1057  v4uf a, b;
1058 
1059  switch (myCompressionType)
1060  {
1061  case COMPRESS_RAW:
1062  case COMPRESS_RAWFULL:
1063  {
1064  T *data = (T *) myData;
1065  int offset = (z * myRes[1] + y) * myRes[0] + x;
1066  int yinc = myRes[0];
1067  int zinc = myRes[0] * myRes[1];
1068 
1069  a = v4uf( data[offset],
1070  data[offset+zinc],
1071  data[offset+yinc],
1072  data[offset+yinc+zinc] );
1073  b = v4uf( data[offset+1],
1074  data[offset+zinc+1],
1075  data[offset+yinc+1],
1076  data[offset+yinc+zinc+1] );
1077  break;
1078  }
1079 
1080  case COMPRESS_CONSTANT:
1081  {
1082  // This is quite trivial to do a trilerp on.
1083  return rawConstVal();
1084  }
1085 
1086  default:
1087  {
1088  UT_VoxelTileCompress<T> *engine;
1089 
1090  engine = getCompressionEngine(myCompressionType);
1091  // Lerp x:x+1, y, z
1092  a = v4uf( engine->getValue(*this, x, y, z),
1093  engine->getValue(*this, x, y, z+1),
1094  engine->getValue(*this, x, y+1, z),
1095  engine->getValue(*this, x, y+1, z+1) );
1096  b = v4uf( engine->getValue(*this, x+1, y, z),
1097  engine->getValue(*this, x+1, y, z+1),
1098  engine->getValue(*this, x+1, y+1, z),
1099  engine->getValue(*this, x+1, y+1, z+1) );
1100  break;
1101  }
1102  }
1103 
1104  v4uf fx, fy, fz;
1105 
1106  fx = frac.swizzle<0, 0, 0, 0>();
1107  fy = frac.swizzle<1, 1, 1, 1>();
1108  fz = frac.swizzle<2, 2, 2, 2>();
1109 
1110  b -= a;
1111  a = madd(b, fx, a);
1112 
1113  b = a.swizzle<2, 3, 0, 1>();
1114  b -= a;
1115  a = madd(b, fy, a);
1116 
1117  b = a.swizzle<1, 2, 3, 0>();
1118  b -= a;
1119  a = madd(b, fz, a);
1120 
1121  return a[0];
1122 }
1123 #endif
1124 
1125 template <typename T>
1126 T *
1127 UT_VoxelTile<T>::fillCacheLine(T *cacheline, int &stride, int x, int y, int z, bool forcecopy, bool strideofone) const
1128 {
1129  UT_ASSERT_P(x >= 0 && y >= 0 && z >= 0);
1130  UT_ASSERT_P(x < myRes[0] && y < myRes[1] && z < myRes[2]);
1131 
1132  T *src;
1133  int i, xres = myRes[0];
1134 
1135  // All the directly handled types exit directly from this switch.
1136  switch (myCompressionType)
1137  {
1138  case COMPRESS_RAW:
1139  stride = 1;
1140  src = (T *)myData + (z * myRes[1] + y) * xres;
1141  if (!forcecopy)
1142  return &src[x];
1143 
1144  for (i = 0; i < xres; i++)
1145  cacheline[i] = src[i];
1146 
1147  return &cacheline[x];
1148 
1149  case COMPRESS_FPREAL16:
1150  {
1151  fpreal16 *src;
1152 
1153  stride = 1;
1154  src = (fpreal16 *)myData + (z * myRes[1] + y) * xres;
1155 
1156  for (i = 0; i < xres; i++)
1157  cacheline[i] = src[i];
1158 
1159  return &cacheline[x];
1160  }
1161 
1162 
1163  case COMPRESS_CONSTANT:
1164  src = rawConstData();
1165  if (!forcecopy && !strideofone)
1166  {
1167  stride = 0;
1168  return src;
1169  }
1170  stride = 1;
1171  // Fill out a constant value
1172  for (i = 0; i < xres; i++)
1173  cacheline[i] = *src;
1174 
1175  return &cacheline[x];
1176 
1177 
1178  case COMPRESS_RAWFULL:
1179  stride = 1;
1180  src = (T *)myData + (z * TILESIZE + y) * TILESIZE;
1181  if (!forcecopy)
1182  return &src[x];
1183 
1184  for (i = 0; i < xres; i++)
1185  cacheline[i] = src[i];
1186 
1187  return &cacheline[x];
1188  }
1189 
1190  // By default use the compression engine.
1191  UT_VoxelTileCompress<T> *engine;
1192 
1193  engine = getCompressionEngine(myCompressionType);
1194 
1195  // We could add support for a direct cacheline fill to accelerate
1196  // this case as well.
1197  stride = 1;
1198  for (i = 0; i < xres; i++)
1199  cacheline[i] = engine->getValue(*this, i, y, z);
1200 
1201  return &cacheline[x];
1202 }
1203 
1204 template <typename T>
1205 void
1206 UT_VoxelTile<T>::writeCacheLine(T *cacheline, int y, int z)
1207 {
1208  UT_ASSERT_P(y >= 0 && z >= 0);
1209  UT_ASSERT_P(y < myRes[1] && z < myRes[2]);
1210 
1211  T *dst, value;
1212  int i, xres = myRes[0];
1213 
1214  // All the directly handled types exit directly from this switch.
1215  switch (myCompressionType)
1216  {
1217  case COMPRESS_RAW:
1218  dst = (T *)myData + (z * myRes[1] + y) * xres;
1219  for (i = 0; i < xres; i++)
1220  *dst++ = *cacheline++;
1221  return;
1222 
1223  case COMPRESS_CONSTANT:
1224  value = rawConstVal();
1225  for (i = 0; i < xres; i++)
1226  if (cacheline[i] != value)
1227  break;
1228  // If everything was equal, our write is trivial.
1229  if (i == xres)
1230  return;
1231 
1232  break;
1233 
1234  case COMPRESS_RAWFULL:
1235  dst = (T *)myData + (z * TILESIZE + y) * TILESIZE;
1236  for (i = 0; i < TILESIZE; i++)
1237  *dst++ = *cacheline++;
1238 
1239  return;
1240  }
1241 
1242  // Switch back to invoking writeThrough. Ideally we'd have
1243  // a version that can handle a whole cache line at once
1244  for (i = 0; i < xres; i++)
1245  if (!writeThrough(i, y, z, cacheline[i]))
1246  break;
1247 
1248  // Determine if we failed to write everything through
1249  if (i != xres)
1250  {
1251  // Uncompress and reinvoke ourselves.
1252  uncompress();
1253  writeCacheLine(cacheline, y, z);
1254  }
1255 }
1256 
1257 template <typename T>
1258 void
1259 UT_VoxelTile<T>::copyFragment(int dstx, int dsty, int dstz,
1260  const UT_VoxelTile<T> &srctile,
1261  int srcx, int srcy, int srcz)
1262 {
1263  int w = SYSmin(xres() - dstx, srctile.xres() - srcx);
1264  int h = SYSmin(yres() - dsty, srctile.yres() - srcy);
1265  int d = SYSmin(zres() - dstz, srctile.zres() - srcz);
1266 
1267 #if 1
1268  if (srctile.isSimpleCompression())
1269  {
1270  T *dst;
1271  const T *src;
1272  int srcinc;
1273 
1274  src = srctile.rawData();
1275  srcinc = srctile.isConstant() ? 0 : 1;
1276 
1277  // Ensure we are easy to write to.
1278  uncompress();
1279 
1280  dst = rawData();
1281  dst += dstx + (dsty + dstz*yres())*xres();
1282 
1283  if (srcinc)
1284  src += srcx + (srcy + srcz*srctile.yres())*srctile.xres();
1285 
1286  for (int z = 0; z < d; z++)
1287  {
1288  for (int y = 0; y < h; y++)
1289  {
1290  if (srcinc)
1291  {
1292  for (int x = 0; x < w; x++)
1293  dst[x] = src[x];
1294  }
1295  else
1296  {
1297  for (int x = 0; x < w; x++)
1298  dst[x] = *src;
1299  }
1300  dst += xres();
1301  if (srcinc)
1302  src += srctile.xres();
1303  }
1304  dst += (yres()-h) * xres();
1305  if (srcinc)
1306  src += (srctile.yres() - h) * srctile.xres();
1307  }
1308 
1309  return;
1310  }
1311 #endif
1312 
1313  // Fail safe approach.
1314  for (int z = 0; z < d; z++)
1315  for (int y = 0; y < h; y++)
1316  for (int x = 0; x < w; x++)
1317  {
1318  setValue(dstx+x, dsty+y, dstz+z,
1319  srctile(srcx+x, srcy+y, srcz+z));
1320  }
1321 }
1322 
1323 template <typename T>
1324 template <typename S>
1325 void
1327 {
1328  int w = xres();
1329  int h = yres();
1330  int d = zres();
1331 
1332  if (isSimpleCompression())
1333  {
1334  const T *src;
1335  int srcinc;
1336 
1337  src = rawData();
1338  srcinc = isConstant() ? 0 : 1;
1339 
1340  if (stride == 1)
1341  {
1342  if (srcinc == 1)
1343  {
1344  // Super trivial!
1345  for (int i = 0; i < w * h * d; i++)
1346  {
1347  *dst++ = T(*src++);
1348  }
1349  }
1350  else
1351  {
1352  // Constant, also trivial!
1353  T cval = T(*src);
1354 
1355  for (int i = 0; i < w * h * d; i++)
1356  {
1357  *dst++ = cval;
1358  }
1359  }
1360  }
1361  else
1362  {
1363  for (int z = 0; z < d; z++)
1364  {
1365  for (int y = 0; y < h; y++)
1366  {
1367  if (srcinc)
1368  {
1369  for (int x = 0; x < w; x++)
1370  {
1371  *dst = S(src[x]);
1372  dst += stride;
1373  }
1374  }
1375  else
1376  {
1377  for (int x = 0; x < w; x++)
1378  {
1379  *dst = S(*src);
1380  dst += stride;
1381  }
1382  }
1383  if (srcinc)
1384  src += w;
1385  }
1386  }
1387  }
1388 
1389  return;
1390  }
1391 
1392  // Fail safe approach.
1393  for (int z = 0; z < d; z++)
1394  for (int y = 0; y < h; y++)
1395  for (int x = 0; x < w; x++)
1396  {
1397  *dst = S((*this)(x, y, z));
1398  dst += stride;
1399  }
1400 }
1401 
1402 template <typename T>
1403 template <typename S>
1404 void
1405 UT_VoxelTile<T>::writeData(const S *srcdata, int srcstride)
1406 {
1407  int w = xres();
1408  int h = yres();
1409  int d = zres();
1410  int i, n = w * h * d;
1411 
1412  // Check if src is constant
1413  S compare = srcdata[0];
1414  int srcoff = srcstride;
1415 
1416  if (srcstride == 0)
1417  {
1418  // Trivially constant
1419  makeConstant(T(compare));
1420  return;
1421  }
1422 
1423  for (i = 1; i < n; i++)
1424  {
1425  if (srcdata[srcoff] != compare)
1426  break;
1427  srcoff += srcstride;
1428  }
1429 
1430  if (i == n)
1431  {
1432  // Constant source!
1433  makeConstant(compare);
1434  return;
1435  }
1436 
1437  // Create a raw tile, expanding constants
1438  uncompress();
1439 
1440  if (srcstride == 1)
1441  {
1442  T *dst = rawData();
1443  for (i = 0; i < n; i++)
1444  {
1445  *dst++ = T(*srcdata++);
1446  }
1447  }
1448  else
1449  {
1450  T *dst = rawData();
1451 
1452  srcoff = 0;
1453  for (i = 0; i < n; i++)
1454  {
1455  dst[i] = T(srcdata[srcoff]);
1456  srcoff += srcstride;
1457  }
1458  }
1459 }
1460 
1461 template <typename T>
1462 void
1463 UT_VoxelTile<T>::setValue(int x, int y, int z, T t)
1464 {
1465  UT_ASSERT_P(x >= 0 && y >= 0 && z >= 0);
1466  UT_ASSERT_P(x < myRes[0] && y < myRes[1] && z < myRes[2]);
1467 
1468  // Determine if assignment is compatible with current
1469  // compression technique.
1470  if (writeThrough(x, y, z, t))
1471  {
1472  return;
1473  }
1474  // Can't write to our current type of tile with the
1475  // given value, so abandon.
1476  uncompress();
1477 
1478  // Attempt to write through again. This must succeed
1479  // as we should now be uncompressed.
1480  UT_VERIFY_P(writeThrough(x, y, z, t));
1481 }
1482 
1483 template <typename T>
1484 void
1486 {
1487  switch (myCompressionType)
1488  {
1489  case COMPRESS_RAW:
1490  // Trivial.
1491  return;
1492 
1493  case COMPRESS_CONSTANT:
1494  {
1495  // Must expand the tile!
1496  T cval;
1497  int i, n;
1498 
1499  cval = rawConstVal();
1500  freeData();
1501 
1502  myCompressionType = COMPRESS_RAW;
1503 
1504  if (myRes[0] == TILESIZE &&
1505  myRes[1] == TILESIZE)
1506  {
1507  // Flag that we can do fast lookups on this tile.
1508  myCompressionType = COMPRESS_RAWFULL;
1509  }
1510 
1511  n = myRes[0] * myRes[1] * myRes[2];
1512  myData = UT_VOXEL_ALLOC(sizeof(T) * n);
1513 
1514  for (i = 0; i < n; i++)
1515  {
1516  ((T *)myData)[i] = cval;
1517  }
1518  return;
1519  }
1520  case COMPRESS_RAWFULL:
1521  {
1522  T *raw;
1523  int x, y, z, i, n;
1524 
1525  if (myRes[0] == TILESIZE &&
1526  myRes[1] == TILESIZE)
1527  {
1528  // Trivial
1529  return;
1530  }
1531 
1532  // Need to contract to the actual compact size.
1533  n = myRes[0] * myRes[1] * myRes[2];
1534  raw = (T *)UT_VOXEL_ALLOC(sizeof(T) * n);
1535  i = 0;
1536  for (z = 0; z < myRes[2]; z++)
1537  {
1538  for (y = 0; y < myRes[1]; y++)
1539  {
1540  for (x = 0; x < myRes[0]; x++)
1541  {
1542  raw[i++] = ((T *)myData)[x+(y+z*TILESIZE)*TILESIZE];
1543  }
1544  }
1545  }
1546 
1547  freeData();
1548  myCompressionType = COMPRESS_RAW;
1549  myData = raw;
1550 
1551  return;
1552  }
1553 
1554  case COMPRESS_FPREAL16:
1555  {
1556  T *raw;
1557  int x, y, z, i, n;
1558  fpreal16 *src = (fpreal16 *) myData;
1559 
1560  n = myRes[0] * myRes[1] * myRes[2];
1561  raw = (T *)UT_VOXEL_ALLOC(sizeof(T) * n);
1562  i = 0;
1563  for (z = 0; z < myRes[2]; z++)
1564  {
1565  for (y = 0; y < myRes[1]; y++)
1566  {
1567  for (x = 0; x < myRes[0]; x++)
1568  {
1569  raw[i] = src[i];
1570  i++;
1571  }
1572  }
1573  }
1574  freeData();
1575  myCompressionType = COMPRESS_RAW;
1576  myData = raw;
1577  return;
1578  }
1579  }
1580 
1581  // Use the compression engine.
1582  UT_VoxelTileCompress<T> *engine;
1583 
1584  engine = getCompressionEngine(myCompressionType);
1585 
1586  // We use the read ability to set our values.
1587  int x, y, z, i;
1588  T *raw;
1589 
1590  raw = (T *) UT_VOXEL_ALLOC(sizeof(T) * myRes[0] * myRes[1] * myRes[2]);
1591  i = 0;
1592  for (z = 0; z < myRes[2]; z++)
1593  {
1594  for (y = 0; y < myRes[1]; y++)
1595  {
1596  for (x = 0; x < myRes[0]; x++)
1597  {
1598  raw[i++] = engine->getValue(*this, x, y, z);
1599  }
1600  }
1601  }
1602 
1603  freeData();
1604 
1605  // Now copy in the result
1606  myCompressionType = COMPRESS_RAW;
1607  if (myRes[0] == TILESIZE &&
1608  myRes[1] == TILESIZE)
1609  {
1610  // Flag that we can do fast lookups on this tile.
1611  myCompressionType = COMPRESS_RAWFULL;
1612  }
1613 
1614  myData = raw;
1615 }
1616 
1617 template <typename T>
1618 void
1620 {
1621  T *raw;
1622  int x, y, z, i;
1623 
1624  if (myCompressionType == COMPRESS_RAWFULL)
1625  return;
1626 
1627  uncompress();
1628 
1629  UT_ASSERT(myCompressionType == COMPRESS_RAW);
1630 
1631  // The RAWFULL format only differs from RAW when the tile dimensions
1632  // are smaller than the maximum tile size.
1633  if (myRes[0] < TILESIZE || myRes[1] < TILESIZE)
1634  {
1635  raw = (T *)UT_VOXEL_ALLOC(sizeof(T) * TILESIZE * TILESIZE * myRes[2]);
1636  i = 0;
1637  for (z = 0; z < myRes[2]; z++)
1638  {
1639  for (y = 0; y < myRes[1]; y++)
1640  {
1641  for (x = 0; x < myRes[0]; x++)
1642  {
1643  raw[x+(y+z*TILESIZE)*TILESIZE] = ((T *)myData)[i++];
1644  }
1645  }
1646  }
1647  freeData();
1648  myData = raw;
1649  }
1650  myCompressionType = COMPRESS_RAWFULL;
1651 }
1652 
1653 template <typename T>
1654 void
1656 {
1657  if (isRaw() || isRawFull())
1658  return;
1659 
1660  freeData();
1661 
1662  if (myRes[0] == TILESIZE && myRes[1] == TILESIZE)
1663  myCompressionType = COMPRESS_RAWFULL;
1664  else
1665  myCompressionType = COMPRESS_RAW;
1666 
1667  myData = UT_VOXEL_ALLOC(sizeof(T) * numVoxels());
1668 }
1669 
1670 template <typename T>
1671 void
1673 {
1674  float irx, iry, irz;
1675 
1676  irx = 1.0 / myRes[0];
1677  iry = 1.0 / myRes[1];
1678  irz = 1.0 / myRes[2];
1679  switch (myCompressionType)
1680  {
1681  case COMPRESS_RAW:
1682  {
1683  int i;
1684  const T *data = (const T*) myData;
1685 
1686  i = 0;
1687  T zavg;
1688  zavg = 0;
1689  for (int z = 0; z < myRes[2]; z++)
1690  {
1691  T yavg;
1692  yavg = 0;
1693  for (int y = 0; y < myRes[1]; y++)
1694  {
1695  T xavg;
1696  xavg = 0;
1697  for (int x = 0; x < myRes[0]; x++)
1698  {
1699  xavg += data[i++];
1700  }
1701  xavg *= irx;
1702  yavg += xavg;
1703  }
1704  yavg *= iry;
1705  zavg += yavg;
1706  }
1707  zavg *= irz;
1708 
1709  avg = zavg;
1710  return;
1711  }
1712 
1713  case COMPRESS_FPREAL16:
1714  {
1715  int i;
1716  const fpreal16 *data = (fpreal16 *) myData;
1717 
1718  i = 0;
1719  fpreal16 zavg = 0;
1720  for (int z = 0; z < myRes[2]; z++)
1721  {
1722  fpreal16 yavg = 0;
1723  for (int y = 0; y < myRes[1]; y++)
1724  {
1725  fpreal16 xavg = 0;
1726  for (int x = 0; x < myRes[0]; x++)
1727  {
1728  xavg += data[i++];
1729  }
1730  xavg *= irx;
1731  yavg += xavg;
1732  }
1733  yavg *= iry;
1734  zavg += yavg;
1735  }
1736  zavg *= irz;
1737 
1738  avg = zavg;
1739  return;
1740  }
1741 
1742  case COMPRESS_CONSTANT:
1743  avg = rawConstVal();
1744  return;
1745 
1746  case COMPRESS_RAWFULL:
1747  {
1748  int offset;
1749  const T *data = (const T*) myData;
1750 
1751  offset = 0;
1752  T zavg;
1753  zavg = 0;
1754  for (int z = 0; z < myRes[2]; z++)
1755  {
1756  T yavg;
1757  yavg = 0;
1758  for (int y = 0; y < myRes[1]; y++)
1759  {
1760  T xavg;
1761  xavg = 0;
1762  for (int x = 0; x < myRes[0]; x++)
1763  {
1764  xavg += data[x+offset];
1765  }
1766  xavg *= irx;
1767  yavg += xavg;
1768  offset += TILESIZE;
1769  }
1770  yavg *= iry;
1771  zavg += yavg;
1772  }
1773  zavg *= irz;
1774 
1775  avg = zavg;
1776  return;
1777  }
1778 
1779  default:
1780  {
1781  T zavg;
1782  zavg = 0;
1783  for (int z = 0; z < myRes[2]; z++)
1784  {
1785  T yavg;
1786  yavg = 0;
1787  for (int y = 0; y < myRes[1]; y++)
1788  {
1789  T xavg;
1790  xavg = 0;
1791  for (int x = 0; x < myRes[0]; x++)
1792  {
1793  xavg += (*this)(x, y, z);
1794  }
1795  xavg *= irx;
1796  yavg += xavg;
1797  }
1798  yavg *= iry;
1799  zavg += yavg;
1800  }
1801  zavg *= irz;
1802 
1803  avg = zavg;
1804  return;
1805  }
1806  }
1807 }
1808 
1809 template <typename T>
1810 void
1812 {
1813  switch (myCompressionType)
1814  {
1815  case COMPRESS_RAW:
1816  {
1817  int n = myRes[0] * myRes[1] * myRes[2];
1818  int i;
1819 
1820  min = max = *(T*)myData;
1821  for (i = 1; i < n; i++)
1822  {
1823  expandMinMax( ((T*)myData)[i], min, max );
1824  }
1825  return;
1826  }
1827 
1828  case COMPRESS_FPREAL16:
1829  {
1830  int n = myRes[0] * myRes[1] * myRes[2];
1831  int i;
1832  fpreal16 *src = (fpreal16 *)myData;
1833 
1834  min = max = *src;
1835  for (i = 1; i < n; i++)
1836  {
1837  T val;
1838  val = src[i];
1839  expandMinMax( val, min, max );
1840  }
1841  return;
1842  }
1843 
1844  case COMPRESS_CONSTANT:
1845  min = max = rawConstVal();
1846  return;
1847 
1848  case COMPRESS_RAWFULL:
1849  {
1850  int x, y, z, offset;
1851 
1852  min = max = *(T*)myData;
1853  offset = 0;
1854  for (z = 0; z < myRes[2]; z++)
1855  {
1856  for (y = 0; y < myRes[1]; y++)
1857  {
1858  for (x = 0; x < myRes[0]; x++)
1859  {
1860  expandMinMax(
1861  ((T*)myData)[x+offset],
1862  min, max );
1863  }
1864  offset += TILESIZE;
1865  }
1866  }
1867  return;
1868  }
1869 
1870  default:
1871  {
1872  // Use the compression engine.
1873  UT_VoxelTileCompress<T> *engine;
1874 
1875  engine = getCompressionEngine(myCompressionType);
1876 
1877  engine->findMinMax(*this, min, max);
1878  return;
1879  }
1880  }
1881 }
1882 
1883 template <typename T>
1884 bool
1886 {
1887  switch (myCompressionType)
1888  {
1889  case COMPRESS_RAW:
1890  case COMPRESS_RAWFULL:
1891  {
1892  int n = myRes[0] * myRes[1] * myRes[2];
1893  int i;
1894 
1895  for (i = 0; i < n; i++)
1896  {
1897  if (SYSisNan( ((T*)myData)[i] ))
1898  return true;
1899  }
1900  return false;
1901  }
1902 
1903  case COMPRESS_FPREAL16:
1904  {
1905  return false;
1906  }
1907 
1908  case COMPRESS_CONSTANT:
1909  if (SYSisNan(rawConstVal()))
1910  return true;
1911  return false;
1912 
1913  default:
1914  {
1915  // Use the compression engine.
1916  UT_VoxelTileCompress<T> *engine;
1917  int x, y, z;
1918 
1919  engine = getCompressionEngine(myCompressionType);
1920 
1921  for (z = 0; z < myRes[2]; z++)
1922  for (y = 0; y < myRes[1]; y++)
1923  for (x = 0; x < myRes[0]; x++)
1924  if (SYSisNan(engine->getValue(*this, x, y, z)))
1925  {
1926  return true;
1927  }
1928 
1929  return false;
1930  }
1931  }
1932 }
1933 
1934 template <typename T>
1935 bool
1937 {
1938  T min, max;
1939  bool losslessonly = (options.myQuantizeTol == 0.0);
1940  int i;
1941  UT_VoxelTileCompress<T> *engine;
1942 
1943  // This is as easy as it gets.
1944  if (myCompressionType == COMPRESS_CONSTANT)
1945  return false;
1946 
1947  findMinMax(min, max);
1948 
1949  // See if we can be made into a constant tile.
1950  if (dist(min, max) <= options.myConstantTol)
1951  {
1952  // Ignore if we are already constant.
1953  if (myCompressionType == COMPRESS_CONSTANT)
1954  return false;
1955 
1956  // We are fully constant!
1957  if (min != max)
1958  {
1959  T avg;
1960  findAverage(avg);
1961  makeConstant(avg);
1962  }
1963  else
1964  {
1965  makeConstant(min);
1966  }
1967  return true;
1968  }
1969 
1970  bool result = false;
1971 
1972  if (options.myAllowFP16)
1973  {
1974  if (myCompressionType == COMPRESS_RAW ||
1975  myCompressionType == COMPRESS_RAWFULL)
1976  {
1977  makeFpreal16();
1978  result = true;
1979  }
1980  }
1981 
1982  for (i = 0; i < getCompressionEngines().entries(); i++)
1983  {
1984  engine = getCompressionEngines()(i);
1985 
1986  // Ignore possibly expensive lossy compressions.
1987  if (losslessonly && !engine->isLossless())
1988  continue;
1989 
1990  if (engine->tryCompress(*this, options, min, max))
1991  {
1992  myCompressionType = i + COMPRESS_ENGINE;
1993  result = true;
1994  // We keep testing in case another compression engine
1995  // can get a better number...
1996  }
1997  }
1998 
1999  // If we are RAW compress, check to see if we could become
2000  // RAWFULL for faster access.
2001  if (myCompressionType == COMPRESS_RAW)
2002  {
2003  if (myRes[0] == TILESIZE && myRes[1] == TILESIZE)
2004  {
2005  myCompressionType = COMPRESS_RAWFULL;
2006  }
2007  }
2008 
2009  // No suitable compression found.
2010  return result;
2011 }
2012 
2013 template <typename T>
2014 void
2016 {
2017  if (!isConstant())
2018  {
2019  freeData();
2020 
2021  if (!inlineConstant())
2022  myData = UT_VOXEL_ALLOC(sizeof(T));
2023  }
2024 
2025  myCompressionType = COMPRESS_CONSTANT;
2026  *rawConstData() = t;
2027 }
2028 
2029 template <typename T>
2030 void
2032 {
2033  if (isConstant())
2034  {
2035  return;
2036  }
2037 
2038  if (myCompressionType == COMPRESS_FPREAL16)
2039  return;
2040 
2041  // Get our new data.
2042  int len = myRes[2] * myRes[1] * myRes[0];
2043  fpreal16 *data = (fpreal16 *)UT_VOXEL_ALLOC(sizeof(fpreal16) * len);
2044 
2045  if (myCompressionType == COMPRESS_RAW ||
2046  myCompressionType == COMPRESS_RAWFULL)
2047  {
2048  for (int i = 0; i < len; i++)
2049  {
2050  data[i] = UTvoxelConvertFP16( ((T*)myData)[i] );
2051  }
2052  }
2053  else
2054  {
2055  // Apply any converters.
2056  int i = 0;
2057 
2058  for (int z = 0; z < myRes[2]; z++)
2059  {
2060  for (int y = 0; y < myRes[1]; y++)
2061  {
2062  for (int x = 0; x < myRes[0]; x++)
2063  {
2064  data[i++] = UTvoxelConvertFP16( (*this)(x, y, z));
2065  }
2066  }
2067  }
2068  }
2069 
2070  freeData();
2071  myData = data;
2072  myCompressionType = COMPRESS_FPREAL16;
2073 }
2074 
2075 template <typename T>
2076 int64
2077 UT_VoxelTile<T>::getMemoryUsage(bool inclusive) const
2078 {
2079  int64 mem = inclusive ? sizeof(*this) : 0;
2080  mem += getDataLength();
2081  return mem;
2082 }
2083 
2084 template <typename T>
2085 int64
2087 {
2088  exint usage;
2089 
2090  switch (myCompressionType)
2091  {
2092  case COMPRESS_RAW:
2093  usage = sizeof(T) * xres() * yres() * zres();
2094  break;
2095  case COMPRESS_FPREAL16:
2096  usage = sizeof(fpreal16) * xres() * yres() * zres();
2097  break;
2098  case COMPRESS_CONSTANT:
2099  if (inlineConstant())
2100  usage = 0;
2101  else
2102  usage = sizeof(T);
2103  break;
2104  case COMPRESS_RAWFULL:
2105  usage = sizeof(T) * TILESIZE * TILESIZE * zres();
2106  break;
2107  default:
2108  {
2109  // Use the compression engine.
2110  UT_VoxelTileCompress<T> *engine;
2111  engine = getCompressionEngine(myCompressionType);
2112  usage = engine->getDataLength(*this);
2113  break;
2114  }
2115  }
2116  return usage;
2117 }
2118 
2119 template <typename T>
2120 void
2121 UT_VoxelTile<T>::weightedSum(int pstart[3], int pend[3],
2122  const float *weights[3], int start[3],
2123  T &result)
2124 {
2125  int ix, iy, iz, i;
2126  int px, py, pz;
2127  int ixstart, ixend;
2128  int tstart[3];
2129  fpreal w, pw;
2130  T psumx, psumy;
2131 
2132  switch (myCompressionType)
2133  {
2134  case COMPRESS_CONSTANT:
2135  {
2136  w = 1;
2137  for (i = 0; i < 3; i++)
2138  {
2139  pw = 0;
2140  for (ix = 0; ix < pend[i]-pstart[i]; ix++)
2141  pw += weights[i][ix+pstart[i]-start[i]];
2142  w *= pw;
2143  }
2144  result += w * rawConstVal();
2145  break;
2146  }
2147 
2148  case COMPRESS_RAW:
2149  {
2150  tstart[0] = pstart[0] & TILEMASK;
2151  tstart[1] = pstart[1] & TILEMASK;
2152  tstart[2] = pstart[2] & TILEMASK;
2153  ixstart = pstart[0]-start[0];
2154  ixend = pend[0]-start[0];
2155  pz = tstart[2];
2156  UT_ASSERT(pz < myRes[2]);
2157  UT_ASSERT(ixend - ixstart <= myRes[0]);
2158  for (iz = pstart[2]; iz < pend[2]; iz++, pz++)
2159  {
2160  psumy = 0;
2161  py = tstart[1];
2162  UT_ASSERT(py < myRes[1]);
2163  for (iy = pstart[1]; iy < pend[1]; iy++, py++)
2164  {
2165  psumx = 0;
2166  px = ((pz * myRes[1]) + py) * myRes[0] + tstart[0];
2167  for (ix = ixstart; ix < ixend; ix++, px++)
2168  {
2169  psumx += weights[0][ix]* ((T*)myData)[px];
2170  }
2171  psumy += weights[1][iy-start[1]] * psumx;
2172  }
2173  result += weights[2][iz-start[2]] * psumy;
2174  }
2175  break;
2176  }
2177 
2178  case COMPRESS_FPREAL16:
2179  {
2180  fpreal16 *src = (fpreal16 *) myData;
2181 
2182  tstart[0] = pstart[0] & TILEMASK;
2183  tstart[1] = pstart[1] & TILEMASK;
2184  tstart[2] = pstart[2] & TILEMASK;
2185  ixstart = pstart[0]-start[0];
2186  ixend = pend[0]-start[0];
2187  pz = tstart[2];
2188  UT_ASSERT(pz < myRes[2]);
2189  UT_ASSERT(ixend - ixstart <= myRes[0]);
2190  for (iz = pstart[2]; iz < pend[2]; iz++, pz++)
2191  {
2192  psumy = 0;
2193  py = tstart[1];
2194  UT_ASSERT(py < myRes[1]);
2195  for (iy = pstart[1]; iy < pend[1]; iy++, py++)
2196  {
2197  psumx = 0;
2198  px = ((pz * myRes[1]) + py) * myRes[0] + tstart[0];
2199  for (ix = ixstart; ix < ixend; ix++, px++)
2200  {
2201  T val;
2202  val = src[px];
2203  psumx += weights[0][ix]* val;
2204  }
2205  psumy += weights[1][iy-start[1]] * psumx;
2206  }
2207  result += weights[2][iz-start[2]] * psumy;
2208  }
2209  break;
2210  }
2211  case COMPRESS_RAWFULL:
2212  {
2213  tstart[0] = pstart[0] & TILEMASK;
2214  tstart[1] = pstart[1] & TILEMASK;
2215  tstart[2] = pstart[2] & TILEMASK;
2216  ixstart = pstart[0]-start[0];
2217  ixend = pend[0]-start[0];
2218  pz = tstart[2];
2219  for (iz = pstart[2]; iz < pend[2]; iz++, pz++)
2220  {
2221  psumy = 0;
2222  py = tstart[1];
2223  for (iy = pstart[1]; iy < pend[1]; iy++, py++)
2224  {
2225  psumx = 0;
2226  px = ((pz * TILESIZE) + py) * TILESIZE + tstart[0];
2227  for (ix = ixstart; ix < ixend; ix++, px++)
2228  {
2229  psumx += weights[0][ix]* ((T*)myData)[px];
2230  }
2231  psumy += weights[1][iy-start[1]] * psumx;
2232  }
2233  result += weights[2][iz-start[2]] * psumy;
2234  }
2235  break;
2236  }
2237 
2238  default:
2239  {
2240  // For all other compression types, we use our
2241  // getValue() accessor rather than trying to
2242  // do anything fancy.
2243  tstart[0] = pstart[0] & TILEMASK;
2244  tstart[1] = pstart[1] & TILEMASK;
2245  tstart[2] = pstart[2] & TILEMASK;
2246  ixstart = pstart[0]-start[0];
2247  ixend = pend[0]-start[0];
2248  pz = tstart[2];
2249  for (iz = pstart[2]; iz < pend[2]; iz++, pz++)
2250  {
2251  psumy = 0;
2252  py = tstart[1];
2253  for (iy = pstart[1]; iy < pend[1]; iy++, py++)
2254  {
2255  psumx = 0;
2256  px = tstart[0];
2257  for (ix = ixstart; ix < ixend; ix++, px++)
2258  {
2259  psumx += weights[0][ix] *
2260  (*this)(px, py, pz);
2261  }
2262  psumy += weights[1][iy-start[1]] * psumx;
2263  }
2264  result += weights[2][iz-start[2]] * psumy;
2265  }
2266  break;
2267  }
2268  }
2269 }
2270 
2271 template <typename T>
2272 void
2274 {
2275  int i;
2276 
2277  // Repalce old copy, if any.
2278  for (i = 0; i < getCompressionEngines().entries(); i++)
2279  {
2280  if (!strcmp(engine->getName(), getCompressionEngines()(i)->getName()))
2281  {
2282  getCompressionEngines()(i) = engine;
2283  return;
2284  }
2285  }
2286 
2287  getCompressionEngines().append(engine);
2288 }
2289 
2290 template <typename T>
2291 int
2293 {
2294  int i;
2295 
2296  if (!name)
2297  return -1;
2298 
2299  for (i = 0; i < getCompressionEngines().entries(); i++)
2300  {
2301  if (!strcmp(name, getCompressionEngines()(i)->getName()))
2302  {
2303  return i + COMPRESS_ENGINE;
2304  }
2305  }
2306 
2307  return -1;
2308 }
2309 
2310 template <typename T>
2313 {
2314  index -= COMPRESS_ENGINE;
2315 
2316  return getCompressionEngines()(index);
2317 }
2318 
2319 template <typename T>
2320 void
2321 UT_VoxelTile<T>::save(std::ostream &os) const
2322 {
2323  // Always save in our native format...
2324  if (myCompressionType >= COMPRESS_ENGINE)
2325  {
2326  UT_VoxelTileCompress<T> *engine = getCompressionEngine(myCompressionType);
2327 
2328  if (engine->canSave())
2329  {
2330  char type = myCompressionType;
2331 
2332  UTwrite(os, &type, 1);
2333  engine->save(os, *this);
2334  return;
2335  }
2336 
2337  // Can't save, must save as raw.
2338  char type = COMPRESS_RAW;
2339  T value;
2340 
2341  UTwrite(os, &type, 1);
2342 
2343  for (int z = 0; z < zres(); z++)
2344  for (int y = 0; y < yres(); y++)
2345  for (int x = 0; x < xres(); x++)
2346  {
2347  value = (*this)(x, y, z);
2348  UTwrite<T>(os, &value, 1);
2349  }
2350  return;
2351  }
2352 
2353  int len;
2354  char type = myCompressionType;
2355 
2356  if (type == COMPRESS_RAWFULL &&
2357  myRes[0] == TILESIZE &&
2358  myRes[1] == TILESIZE &&
2359  myRes[2] != TILESIZE)
2360  {
2361  // Forbid saving this as a raw full as that will confuse
2362  // older versions of Houdini.
2363  type = COMPRESS_RAW;
2364  }
2365 
2366  UT_ASSERT(type >= 0 && type < COMPRESS_ENGINE);
2367 
2368  UTwrite(os, &type, 1);
2369 
2370  switch ((CompressionType) type)
2371  {
2372  case COMPRESS_RAW:
2373  len = myRes[2] * myRes[1] * myRes[0];
2374  UTwrite<T>(os, (T *) myData, len);
2375  break;
2376 
2377  case COMPRESS_FPREAL16:
2378  len = myRes[2] * myRes[1] * myRes[0];
2379  UTwrite(os, (int16 *) myData, len);
2380  break;
2381 
2382  case COMPRESS_RAWFULL:
2383  len = TILESIZE * TILESIZE * TILESIZE;
2384  UTwrite<T>(os, (T *) myData, len);
2385  break;
2386 
2387  case COMPRESS_CONSTANT:
2388  UTwrite<T>(os, rawConstData(), 1);
2389  break;
2390 
2391  case COMPRESS_ENGINE:
2392  UT_ASSERT(!"Invalid compression type");
2393  break;
2394  }
2395 }
2396 
2397 template <typename T>
2398 void
2400 {
2401  char type, otype;
2402  int len;
2403 
2404  is.readChar(type);
2405 
2406  // Perform the indirection to find out our native type.
2407  if (type >= 0 && type < compress.entries())
2408  {
2409  otype = type;
2410  type = compress(type);
2411  if (type == -1)
2412  {
2413  std::cerr << "Missing compression engine " << (int) otype << "\n";
2414  }
2415  }
2416 
2417  if (type >= COMPRESS_ENGINE)
2418  {
2419  freeData();
2420  myCompressionType = type;
2421 
2422  if (type - COMPRESS_ENGINE >= getCompressionEngines().entries())
2423  {
2424  // Invalid type!
2425  std::cerr << "Invalid compression engine " << (int) otype << "\n";
2426  return;
2427  }
2428 
2429  UT_VoxelTileCompress<T> *engine = getCompressionEngine(myCompressionType);
2430 
2431  engine->load(is, *this);
2432 
2433  return;
2434  }
2435 
2436  UT_ASSERT(type >= 0);
2437  if (type < 0)
2438  {
2439  return;
2440  }
2441 
2442  freeData();
2443 
2444  myCompressionType = type;
2445 
2446  switch ((CompressionType) myCompressionType)
2447  {
2448  case COMPRESS_RAW:
2449  len = myRes[2] * myRes[1] * myRes[0];
2450  myData = UT_VOXEL_ALLOC(sizeof(T) * len);
2451  is.read<T>((T *) myData, len);
2452  break;
2453 
2454  case COMPRESS_FPREAL16:
2455  len = myRes[2] * myRes[1] * myRes[0];
2456  myData = UT_VOXEL_ALLOC(sizeof(fpreal16) * len);
2457  is.read((int16 *) myData, len);
2458  break;
2459 
2460  case COMPRESS_RAWFULL:
2461  len = TILESIZE * TILESIZE * TILESIZE;
2462  myData = UT_VOXEL_ALLOC(sizeof(T) * len);
2463  is.read<T>((T *) myData, len);
2464  break;
2465 
2466  case COMPRESS_CONSTANT:
2467  if (!inlineConstant())
2468  myData = UT_VOXEL_ALLOC(sizeof(T));
2469  is.read<T>(rawConstData(), 1);
2470  break;
2471 
2472  case COMPRESS_ENGINE:
2473  UT_ASSERT(!"Invalid compression type");
2474  break;
2475  }
2476 
2477  // Recompress raw full that are not complete in z
2478  if (myCompressionType == COMPRESS_RAW &&
2479  myRes[0] == TILESIZE &&
2480  myRes[1] == TILESIZE)
2481  myCompressionType = COMPRESS_RAWFULL;
2482 }
2483 
2484 template <typename T>
2485 bool
2487 {
2488  bool ok = true;
2489 
2490  // Always save in our native format...
2491  if (myCompressionType >= COMPRESS_ENGINE)
2492  {
2493  UT_VoxelTileCompress<T> *engine = getCompressionEngine(myCompressionType);
2494 
2495  if (engine->canSave())
2496  {
2497  char type = myCompressionType;
2498 
2499  ok = ok && w.jsonBeginArray();
2500  ok = ok && w.jsonKey(UT_VoxelArrayJSON::getToken(
2502  ok = ok && w.jsonInt(type);
2503  ok = ok && w.jsonKey(UT_VoxelArrayJSON::getToken(
2505  ok = ok && engine->save(w, *this);
2506  ok = ok && w.jsonEndArray();
2507  return ok;
2508  }
2509 
2510  // Can't save, must save as raw.
2511  char type = COMPRESS_RAW;
2512  T value;
2513 
2514  ok = ok && w.jsonBeginArray();
2515  ok = ok && w.jsonKey(UT_VoxelArrayJSON::getToken(
2517  ok = ok && w.jsonInt(type);
2518  ok = ok && w.jsonKey(UT_VoxelArrayJSON::getToken(
2520 
2521  ok = ok && w.beginUniformArray(xres()*yres()*zres(),
2522  w.jidFromValue(&value));
2523  for (int z = 0; z < zres(); z++)
2524  for (int y = 0; y < yres(); y++)
2525  for (int x = 0; x < xres(); x++)
2526  {
2527  value = (*this)(x, y, z);
2528  ok = ok && w.uniformWrite(value);
2529  }
2530  ok = ok && w.endUniformArray();
2531  ok = ok && w.jsonEndArray();
2532  return ok;
2533  }
2534 
2535  int len;
2536  char type = myCompressionType;
2537 
2538  UT_ASSERT(type >= 0 && type < COMPRESS_ENGINE);
2539  ok = ok && w.jsonBeginArray();
2540  ok = ok && w.jsonKey(UT_VoxelArrayJSON::getToken(
2542  ok = ok && w.jsonInt(type);
2543  ok = ok && w.jsonKey(UT_VoxelArrayJSON::getToken(
2545 
2546  switch (myCompressionType)
2547  {
2548  case COMPRESS_RAW:
2549  len = myRes[2] * myRes[1] * myRes[0];
2550  ok = ok && w.jsonUniformArray(len, (T *)myData);
2551  break;
2552 
2553  case COMPRESS_FPREAL16:
2554  len = myRes[2] * myRes[1] * myRes[0];
2555  ok = ok && w.jsonUniformArray(len, (fpreal16 *)myData);
2556  break;
2557 
2558  case COMPRESS_RAWFULL:
2559  len = TILESIZE * TILESIZE * myRes[2];
2560  ok = ok && w.jsonUniformArray(len, (T *)myData);
2561  break;
2562 
2563  case COMPRESS_CONSTANT:
2564  ok = ok && w.jsonValue(rawConstVal());
2565  break;
2566  }
2567 
2568  ok = ok && w.jsonEndArray();
2569  return ok;
2570 }
2571 
2572 template <typename T>
2573 bool
2575 {
2577  UT_WorkBuffer key;
2578  int8 type, otype;
2579  int len;
2580 
2581  it = p.beginArray();
2582  if (it.atEnd() || !it.getLowerKey(key))
2583  return false;
2584  if (UT_VoxelArrayJSON::getTileID(key.buffer()) !=
2586  return false;
2587  if (!p.parseNumber(type))
2588  return false;
2589 
2590  ++it;
2591 
2592  if (it.atEnd() || !it.getLowerKey(key))
2593  return false;
2594  if (UT_VoxelArrayJSON::getTileID(key.buffer()) !=
2596  return false;
2597 
2598  // Perform the indirection to find out our native type.
2599  if (type >= 0 && type < compress.entries())
2600  {
2601  otype = type;
2602  type = compress(type);
2603  if (type == -1)
2604  {
2605  std::cerr << "Missing compression engine " << (int) otype << "\n";
2606  }
2607  }
2608 
2609  if (type >= COMPRESS_ENGINE)
2610  {
2611  freeData();
2612  myCompressionType = type;
2613 
2614  if (type - COMPRESS_ENGINE >= getCompressionEngines().entries())
2615  {
2616  // Invalid type!
2617  std::cerr << "Invalid compression engine " << (int) otype << "\n";
2618  return false;
2619  }
2620 
2621  UT_VoxelTileCompress<T> *engine = getCompressionEngine(myCompressionType);
2622 
2623  bool engine_ok = engine->load(p, *this);
2624  ++it;
2625  return it.atEnd() && engine_ok;
2626  }
2627 
2628  UT_ASSERT(type >= 0);
2629  if (type < 0)
2630  {
2631  return false;
2632  }
2633 
2634  freeData();
2635 
2636  myCompressionType = type;
2637 
2638  switch (myCompressionType)
2639  {
2640  case COMPRESS_RAW:
2641  len = myRes[2] * myRes[1] * myRes[0];
2642  myData = UT_VOXEL_ALLOC(sizeof(T) * len);
2643  if (p.parseUniformArray((T *)myData, len) != len)
2644  return false;
2645  break;
2646 
2647  case COMPRESS_FPREAL16:
2648  len = myRes[2] * myRes[1] * myRes[0];
2649  myData = UT_VOXEL_ALLOC(sizeof(fpreal16) * len);
2650  if (p.parseUniformArray((fpreal16 *)myData, len) != len)
2651  return false;
2652  break;
2653 
2654  case COMPRESS_RAWFULL:
2655  len = TILESIZE * TILESIZE * myRes[2];
2656  myData = UT_VOXEL_ALLOC(sizeof(T) * len);
2657  if (p.parseUniformArray((T *)myData, len) != len)
2658  return false;
2659  break;
2660 
2661  case COMPRESS_CONSTANT:
2662  if (!inlineConstant())
2663  myData = UT_VOXEL_ALLOC(sizeof(T));
2664  if (!p.parseNumber(*rawConstData()))
2665  return false;
2666  break;
2667  }
2668 
2669  ++it;
2670  return it.atEnd();
2671 }
2672 
2673 template <typename T>
2674 void
2676 {
2677  int16 ntype = getCompressionEngines().entries();
2678  int i;
2679 
2680  ntype += COMPRESS_ENGINE;
2681 
2682  UTwrite(os, &ntype);
2683 
2684  UT_ASSERT(COMPRESS_ENGINE == 4); // Change lower lines!
2686  UTsaveStringBinary(os, "rawfull", UT_STRING_8BIT_IO);
2687  UTsaveStringBinary(os, "constant", UT_STRING_8BIT_IO);
2688  UTsaveStringBinary(os, "fpreal16", UT_STRING_8BIT_IO);
2689 
2690  ntype -= COMPRESS_ENGINE;
2691  for (i = 0; i < ntype; i++)
2692  {
2693  UTsaveStringBinary(os, getCompressionEngines()(i)->getName(), UT_STRING_8BIT_IO);
2694  }
2695 }
2696 
2697 template <typename T>
2698 void
2700 {
2701  int16 ntype;
2702  int i, idx;
2703 
2704  compress.entries(0);
2705 
2706  is.read(&ntype);
2707 
2708  for (i = 0; i < ntype; i++)
2709  {
2710  UT_String name;
2711 
2713  if (name == "raw")
2714  compress.append(COMPRESS_RAW);
2715  else if (name == "rawfull")
2716  compress.append(COMPRESS_RAWFULL);
2717  else if (name == "constant")
2718  compress.append(COMPRESS_CONSTANT);
2719  else if (name == "fpreal16")
2720  compress.append(COMPRESS_FPREAL16);
2721  else
2722  {
2723  idx = lookupCompressionEngine(name);
2724 
2725  // -1 means a failure to find it in our set..
2726  // this is only a bad thing if a tile actually uses this engine.
2727 
2728  compress.append(idx);
2729  }
2730  }
2731 }
2732 
2733 template <typename T>
2734 bool
2736 {
2737  int16 ntype = getCompressionEngines().entries();
2738  int i;
2739  bool ok = true;
2740 
2741  ntype += COMPRESS_ENGINE;
2742 
2743  UT_ASSERT(COMPRESS_ENGINE == 4); // Change lower lines!
2744  ok = ok && w.beginUniformArray(ntype, UT_JID_STRING);
2745  ok = ok && w.uniformWrite("raw");
2746  ok = ok && w.uniformWrite("rawfull");
2747  ok = ok && w.uniformWrite("constant");
2748  ok = ok && w.uniformWrite("fpreal16");
2749 
2750  ntype -= COMPRESS_ENGINE;
2751  for (i = 0; i < ntype; i++)
2752  {
2753  ok = ok && w.uniformWrite(getCompressionEngines()(i)->getName());
2754  }
2755 
2756  ok = ok && w.endUniformArray();
2757  return ok;
2758 }
2759 
2760 template <typename T>
2761 bool
2763 {
2766  int idx;
2767 
2768  compress.entries(0);
2769  for (it = p.beginArray(); !it.atEnd(); ++it)
2770  {
2771  if (!p.parseString(buffer))
2772  return false;
2773 
2774  if (!buffer.strcmp("raw"))
2775  compress.append(COMPRESS_RAW);
2776  else if (!buffer.strcmp("rawfull"))
2777  compress.append(COMPRESS_RAWFULL);
2778  else if (!buffer.strcmp("constant"))
2779  compress.append(COMPRESS_CONSTANT);
2780  else if (!buffer.strcmp("fpreal16"))
2781  compress.append(COMPRESS_FPREAL16);
2782  else
2783  {
2784  idx = lookupCompressionEngine(buffer.buffer());
2785 
2786  // -1 means a failure to find it in our set..
2787  // this is only a bad thing if a tile actually uses this engine.
2788 
2789  compress.append(idx);
2790  }
2791  }
2792  return true;
2793 }
2794 
2795 
2796 
2797 //
2798 // VoxelArray definitions.
2799 //
2800 
2801 template <typename T>
2803 {
2804  myRes[0] = 0;
2805  myRes[1] = 0;
2806  myRes[2] = 0;
2807  myTileRes[0] = 0;
2808  myTileRes[1] = 0;
2809  myTileRes[2] = 0;
2810 
2811  myTiles = 0;
2812 
2813  myInvRes = 1;
2814 
2815  myBorderValue = 0;
2816  myBorderScale[0] = 0;
2817  myBorderScale[1] = 0;
2818  myBorderScale[2] = 0;
2819  myBorderType = UT_VOXELBORDER_STREAK;
2820 
2821  mySharedMem = 0;
2822  mySharedMemView = 0;
2823 }
2824 
2825 template <typename T>
2827 {
2828  deleteVoxels();
2829 }
2830 
2831 template <typename T>
2833 {
2834  myRes[0] = 0;
2835  myRes[1] = 0;
2836  myRes[2] = 0;
2837  myInvRes = 1;
2838  myTileRes[0] = 0;
2839  myTileRes[1] = 0;
2840  myTileRes[2] = 0;
2841 
2842  myTiles = 0;
2843 
2844  mySharedMem = 0;
2845  mySharedMemView = 0;
2846 
2847  *this = src;
2848 }
2849 
2850 template <typename T>
2851 void
2853  const UT_JobInfo &info)
2854 {
2855  UT_ASSERT(isMatching(src));
2856  UT_VoxelArrayIterator<T> vit(this);
2857  vit.splitByTile(info);
2858  for (vit.rewind(); !vit.atEnd(); vit.advanceTile())
2859  {
2860  int i = vit.getLinearTileNum();
2861  myTiles[i] = src.myTiles[i];
2862  }
2863 }
2864 
2865 template <typename T>
2866 const UT_VoxelArray<T> &
2868 {
2869  // Paranoid code:
2870  if (&src == this)
2871  return *this;
2872 
2873  // Clear out my own data (this is done again by size, but not
2874  // maybe with resize?)
2875  deleteVoxels();
2876 
2877  myBorderScale[0] = src.myBorderScale[0];
2878  myBorderScale[1] = src.myBorderScale[1];
2879  myBorderScale[2] = src.myBorderScale[2];
2880  myBorderValue = src.myBorderValue;
2881  myBorderType = src.myBorderType;
2882 
2883  myCompressionOptions = src.myCompressionOptions;
2884 
2885  // Allocate proper size.
2886  size(src.myRes[0], src.myRes[1], src.myRes[2]);
2887 
2888  copyData(src);
2889 
2890  return *this;
2891 }
2892 
2893 template <typename T>
2894 void
2896 {
2897  delete [] myTiles;
2898  myTiles = 0;
2899 
2900  myTileRes[0] = myTileRes[1] = myTileRes[2] = 0;
2901  myRes[0] = myRes[1] = myRes[2] = 0;
2902 
2903  delete mySharedMemView;
2904  mySharedMemView = 0;
2905 
2906  delete mySharedMem;
2907  mySharedMem = 0;
2908 }
2909 
2910 template <typename T>
2911 void
2912 UT_VoxelArray<T>::size(int xres, int yres, int zres)
2913 {
2914  deleteVoxels();
2915 
2916  // Initialize our tiles.
2917  int tile_res[3];
2918  tile_res[0] = (xres + TILEMASK) >> TILEBITS;
2919  tile_res[1] = (yres + TILEMASK) >> TILEBITS;
2920  tile_res[2] = (zres + TILEMASK) >> TILEBITS;
2921 
2922  exint ntiles = ((exint)tile_res[0]) * tile_res[1] * tile_res[2];
2923  if (ntiles)
2924  {
2925  myTiles = new UT_VoxelTile<T> [ntiles];
2926 
2927  // Only set resolutions *AFTER* we successfully allocate
2928  myRes[0] = xres;
2929  myRes[1] = yres;
2930  myRes[2] = zres;
2931 
2932  myInvRes = 1;
2933  if (xres)
2934  myInvRes.x() = 1.0f / myRes[0];
2935  if (yres)
2936  myInvRes.y() = 1.0f / myRes[1];
2937  if (zres)
2938  myInvRes.z() = 1.0f / myRes[2];
2939 
2940  myTileRes[0] = tile_res[0];
2941  myTileRes[1] = tile_res[1];
2942  myTileRes[2] = tile_res[2];
2943 
2944  int i = 0;
2945  for (int tz = 0; tz < myTileRes[2]; tz++)
2946  {
2947  int zr;
2948  if (tz < myTileRes[2]-1)
2949  zr = TILESIZE;
2950  else
2951  zr = zres - tz * TILESIZE;
2952 
2953  for (int ty = 0; ty < myTileRes[1]; ty++)
2954  {
2955  int yr;
2956  if (ty < myTileRes[1]-1)
2957  yr = TILESIZE;
2958  else
2959  yr = yres - ty * TILESIZE;
2960 
2961  int tx, xr = TILESIZE;
2962  for (tx = 0; tx < myTileRes[0]-1; tx++)
2963  {
2964  myTiles[i].setRes(xr, yr, zr);
2965 
2966  i++;
2967  }
2968  xr = xres - tx * TILESIZE;
2969  myTiles[i].setRes(xr, yr, zr);
2970  i++;
2971  }
2972  }
2973  }
2974  else
2975  myTiles = 0;
2976 }
2977 
2978 template <typename T>
2979 void
2981 {
2982  // Only resize if our res is different.
2983  if (src.getXRes() != getXRes() ||
2984  src.getYRes() != getYRes() ||
2985  src.getZRes() != getZRes())
2986  {
2987  size(src.getXRes(), src.getYRes(), src.getZRes());
2988  }
2989 
2990  // Update border conditions.
2991  myBorderType = src.myBorderType;
2992  myBorderScale[0] = src.myBorderScale[0];
2993  myBorderScale[1] = src.myBorderScale[1];
2994  myBorderScale[2] = src.myBorderScale[2];
2995  myBorderValue = src.myBorderValue;
2996 
2997  // Match our compression tolerances
2998  myCompressionOptions = src.myCompressionOptions;
2999 }
3000 
3001 template <typename T>
3002 int64
3004 {
3005  int64 mem = inclusive ? sizeof(*this) : 0;
3006 
3007  int ntiles = numTiles();
3008  for (int i = 0; i < ntiles; i++)
3009  mem += myTiles[i].getMemoryUsage(true);
3010 
3011  if (mySharedMem)
3012  mem += mySharedMem->getMemoryUsage(true);
3013 
3014  if (mySharedMemView)
3015  mem += mySharedMemView->getMemoryUsage(true);
3016 
3017  return mem;
3018 }
3019 
3020 template <typename T>
3021 void
3023 {
3024  UT_VoxelArrayIterator<T> vit(this);
3025  vit.splitByTile(info);
3026  for (vit.rewind(); !vit.atEnd(); vit.advanceTile())
3027  {
3028  int i = vit.getLinearTileNum();
3029  myTiles[i].makeConstant(t);
3030  }
3031 }
3032 
3033 template <typename T>
3034 bool
3036 {
3037  int i, ntiles;
3038  T cval = 0;
3039  const fpreal tol = SYS_FTOLERANCE_R;
3040 
3041  ntiles = numTiles();
3042  for (i = 0; i < ntiles; i++)
3043  {
3044  if (!myTiles[i].isConstant())
3045  {
3046  return false;
3047  }
3048 
3049  if (!i)
3050  {
3051  // First tile, get the constant value.
3052  cval = myTiles[i].rawConstVal();
3053  }
3054  else
3055  {
3056  // See if we have deviated too much.
3057  if (UT_VoxelTile<T>::dist(cval, myTiles[i].rawConstVal()) > tol)
3058  {
3059  return false;
3060  }
3061  }
3062  }
3063 
3064  // All tiles are both constant and within tolerance of each
3065  // other. Write out our constant value and return true.
3066  if (t)
3067  *t = cval;
3068 
3069  return true;
3070 }
3071 
3072 template <typename T>
3073 bool
3075 {
3076  int i, ntiles;
3077 
3078  ntiles = numTiles();
3079  for (i = 0; i < ntiles; i++)
3080  {
3081  if (myTiles[i].hasNan())
3082  {
3083  return true;
3084  }
3085  }
3086 
3087  return false;
3088 }
3089 
3090 template <typename T>
3091 T
3093 {
3094  // We go from the position in the unit cube into the index.
3095  // The center of cells must map to the exact integer indices.
3096  pos.x() *= myRes[0];
3097  pos.y() *= myRes[1];
3098  pos.z() *= myRes[2];
3099  pos.x() -= 0.5;
3100  pos.y() -= 0.5;
3101  pos.z() -= 0.5;
3102 
3103  return lerpVoxelCoord(pos);
3104 }
3105 
3106 template <typename T>
3107 T
3109 {
3110  int x, y, z;
3111  // Yes, these have to be 32 becaues split float requires 32!
3112  fpreal32 fx, fy, fz;
3113 
3114  splitVoxelCoord(pos, x, y, z, fx, fy, fz);
3115 
3116  return lerpVoxel(x, y, z, fx, fy, fz);
3117 }
3118 
3119 template <typename T>
3120 template <int AXIS2D>
3121 T
3123 {
3124  int x, y, z;
3125  // Yes, these have to be 32 becaues split float requires 32!
3126  fpreal32 fx, fy, fz;
3127 
3128  splitVoxelCoordAxis<AXIS2D>(pos, x, y, z, fx, fy, fz);
3129 
3130  return lerpVoxelAxis<AXIS2D>(x, y, z, fx, fy, fz);
3131 }
3132 
3133 template <typename T>
3134 T
3135 UT_VoxelArray<T>::lerpVoxel(int x, int y, int z,
3136  float fx, float fy, float fz) const
3137 {
3138  // Do trilinear interpolation.
3139  T vx, vx1, vy, vy1, vz;
3140 
3141  // Optimization for common cases (values are within the voxel range and
3142  // are all within the same tile)
3143  if ( !((x | y | z) < 0) &&
3144  (((x - myRes[0]+1) & (y - myRes[1]+1) & (z - myRes[2]+1)) < 0))
3145 
3146  // (x > 0) && (y > 0) && (z > 0) &&
3147  // Do not use x+1 < foo as if x is MAX_INT that will falsely
3148  // report in bounds!
3149  // (x < myRes[0]-1) && (y < myRes[1]-1) && (z < myRes[2]-1) )
3150  {
3151  int xm, ym, zm;
3152 
3153  xm = x & TILEMASK;
3154  ym = y & TILEMASK;
3155  zm = z & TILEMASK;
3156 
3157  if ((xm != TILEMASK) && (ym != TILEMASK) && (zm != TILEMASK))
3158  {
3159  const UT_VoxelTile<T> *tile =
3160  getTile(x >> TILEBITS, y >> TILEBITS, z >> TILEBITS);
3161 
3162  vz = tile->lerp(xm, ym, zm, fx, fy, fz);
3163  }
3164  else
3165  {
3166  // We cross tile boundaries but we remain within
3167  // the voxel grid. We can thus avoid any bounds checking
3168  // and use operator() rather than getValue.
3169 
3170  // Lerp x:x+1, y, z
3171  vx = UT_VoxelTile<T>::lerpValues((*this)(x, y, z),
3172  (*this)(x+1, y, z),
3173  fx);
3174  // Lerp x:x+1, y+1, z
3175  vx1= UT_VoxelTile<T>::lerpValues((*this)(x, y+1, z),
3176  (*this)(x+1, y+1, z),
3177  fx);
3178  // Lerp x:x+1, y:y+1, z
3179  vy = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);
3180 
3181  // Lerp x:x+1, y, z+1
3182  vx = UT_VoxelTile<T>::lerpValues((*this)(x, y, z+1),
3183  (*this)(x+1, y, z+1),
3184  fx);
3185  // Lerp x:x+1, y+1, z+1
3186  vx1= UT_VoxelTile<T>::lerpValues((*this)(x, y+1, z+1),
3187  (*this)(x+1, y+1, z+1),
3188  fx);
3189 
3190  // Lerp x:x+1, y:y+1, z+1
3191  vy1 = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);
3192 
3193  // Lerp x:x+1, y:y+1, z:z+1
3194  vz = UT_VoxelTile<T>::lerpValues(vy, vy1, fz);
3195  }
3196  }
3197  else
3198  {
3199  // Lerp x:x+1, y, z
3200  vx = UT_VoxelTile<T>::lerpValues(getValue(x, y, z),
3201  getValue(x+1, y, z),
3202  fx);
3203  // Lerp x:x+1, y+1, z
3204  vx1= UT_VoxelTile<T>::lerpValues(getValue(x, y+1, z),
3205  getValue(x+1, y+1, z),
3206  fx);
3207 
3208  // Lerp x:x+1, y:y+1, z
3209  vy = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);
3210 
3211  // Lerp x:x+1, y, z+1
3212  vx = UT_VoxelTile<T>::lerpValues(getValue(x, y, z+1),
3213  getValue(x+1, y, z+1),
3214  fx);
3215  // Lerp x:x+1, y+1, z+1
3216  vx1= UT_VoxelTile<T>::lerpValues(getValue(x, y+1, z+1),
3217  getValue(x+1, y+1, z+1),
3218  fx);
3219 
3220  // Lerp x:x+1, y:y+1, z+1
3221  vy1 = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);
3222 
3223  // Lerp x:x+1, y:y+1, z:z+1
3224  vz = UT_VoxelTile<T>::lerpValues(vy, vy1, fz);
3225  }
3226 
3227  return vz;
3228 }
3229 
3230 template <typename T>
3231 template <int AXIS2D>
3232 T
3234  float fx, float fy, float fz) const
3235 {
3236  // Do bilinear interpolation.
3237  T vx, vx1, vy, vy1, vz;
3238 
3239  int lesscomp = 0, greatercomp = -1;
3240 
3241  if (AXIS2D != 0)
3242  {
3243  lesscomp |= x;
3244  greatercomp &= (x - myRes[0]+1);
3245  }
3246  if (AXIS2D != 1)
3247  {
3248  lesscomp |= y;
3249  greatercomp &= (y - myRes[1]+1);
3250  }
3251  if (AXIS2D != 2)
3252  {
3253  lesscomp |= z;
3254  greatercomp &= (z - myRes[2]+1);
3255  }
3256 
3257  // Optimization for common cases (values are within the voxel range and
3258  // are all within the same tile)
3259  if ( !(lesscomp < 0) && (greatercomp < 0) )
3260  {
3261  int xm, ym, zm;
3262 
3263  xm = x & TILEMASK;
3264  ym = y & TILEMASK;
3265  zm = z & TILEMASK;
3266 
3267  if ((AXIS2D == 0 || xm != TILEMASK) &&
3268  (AXIS2D == 1 || ym != TILEMASK) &&
3269  (AXIS2D == 2 || zm != TILEMASK))
3270  {
3271  const UT_VoxelTile<T> *tile =
3272  getTile( (AXIS2D == 0) ? 0 : (x >> TILEBITS),
3273  (AXIS2D == 1) ? 0 : (y >> TILEBITS),
3274  (AXIS2D == 2) ? 0 : (z >> TILEBITS) );
3275 
3276  vz = tile->template lerpAxis<AXIS2D>(xm, ym, zm, fx, fy, fz);
3277  }
3278  else
3279  {
3280  // We cross tile boundaries but we remain within
3281  // the voxel grid. We can thus avoid any bounds checking
3282  // and use operator() rather than getValue.
3283 
3284  // Lerp x:x+1, y, z
3285  if (AXIS2D != 0)
3286  vx = UT_VoxelTile<T>::lerpValues((*this)(x, y, z),
3287  (*this)(x+1, y, z),
3288  fx);
3289  else
3290  vx = (*this)(x, y, z);
3291 
3292  if (AXIS2D != 1)
3293  {
3294  // Lerp x:x+1, y+1, z
3295  if (AXIS2D != 0)
3296  vx1= UT_VoxelTile<T>::lerpValues((*this)(x, y+1, z),
3297  (*this)(x+1, y+1, z),
3298  fx);
3299  else
3300  vx1 = (*this)(x, y+1, z);
3301  // Lerp x:x+1, y:y+1, z
3302  vy = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);
3303  }
3304  else
3305  vy = vx;
3306 
3307  if (AXIS2D != 2)
3308  {
3309  // Lerp x:x+1, y, z+1
3310  if (AXIS2D != 0)
3311  vx = UT_VoxelTile<T>::lerpValues((*this)(x, y, z+1),
3312  (*this)(x+1, y, z+1),
3313  fx);
3314  else
3315  vx = (*this)(x, y, z+1);
3316 
3317  if (AXIS2D != 1)
3318  {
3319  // Lerp x:x+1, y+1, z+1
3320  if (AXIS2D != 0)
3321  vx1= UT_VoxelTile<T>::lerpValues((*this)(x, y+1, z+1),
3322  (*this)(x+1, y+1, z+1),
3323  fx);
3324  else
3325  vx1 = (*this)(x, y+1, z+1);
3326  // Lerp x:x+1, y:y+1, z+1
3327  vy1 = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);
3328  }
3329  else
3330  vy1 = vx;
3331 
3332  // Lerp x:x+1, y:y+1, z:z+1
3333  vz = UT_VoxelTile<T>::lerpValues(vy, vy1, fz);
3334  }
3335  else
3336  vz = vy;
3337  }
3338  }
3339  else
3340  {
3341  // Lerp x:x+1, y, z
3342  if (AXIS2D != 0)
3343  vx = UT_VoxelTile<T>::lerpValues(getValue(x, y, z),
3344  getValue(x+1, y, z),
3345  fx);
3346  else
3347  vx = getValue(x, y, z);
3348 
3349  if (AXIS2D != 1)
3350  {
3351  // Lerp x:x+1, y+1, z
3352  if (AXIS2D != 0)
3353  vx1= UT_VoxelTile<T>::lerpValues(getValue(x, y+1, z),
3354  getValue(x+1, y+1, z),
3355  fx);
3356  else
3357  vx1 = getValue(x, y+1, z);
3358  // Lerp x:x+1, y:y+1, z
3359  vy = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);
3360  }
3361  else
3362  vy = vx;
3363 
3364  if (AXIS2D != 2)
3365  {
3366  // Lerp x:x+1, y, z+1
3367  if (AXIS2D != 0)
3368  vx = UT_VoxelTile<T>::lerpValues(getValue(x, y, z+1),
3369  getValue(x+1, y, z+1),
3370  fx);
3371  else
3372  vx = getValue(x, y, z+1);
3373 
3374  if (AXIS2D != 1)
3375  {
3376  // Lerp x:x+1, y+1, z+1
3377  if (AXIS2D != 0)
3378  vx1= UT_VoxelTile<T>::lerpValues(getValue(x, y+1, z+1),
3379  getValue(x+1, y+1, z+1),
3380  fx);
3381  else
3382  vx1 = getValue(x, y+1, z+1);
3383  // Lerp x:x+1, y:y+1, z+1
3384  vy1 = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);
3385  }
3386  else
3387  vy1 = vx;
3388 
3389  // Lerp x:x+1, y:y+1, z:z+1
3390  vz = UT_VoxelTile<T>::lerpValues(vy, vy1, fz);
3391  }
3392  else
3393  vz = vy;
3394  }
3395 
3396  return vz;
3397 }
3398 
3399 template <typename T>
3400 void
3402  UT_Vector3F pos) const
3403 {
3404  // We go from the position in the unit cube into the index.
3405  // The center of cells must map to the exact integer indices.
3406  pos.x() *= myRes[0];
3407  pos.y() *= myRes[1];
3408  pos.z() *= myRes[2];
3409  pos.x() -= 0.5;
3410  pos.y() -= 0.5;
3411  pos.z() -= 0.5;
3412 
3413  lerpVoxelCoordMinMax(lerp, lmin, lmax, pos);
3414 }
3415 
3416 template <typename T>
3417 void
3419  UT_Vector3F pos) const
3420 {
3421  int x, y, z;
3422  // Yes, these have to be 32 becaues split float requires 32!
3423  fpreal32 fx, fy, fz;
3424 
3425  splitVoxelCoord(pos, x, y, z, fx, fy, fz);
3426 
3427  lerpVoxelMinMax(lerp, lmin, lmax, x, y, z, fx, fy, fz);
3428 }
3429 
3430 
3431 template <typename T>
3432 template <int AXIS2D>
3433 void
3435  UT_Vector3F pos) const
3436 {
3437  int x, y, z;
3438  // Yes, these have to be 32 becaues split float requires 32!
3439  fpreal32 fx, fy, fz;
3440 
3441  splitVoxelCoordAxis<AXIS2D>(pos, x, y, z, fx, fy, fz);
3442 
3443  lerpVoxelMinMaxAxis<AXIS2D>(lerp, lmin, lmax, x, y, z, fx, fy, fz);
3444 }
3445 
3446 
3447 template <typename T>
3448 void
3450  T &lerp, T &smin, T &smax,
3451  int x, int y, int z,
3452  float fx, float fy, float fz) const
3453 {
3454  T samples[8];
3455 
3456  if (extractSample(x, y, z, samples))
3457  {
3458  lerp = smin = smax = samples[0];
3459  return;
3460  }
3461 
3462  lerp = lerpSample(samples, fx, fy, fz);
3463 
3464  T smin1, smax1;
3465 
3466  SYSminmax(samples[0], samples[1], samples[2], samples[3],
3467  smin, smax);
3468  SYSminmax(samples[4+0], samples[4+1], samples[4+2], samples[4+3],
3469  smin1, smax1);
3470 
3471  smin = SYSmin(smin, smin1);
3472  smax = SYSmax(smax, smax1);
3473 }
3474 
3475 template <typename T>
3476 template <int AXIS2D>
3477 void
3479  T &lerp, T &smin, T &smax,
3480  int x, int y, int z,
3481  float fx, float fy, float fz) const
3482 {
3483  T samples[8];
3484 
3485  if (extractSampleAxis<AXIS2D>(x, y, z, samples))
3486  {
3487  lerp = smin = smax = samples[0];
3488  return;
3489  }
3490 
3491  lerp = lerpSampleAxis<AXIS2D>(samples, fx, fy, fz);
3492 
3493  if (AXIS2D == 0)
3494  SYSminmax(samples[0], samples[2], samples[4], samples[6],
3495  smin, smax);
3496  else if (AXIS2D == 1)
3497  SYSminmax(samples[0], samples[1], samples[4], samples[5],
3498  smin, smax);
3499  else if (AXIS2D == 2)
3500  SYSminmax(samples[0], samples[1], samples[2], samples[3],
3501  smin, smax);
3502  else
3503  {
3504  T smin1, smax1;
3505 
3506  SYSminmax(samples[0], samples[1], samples[2], samples[3],
3507  smin, smax);
3508  SYSminmax(samples[4+0], samples[4+1], samples[4+2], samples[4+3],
3509  smin1, smax1);
3510 
3511  smin = SYSmin(smin, smin1);
3512  smax = SYSmax(smax, smax1);
3513  }
3514 }
3515 
3516 template <typename T>
3517 bool
3519  T *samples) const
3520 {
3521  // Optimization for common cases (values are within the voxel range and
3522  // are all within the same tile)
3523  if ( !((x | y | z) < 0) &&
3524  (((x - myRes[0]+1) & (y - myRes[1]+1) & (z - myRes[2]+1)) < 0))
3525 
3526  // (x > 0) && (y > 0) && (z > 0) &&
3527  // Do not use x+1 < foo as if x is MAX_INT that will falsely
3528  // report in bounds!
3529  // (x < myRes[0]-1) && (y < myRes[1]-1) && (z < myRes[2]-1) )
3530  {
3531  int xm, ym, zm;
3532 
3533  xm = x & TILEMASK;
3534  ym = y & TILEMASK;
3535  zm = z & TILEMASK;
3536 
3537  if ((xm != TILEMASK) && (ym != TILEMASK) && (zm != TILEMASK))
3538  {
3539  const UT_VoxelTile<T> *tile =
3540  getTile(x >> TILEBITS, y >> TILEBITS, z >> TILEBITS);
3541 
3542  return tile->extractSample(xm, ym, zm, samples);
3543  }
3544  else
3545  {
3546  // We cross tile boundaries but we remain within
3547  // the voxel grid. We can thus avoid any bounds checking
3548  // and use operator() rather than getValue.
3549  samples[0] = (*this)(x, y, z);
3550  samples[1] = (*this)(x+1, y, z);
3551  samples[2+0] = (*this)(x, y+1, z);
3552  samples[2+1] = (*this)(x+1, y+1, z);
3553  samples[4+0] = (*this)(x, y, z+1);
3554  samples[4+1] = (*this)(x+1, y, z+1);
3555  samples[4+2+0] = (*this)(x, y+1, z+1);
3556  samples[4+2+1] = (*this)(x+1, y+1, z+1);
3557  }
3558  }
3559  else
3560  {
3561  samples[0] = getValue(x, y, z);
3562  samples[1] = getValue(x+1, y, z);
3563  samples[2+0] = getValue(x, y+1, z);
3564  samples[2+1] = getValue(x+1, y+1, z);
3565  samples[4+0] = getValue(x, y, z+1);
3566  samples[4+1] = getValue(x+1, y, z+1);
3567  samples[4+2+0] = getValue(x, y+1, z+1);
3568  samples[4+2+1] = getValue(x+1, y+1, z+1);
3569  }
3570 
3571  return false;
3572 }
3573 
3574 
3575 template <typename T>
3576 template <int AXIS2D>
3577 bool
3579  T *samples) const
3580 {
3581  // Optimization for common cases (values are within the voxel range and
3582  // are all within the same tile)
3583  int lesscomp = 0, greatercomp = -1;
3584 
3585  if (AXIS2D != 0)
3586  {
3587  lesscomp |= x;
3588  greatercomp &= (x - myRes[0]+1);
3589  }
3590  if (AXIS2D != 1)
3591  {
3592  lesscomp |= y;
3593  greatercomp &= (y - myRes[1]+1);
3594  }
3595  if (AXIS2D != 2)
3596  {
3597  lesscomp |= z;
3598  greatercomp &= (z - myRes[2]+1);
3599  }
3600 
3601  // Optimization for common cases (values are within the voxel range and
3602  // are all within the same tile)
3603  if ( !(lesscomp < 0) && (greatercomp < 0) )
3604  {
3605  int xm, ym, zm;
3606 
3607  xm = x & TILEMASK;
3608  ym = y & TILEMASK;
3609  zm = z & TILEMASK;
3610 
3611  if ((AXIS2D == 0 || xm != TILEMASK) &&
3612  (AXIS2D == 1 || ym != TILEMASK) &&
3613  (AXIS2D == 2 || zm != TILEMASK))
3614  {
3615  const UT_VoxelTile<T> *tile =
3616  getTile( (AXIS2D == 0) ? 0 : (x >> TILEBITS),
3617  (AXIS2D == 1) ? 0 : (y >> TILEBITS),
3618  (AXIS2D == 2) ? 0 : (z >> TILEBITS) );
3619 
3620  return tile->template extractSampleAxis<AXIS2D>(xm, ym, zm, samples);
3621  }
3622  else
3623  {
3624  // We cross tile boundaries but we remain within
3625  // the voxel grid. We can thus avoid any bounds checking
3626  // and use operator() rather than getValue.
3627  samples[0] = (*this)(x, y, z);
3628  if (AXIS2D != 0)
3629  samples[1] = (*this)(x+1, y, z);
3630  if (AXIS2D != 1)
3631  {
3632  samples[2+0] = (*this)(x, y+1, z);
3633  if (AXIS2D != 0)
3634  samples[2+1] = (*this)(x+1, y+1, z);
3635  }
3636  if (AXIS2D != 2)
3637  {
3638  samples[4+0] = (*this)(x, y, z+1);
3639  if (AXIS2D != 0)
3640  samples[4+1] = (*this)(x+1, y, z+1);
3641  if (AXIS2D != 1)
3642  {
3643  samples[4+2+0] = (*this)(x, y+1, z+1);
3644  if (AXIS2D != 0)
3645  samples[4+2+1] = (*this)(x+1, y+1, z+1);
3646  }
3647  }
3648  }
3649  }
3650  else
3651  {
3652  samples[0] = getValue(x, y, z);
3653  if (AXIS2D != 0)
3654  samples[1] = getValue(x+1, y, z);
3655  if (AXIS2D != 1)
3656  {
3657  samples[2+0] = getValue(x, y+1, z);
3658  if (AXIS2D != 0)
3659  samples[2+1] = getValue(x+1, y+1, z);
3660  }
3661  if (AXIS2D != 2)
3662  {
3663  samples[4+0] = getValue(x, y, z+1);
3664  if (AXIS2D != 0)
3665  samples[4+1] = getValue(x+1, y, z+1);
3666  if (AXIS2D != 1)
3667  {
3668  samples[4+2+0] = getValue(x, y+1, z+1);
3669  if (AXIS2D != 0)
3670  samples[4+2+1] = getValue(x+1, y+1, z+1);
3671  }
3672  }
3673  }
3674 
3675  return false;
3676 }
3677 
3678 template <typename T>
3679 bool
3681  T *samples) const
3682 {
3683  // Optimization for common cases (values are within the voxel range and
3684  // are all within the same tile)
3685  if ( !(((x-1) | (y-1) | (z-1)) < 0) &&
3686  (((x - myRes[0]+1) & (y - myRes[1]+1) & (z - myRes[2]+1)) < 0))
3687 
3688  // (x > 0) && (y > 0) && (z > 0) &&
3689  // Do not use x+1 < foo as if x is MAX_INT that will falsely
3690  // report in bounds!
3691  // (x < myRes[0]-1) && (y < myRes[1]-1) && (z < myRes[2]-1) )
3692  {
3693  int xm, ym, zm;
3694 
3695  xm = x & TILEMASK;
3696  ym = y & TILEMASK;
3697  zm = z & TILEMASK;
3698 
3699  if (xm && ym && zm && (xm != TILEMASK) && (ym != TILEMASK) && (zm != TILEMASK))
3700  {
3701  const UT_VoxelTile<T> *tile =
3702  getTile(x >> TILEBITS, y >> TILEBITS, z >> TILEBITS);
3703 
3704  return tile->extractSamplePlus(xm, ym, zm, samples);
3705  }
3706  else
3707  {
3708  // We cross tile boundaries but we remain within
3709  // the voxel grid. We can thus avoid any bounds checking
3710  // and use operator() rather than getValue.
3711  samples[0] = (*this)(x-1, y, z);
3712  samples[1] = (*this)(x+1, y, z);
3713  samples[2+0] = (*this)(x, y-1, z);
3714  samples[2+1] = (*this)(x, y+1, z);
3715  samples[4+0] = (*this)(x, y, z-1);
3716  samples[4+1] = (*this)(x, y, z+1);
3717  samples[6] = (*this)(x, y, z);
3718  }
3719  }
3720  else
3721  {
3722  samples[0] = getValue(x-1, y, z);
3723  samples[1] = getValue(x+1, y, z);
3724  samples[2+0] = getValue(x, y-1, z);
3725  samples[2+1] = getValue(x, y+1, z);
3726  samples[4+0] = getValue(x, y, z-1);
3727  samples[4+1] = getValue(x, y, z+1);
3728  samples[6] = getValue(x, y, z);
3729  }
3730 
3731  return false;
3732 }
3733 
3734 template <typename T>
3735 bool
3737  T *samples) const
3738 {
3739  // Optimization for common cases (values are within the voxel range and
3740  // are all within the same tile)
3741  if ( !(((x-1) | (y-1) | (z-1)) < 0) &&
3742  (((x - myRes[0]+1) & (y - myRes[1]+1) & (z - myRes[2]+1)) < 0))
3743 
3744  // (x > 0) && (y > 0) && (z > 0) &&
3745  // Do not use x+1 < foo as if x is MAX_INT that will falsely
3746  // report in bounds!
3747  // (x < myRes[0]-1) && (y < myRes[1]-1) && (z < myRes[2]-1) )
3748  {
3749  int xm, ym, zm;
3750 
3751  xm = x & TILEMASK;
3752  ym = y & TILEMASK;
3753  zm = z & TILEMASK;
3754 
3755  if (xm && ym && zm && (xm != TILEMASK) && (ym != TILEMASK) && (zm != TILEMASK))
3756  {
3757  const UT_VoxelTile<T> *tile =
3758  getTile(x >> TILEBITS, y >> TILEBITS, z >> TILEBITS);
3759 
3760  return tile->extractSampleCube(xm, ym, zm, samples);
3761  }
3762  else
3763  {
3764  // We cross tile boundaries but we remain within
3765  // the voxel grid. We can thus avoid any bounds checking
3766  // and use operator() rather than getValue.
3767  int sampidx = 0;
3768  for (int dz = -1; dz <= 1; dz++)
3769  {
3770  for (int dy = -1; dy <= 1; dy++)
3771  {
3772  for (int dx = -1; dx <= 1; dx++)
3773  {
3774  samples[sampidx++] = (*this)(x+dx, y+dy, z+dz);
3775  }
3776  }
3777  }
3778  }
3779  }
3780  else
3781  {
3782  int sampidx = 0;
3783  for (int dz = -1; dz <= 1; dz++)
3784  {
3785  for (int dy = -1; dy <= 1; dy++)
3786  {
3787  for (int dx = -1; dx <= 1; dx++)
3788  {
3789  samples[sampidx++] = getValue(x+dx, y+dy, z+dz);
3790  }
3791  }
3792  }
3793  }
3794 
3795  return false;
3796 }
3797 template <typename T>
3798 T
3800  float fx, float fy, float fz) const
3801 {
3802 #if 1
3803  // Do trilinear interpolation.
3804  T vx, vx1, vy, vy1, vz;
3805 
3806  // Lerp x:x+1, y, z
3807  vx = UT_VoxelTile<T>::lerpValues(samples[0],
3808  samples[1],
3809  fx);
3810  // Lerp x:x+1, y+1, z
3811  vx1= UT_VoxelTile<T>::lerpValues(samples[2],
3812  samples[2+1],
3813  fx);
3814  // Lerp x:x+1, y:y+1, z
3815  vy = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);
3816 
3817  // Lerp x:x+1, y, z+1
3818  vx = UT_VoxelTile<T>::lerpValues(samples[4],
3819  samples[4+1],
3820  fx);
3821  // Lerp x:x+1, y+1, z+1
3822  vx1= UT_VoxelTile<T>::lerpValues(samples[4+2],
3823  samples[4+2+1],
3824  fx);
3825 
3826  // Lerp x:x+1, y:y+1, z+1
3827  vy1 = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);
3828 
3829  // Lerp x:x+1, y:y+1, z:z+1
3830  vz = UT_VoxelTile<T>::lerpValues(vy, vy1, fz);
3831 
3832  return vz;
3833 #else
3834  v4uf a, b, vfx, vfy, vfz;
3835 
3836  a = v4uf(&samples[0]);
3837  b = v4uf(&samples[4]);
3838 
3839  vfx = v4uf(fx);
3840  vfy = v4uf(fy);
3841  vfz = v4uf(fz);
3842 
3843  b -= a;
3844  a = madd(b, fz, a);
3845 
3846  b = a.swizzle<2, 3, 0, 1>();
3847  b -= a;
3848  a = madd(b, fy, a);
3849 
3850  b = a.swizzle<1, 2, 3, 0>();
3851  b -= a;
3852  a = madd(b, fx, a);
3853 
3854  return a[0];
3855 #endif
3856 }
3857 
3858 template <typename T>
3859 template <int AXIS2D>
3860 T
3862  float fx, float fy, float fz) const
3863 {
3864  // Do trilinear interpolation.
3865  T vx, vx1, vy, vy1, vz;
3866 
3867  // Lerp x:x+1, y, z
3868  if (AXIS2D != 0)
3869  vx = UT_VoxelTile<T>::lerpValues(samples[0], samples[1], fx);
3870  else
3871  vx = samples[0];
3872 
3873  if (AXIS2D != 1)
3874  {
3875  // Lerp x:x+1, y+1, z
3876  if (AXIS2D != 0)
3877  vx1= UT_VoxelTile<T>::lerpValues(samples[2], samples[2+1], fx);
3878  else
3879  vx1= samples[2];
3880 
3881  // Lerp x:x+1, y:y+1, z
3882  vy = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);
3883  }
3884  else
3885  vy = vx;
3886 
3887  if (AXIS2D != 2)
3888  {
3889  // Lerp x:x+1, y, z+1
3890  if (AXIS2D != 0)
3891  vx = UT_VoxelTile<T>::lerpValues(samples[4], samples[4+1], fx);
3892  else
3893  vx = samples[4];
3894 
3895  if (AXIS2D != 1)
3896  {
3897  // Lerp x:x+1, y+1, z+1
3898  if (AXIS2D != 0)
3899  vx1= UT_VoxelTile<T>::lerpValues(samples[4+2], samples[4+2+1], fx);
3900  else
3901  vx1= samples[4+2];
3902 
3903  // Lerp x:x+1, y:y+1, z+1
3904  vy1 = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);
3905  }
3906  else
3907  vy1 = vx;
3908 
3909  // Lerp x:x+1, y:y+1, z:z+1
3910  vz = UT_VoxelTile<T>::lerpValues(vy, vy1, fz);
3911  }
3912  else
3913  vz = vy;
3914 
3915  return vz;
3916 }
3917 
3918 template <typename T>
3919 T
3921 {
3922  int x, y, z;
3923  // Yes, these have to be 32 becaues split float requires 32!
3924  fpreal32 fx, fy, fz;
3925 
3926  // We go from the position in the unit cube into the index.
3927  // The center of cells must map to the exact integer indices.
3928  pos.x() *= myRes[0];
3929  pos.y() *= myRes[1];
3930  pos.z() *= myRes[2];
3931  pos.x() -= 0.5;
3932  pos.y() -= 0.5;
3933  pos.z() -= 0.5;
3934 
3935  // Determine integer & fractional components.
3936  fx = pos.x();
3937  SYSfastSplitFloat(fx, x);
3938  fy = pos.y();
3939  SYSfastSplitFloat(fy, y);
3940  fz = pos.z();
3941  SYSfastSplitFloat(fz, z);
3942 
3943  // Do trilinear interpolation.
3944  T vx, vx1, vy, vy1, vz;
3945 
3946  // NOTE:
3947  // If you get a crash reading one of these getValues, note
3948  // that if you are reading from the same voxel array that you
3949  // are writing to, and doing so in a multi-threaded, tile-by-tile
3950  // approach, just because your positions match doesn't mean you
3951  // won't access neighbouring tiles. To avoid this, perform
3952  // the obvious optimization of checking for aligned tiles
3953  // and using direct indices.
3954  // (I have avoided using if (!fx && !fy && !fz) as the test
3955  // as that is both potentially inaccurate (round off in conversion
3956  // may mess us up) and said optimzation is a useful one anyways)
3957 
3958  // Optimization for common cases (values are within the voxel range and
3959  // are all within the same tile)
3960  if ( !((x | y | z) < 0) &&
3961  (((x - myRes[0]+1) & (y - myRes[1]+1) & (z - myRes[2]+1)) < 0))
3962  // if ( (x > 0) && (y > 0) && (z > 0) &&
3963  // Do not use x+1 < foo as if x is MAX_INT that will falsely
3964  // report in bounds!
3965  // (x < myRes[0]-1) && (y < myRes[1]-1) && (z < myRes[2]-1) )
3966  {
3967  int xm, ym, zm;
3968 
3969  xm = x & TILEMASK;
3970  ym = y & TILEMASK;
3971  zm = z & TILEMASK;
3972 
3973  if ((xm != TILEMASK) && (ym != TILEMASK) && (zm != TILEMASK))
3974  {
3975  const UT_VoxelTile<T> *tile =
3976  getTile(x >> TILEBITS, y >> TILEBITS, z >> TILEBITS);
3977 
3978  vz = tile->lerp(xm, ym, zm, fx, fy, fz);
3979  }
3980  else
3981  {
3982  // We cross tile boundaries but we remain within
3983  // the voxel grid. We can thus avoid any bounds checking
3984  // and use operator() rather than getValue.
3985 
3986  // Lerp x:x+1, y, z
3987  vx = UT_VoxelTile<T>::lerpValues((*this)(x, y, z),
3988  (*this)(x+1, y, z),
3989  fx);
3990  // Lerp x:x+1, y+1, z
3991  vx1= UT_VoxelTile<T>::lerpValues((*this)(x, y+1, z),
3992  (*this)(x+1, y+1, z),
3993  fx);
3994  // Lerp x:x+1, y:y+1, z
3995  vy = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);
3996 
3997  // Lerp x:x+1, y, z+1
3998  vx = UT_VoxelTile<T>::lerpValues((*this)(x, y, z+1),
3999  (*this)(x+1, y, z+1),
4000  fx);
4001  // Lerp x:x+1, y+1, z+1
4002  vx1= UT_VoxelTile<T>::lerpValues((*this)(x, y+1, z+1),
4003  (*this)(x+1, y+1, z+1),
4004  fx);
4005 
4006  // Lerp x:x+1, y:y+1, z+1
4007  vy1 = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);
4008 
4009  // Lerp x:x+1, y:y+1, z:z+1
4010  vz = UT_VoxelTile<T>::lerpValues(vy, vy1, fz);
4011  }
4012  }
4013  else
4014  {
4015  // Lerp x:x+1, y, z
4016  vx = UT_VoxelTile<T>::lerpValues(getValue(x, y, z),
4017  getValue(x+1, y, z),
4018  fx);
4019  // Lerp x:x+1, y+1, z
4020  vx1= UT_VoxelTile<T>::lerpValues(getValue(x, y+1, z),
4021  getValue(x+1, y+1, z),
4022  fx);
4023 
4024  // Lerp x:x+1, y:y+1, z
4025  vy = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);
4026 
4027  // Lerp x:x+1, y, z+1
4028  vx = UT_VoxelTile<T>::lerpValues(getValue(x, y, z+1),
4029  getValue(x+1, y, z+1),
4030  fx);
4031  // Lerp x:x+1, y+1, z+1
4032  vx1= UT_VoxelTile<T>::lerpValues(getValue(x, y+1, z+1),
4033  getValue(x+1, y+1, z+1),
4034  fx);
4035 
4036  // Lerp x:x+1, y:y+1, z+1
4037  vy1 = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);
4038 
4039  // Lerp x:x+1, y:y+1, z:z+1
4040  vz = UT_VoxelTile<T>::lerpValues(vy, vy1, fz);
4041  }
4042 
4043  return vz;
4044 }
4045 
4046 #if 0
4047 template <typename T>
4048 T
4050 {
4051  int x, y, z;
4052 
4053  v4ui idx;
4054 
4055  // We go from the position in the unit cube into the index.
4056  // The center of cells must map to the exact integer indices.
4057  pos *= v4uf((float) myRes[0], (float) myRes[1], (float) myRes[2], 0.0f);
4058 
4059  // Because we use truncation we'll get inaccurate results
4060  // at the zero boundary, we thus bump up by two
4061  pos += 1.5;
4062 
4063  idx = pos.splitFloat();
4064 
4065  // And counteract that bump
4066  idx -= 2;
4067 
4068  x = idx[0];
4069  y = idx[1];
4070  z = idx[2];
4071 
4072  // Do trilinear interpolation.
4073  // A | B
4074  // 0 +z +y +yz | +x +xz +xy +xyz
4075 
4076  v4uf a, b, fx, fy, fz;
4077 
4078  int xm, ym, zm;
4079 
4080  xm = x & TILEMASK;
4081  ym = y & TILEMASK;
4082  zm = z & TILEMASK;
4083 
4084  // Optimization for common case (values are within the voxel range and
4085  // are all within the same tile)
4086  if ( (x > 0) & (y > 0) & (z > 0) &
4087  (x < myRes[0]-1) & (y < myRes[1]-1) & (z < myRes[2]-1) &
4088  (xm != TILEMASK) & (ym != TILEMASK) & (zm != TILEMASK) )
4089 /*
4090  if (isValidIndex(x, y, z) && isValidIndex(x+1, y+1, z+1) &&
4091  (x & TILEMASK) != TILEMASK &&
4092  (y & TILEMASK) != TILEMASK &&
4093  (z & TILEMASK) != TILEMASK)
4094 */
4095  {
4096  const UT_VoxelTile<T> *tile =
4097  getTile(x >> TILEBITS, y >> TILEBITS, z >> TILEBITS);
4098 
4099  return tile->lerp(pos, xm, ym, zm);
4100  }
4101  else
4102  {
4103  a = v4uf( getValue(x, y, z),
4104  getValue(x, y, z+1),
4105  getValue(x, y+1, z),
4106  getValue(x, y+1, z+1) );
4107  b = v4uf( getValue(x+1, y, z),
4108  getValue(x+1, y, z+1),
4109  getValue(x+1, y+1, z),
4110  getValue(x+1, y+1, z+1) );
4111  }
4112 
4113  fx = pos.swizzle<0, 0, 0, 0>();
4114  fy = pos.swizzle<1, 1, 1, 1>();
4115  fz = pos.swizzle<2, 2, 2, 2>();
4116 
4117  b -= a;
4118  a = madd(b, fx, a);
4119 
4120  b = a.swizzle<2, 3, 0, 1>();
4121  b -= a;
4122  a = madd(b, fy, a);
4123 
4124  b = a.swizzle<1, 2, 3, 0>();
4125  b -= a;
4126  a = madd(b, fz, a);
4127 
4128  return a[0];
4129 }
4130 #endif
4131 
4132 static inline int
4133 firstTile(int &start, int &end, int res)
4134 {
4135  if (start < 0)
4136  {
4137  start += res;
4138  end += res;
4139  }
4140  return start;
4141 }
4142 
4143 static inline int
4144 nextTile(int &tile, int &pstart, int &start, int &end, int res)
4145 {
4146  int pend;
4147 
4148  if (pstart >= res)
4149  {
4150  pstart -= res;
4151  start -= res;
4152  end -= res;
4153  }
4154  tile = pstart >> TILEBITS;
4155  pend = SYSmin((tile+1) * TILESIZE, end, res);
4156 
4157  UT_ASSERT(pstart >= 0 && pstart < res);
4158  UT_ASSERT(pend > 0 && pstart <= res);
4159  UT_ASSERT(pend-pstart > 0);
4160  UT_ASSERT(pend-pstart <= TILESIZE);
4161 
4162  return pend;
4163 }
4164 
4165 template <typename T>
4166 T
4168  fpreal radius, int clampaxis) const
4169 {
4170  UT_Vector3 tpos;
4171  UT_FilterWindow win[3];
4172  UT_FilterWrap wrap[3];
4173  UT_VoxelTile<T> *tile;
4174  const float *weights[3];
4175  fpreal visible;
4176  int start[3], end[3], size[3];
4177  int pstart[3], pend[3];
4178  int tx, ty, tz, i;
4179  T result;
4180 
4181  if (!myTiles)
4182  {
4183  // If the array is empty, just use the border value.
4184  return myBorderValue;
4185  }
4186 
4187  UT_FilterWrap basewrap;
4188  switch (myBorderType)
4189  {
4190  case UT_VOXELBORDER_CONSTANT: basewrap = UT_WRAP_BORDER; break;
4191  case UT_VOXELBORDER_REPEAT: basewrap = UT_WRAP_REPEAT; break;
4192  case UT_VOXELBORDER_STREAK: basewrap = UT_WRAP_CLAMP; break;
4193 
4194  // NOTE: This is incorrect!
4195  case UT_VOXELBORDER_EXTRAP: basewrap = UT_WRAP_CLAMP; break;
4196  default: basewrap = UT_WRAP_CLAMP; break;
4197  }
4198  for (i = 0; i < 3; i++)
4199  {
4200  if (i == clampaxis)
4201  wrap[i] = UT_WRAP_CLAMP;
4202  else
4203  wrap[i] = basewrap;
4204  }
4205 
4206  memset(&result, 0, sizeof(result));
4207 
4208  radius *= filter.getSupport();
4209 
4210  // Make a local copy of the position so that we can modify it.
4211  tpos = pos;
4212  visible = 1;
4213  for (i = 0; i < 3; i++)
4214  {
4215  tpos[i] = tpos[i]*myRes[i];
4216  if (!win[i].setWeights(filter, tpos[i], radius, myRes[i], wrap[i]))
4217  return result;
4218 
4219  weights[i] = win[i].getWeights();
4220  start[i] = win[i].getStart() % myRes[i];
4221  size[i] = win[i].getSize();
4222 
4223  UT_ASSERT(start[i] >= 0);
4224  UT_ASSERT(size[i] <= myRes[i]);
4225 
4226  end[i] = start[i] + size[i];
4227  visible *= win[i].getVisible();
4228  }
4229 
4230  // Accumulate filtered results
4231  pstart[2] = firstTile(start[2], end[2], myRes[2]);
4232  while (pstart[2] < end[2])
4233  {
4234  pend[2] = nextTile(tz, pstart[2], start[2], end[2], myRes[2]);
4235  pstart[1] = firstTile(start[1], end[1], myRes[1]);
4236  while (pstart[1] < end[1])
4237  {
4238  pend[1] = nextTile(ty, pstart[1], start[1], end[1], myRes[1]);
4239  pstart[0] = firstTile(start[0], end[0], myRes[0]);
4240  while (pstart[0] < end[0])
4241  {
4242  pend[0] = nextTile(tx, pstart[0], start[0], end[0], myRes[0]);
4243  tile = getTile(tx, ty, tz);
4244  UT_ASSERT(tile);
4245  tile->weightedSum(pstart, pend, weights, start, result);
4246  pstart[0] = pend[0];
4247  }
4248  pstart[1] = pend[1];
4249  }
4250  pstart[2] = pend[2];
4251  }
4252 
4253  if (visible < 1)
4254  {
4255  result += (1-visible)*myBorderValue;
4256  }
4257 
4258  return result;
4259 }
4260 
4261 template <typename T>
4262 void
4264  UT_FilterType filtertype,
4265  float filterwidthscale,
4266  int clampaxis)
4267 {
4268  fpreal radius;
4269  UT_Filter *filter;
4270 
4271  filter = UT_Filter::getFilter(filtertype);
4272 
4273  radius = SYSmax( src.getXRes() / (fpreal)getXRes(),
4274  src.getYRes() / (fpreal)getYRes(),
4275  src.getZRes() / (fpreal)getZRes(),
4276  1.0f );
4277  radius *= 0.5 * filterwidthscale;
4278 
4279  resamplethread(src, filter, radius, clampaxis);
4280 
4281  UT_Filter::releaseFilter(filter);
4282 }
4283 
4284 template <typename T>
4285 void
4286 UT_VoxelArray<T>::flattenPartial(T *flatarray, exint ystride, exint zstride,
4287  const UT_JobInfo &info) const
4288 {
4290 
4291  // Const cast.
4292  vit.setArray((UT_VoxelArray<T> *)this);
4293  vit.splitByTile(info);
4294 
4295  for (vit.rewind(); !vit.atEnd(); vit.advance())
4296  {
4297  flatarray[vit.x() + vit.y()*ystride + vit.z()*zstride] = vit.getValue();
4298  }
4299 }
4300 
4301 
4302 
4303 template <>
4304 inline void
4306  exint ystride, exint zstride,
4307  UT_Vector4F,
4308  const UT_JobInfo &info) const
4309 {
4311 
4312  // Const cast.
4313  vit.setArray(SYSconst_cast(this));
4314  vit.splitByTile(info);
4315 
4316  for (vit.rewind(); !vit.atEnd(); vit.advance())
4317  {
4318  UT_Vector4F v = vit.getValue();
4319  int idx = (vit.x() + vit.y()*ystride + vit.z()*zstride) * 4;
4320 
4321  flatarray[idx] = (uint8) SYSclamp(v.x() * 255.0f, 0.0f, 255.0f);
4322  flatarray[idx+1] = (uint8) SYSclamp(v.y() * 255.0f, 0.0f, 255.0f);
4323  flatarray[idx+2] = (uint8) SYSclamp(v.z() * 255.0f, 0.0f, 255.0f);
4324  flatarray[idx+3] = (uint8) SYSclamp(v.w() * 255.0f, 0.0f, 255.0f);
4325  }
4326 }
4327 
4328 template <typename T>
4329 void
4331  exint ystride, exint zstride,
4332  T, const UT_JobInfo &info) const
4333 {
4334  UT_ASSERT(!"This template requires specific instantiations.");
4335 }
4336 
4337 template <>
4338 inline void
4340  exint ystride, exint zstride,
4341  UT_Vector4F,
4342  const UT_JobInfo &info) const
4343 {
4345 
4346  // Const cast.
4347  vit.setArray(SYSconst_cast(this));
4348  vit.splitByTile(info);
4349 
4350  for (vit.rewind(); !vit.atEnd(); vit.advance())
4351  flatarray[vit.x() + vit.y()*ystride + vit.z()*zstride] = vit.getValue();
4352 }
4353 
4354 template <typename T>
4355 void
4357  exint ystride, exint zstride,
4358  T, const UT_JobInfo &info) const
4359 {
4360  UT_ASSERT(!"This template requires specific instantiations.");
4361 }
4362 
4363 template <>
4364 inline void
4366  exint ystride, exint zstride,
4367  UT_Vector4F,
4368  const UT_JobInfo &info) const
4369 {
4371 
4372  // Const cast.
4373  vit.setArray(SYSconst_cast(this));
4374  vit.splitByTile(info);
4375 
4376  for (vit.rewind(); !vit.atEnd(); vit.advance())
4377  {
4378  UT_Vector4F v = vit.getValue();
4379 
4380  // NOTE: This works around an Nvidia driver bug on OSX, and older
4381  // Nvidia drivers on other platforms. The problem was that very
4382  // small FP values were causing huge random values to appear in
4383  // the 3D Texture.
4384  if(SYSabs(v.x()) < 1e-9)
4385  v.x() = 0.0;
4386  if(SYSabs(v.y()) < 1e-9)
4387  v.y() = 0.0;
4388  if(SYSabs(v.z()) < 1e-9)
4389  v.z() = 0.0;
4390  if(SYSabs(v.w()) < 1e-9)
4391  v.w() = 0.0;
4392 
4393  flatarray[vit.x() + vit.y()*ystride + vit.z()*zstride] = v;
4394  }
4395 }
4396 
4397 template <typename T>
4398 void
4400  exint ystride, exint zstride,
4401  T,const UT_JobInfo &info) const
4402 {
4403  UT_ASSERT(!"This template requires specific instantiations.");
4404 }
4405 
4406 
4407 template <typename T>
4408 void
4410  exint ystride, exint zstride,
4411  const UT_JobInfo &info)
4412 {
4414 
4415  vit.setArray(this);
4416  vit.splitByTile(info);
4417  vit.setCompressOnExit(true);
4418 
4419  for (vit.rewind(); !vit.atEnd(); vit.advance())
4420  {
4421  vit.setValue(flatarray[vit.x() + vit.y()*ystride + vit.z()*zstride]);
4422  }
4423 }
4424 
4425 template <typename T>
4426 template <typename S>
4427 S *
4429  const UT_IntArray &tilelist) const
4430 {
4431  for (int i = 0; i < tilelist.entries(); i++)
4432  {
4433  UT_ASSERT(tilelist(i) >= 0 && tilelist(i) < numTiles());
4434  const UT_VoxelTile<T> *tile = getLinearTile(tilelist(i));
4435 
4436  tile->flatten(dstdata, stride);
4437  dstdata += tile->numVoxels() * stride;
4438  }
4439  return dstdata;
4440 }
4441 
4442 template <typename T>
4443 template <typename S>
4444 const S *
4446  const UT_IntArray &tilelist)
4447 {
4448  for (int i = 0; i < tilelist.entries(); i++)
4449  {
4450  UT_ASSERT(tilelist(i) >= 0 && tilelist(i) < numTiles());
4451  UT_VoxelTile<T> *tile = getLinearTile(tilelist(i));
4452 
4453  tile->writeData(srcdata, stride);
4454  srcdata += tile->numVoxels() * stride;
4455  }
4456  return srcdata;
4457 }
4458 
4459 template <typename T>
4460 void
4462  int offx, int offy, int offz)
4463 {
4464  UT_VoxelBorderType srcborder;
4465  T srcborderval;
4466 
4467  srcborder = src.getBorder();
4468  srcborderval = src.getBorderValue();
4469 
4470  if (srcborder != UT_VOXELBORDER_EXTRAP)
4471  {
4472  // These borders may be constant friendly
4473  T srcval;
4474 
4475  if (src.isConstant(&srcval))
4476  {
4477  if (srcborder != UT_VOXELBORDER_CONSTANT ||
4478  SYSisEqual(srcborderval, srcval))
4479  {
4480  // We can trivially make ourselves constant.
4481  constant(srcval);
4482  return;
4483  }
4484  }
4485  }
4486 
4487  copyWithOffsetInternal(src, offx, offy, offz);
4488 }
4489 
4490 template <typename T>
4491 void
4493  int offx, int offy, int offz,
4494  const UT_JobInfo &info)
4495 {
4497  UT_Vector3I off(offx, offy, offz);
4498 
4499  vit.setArray(this);
4500  vit.splitByTile(info);
4501  vit.setCompressOnExit(true);
4502 
4503  // Iterate over all tiles
4504  for (vit.rewind(); !vit.atEnd(); vit.advanceTile())
4505  {
4506  // Check out this tile
4507  UT_VoxelTile<T> *tile = vit.getTile();
4508  UT_VoxelTile<T> *srctile;
4509  int tx, ty, tz;
4510 
4511  UT_Vector3I start, end, srctileidx, srctileoff, srcend, srcendtile;
4512  vit.getTileVoxels(start, end);
4513 
4514  srctileidx = start;
4515  srctileidx += off;
4516  srctileidx.x() >>= TILEBITS;
4517  srctileidx.y() >>= TILEBITS;
4518  srctileidx.z() >>= TILEBITS;
4519  srctileoff.x() = off.x() & TILEMASK;
4520  srctileoff.y() = off.y() & TILEMASK;
4521  srctileoff.z() = off.z() & TILEMASK;
4522 
4523  srcend = start;
4524  // We are very careful here to be inclusive, so we don't trigger
4525  // the next tile. This is the largest index we expect to index.
4526  srcend.x() += tile->xres() - 1;
4527  srcend.y() += tile->yres() - 1;
4528  srcend.z() += tile->zres() - 1;
4529  srcend += off;
4530  srcendtile = srcend;
4531  srcendtile.x() >>= TILEBITS;
4532  srcendtile.y() >>= TILEBITS;
4533  srcendtile.z() >>= TILEBITS;
4534 
4535  UT_ASSERT(srcendtile.x() == srctileidx.x() ||
4536  srcendtile.x() == srctileidx.x()+1);
4537  UT_ASSERT(srcendtile.y() == srctileidx.y() ||
4538  srcendtile.y() == srctileidx.y()+1);
4539  UT_ASSERT(srcendtile.z() == srctileidx.z() ||
4540  srcendtile.z() == srctileidx.z()+1);
4541 
4542  // Check if we are fully in bounds.
4543  if (srctileidx.x() >= 0 &&
4544  srctileidx.y() >= 0 &&
4545  srctileidx.z() >= 0 &&
4546  srcendtile.x() < src.getTileRes(0) &&
4547  srcendtile.y() < src.getTileRes(1) &&
4548  srcendtile.z() < src.getTileRes(2) &&
4549  srcend.x() < src.getXRes() &&
4550  srcend.y() < src.getYRes() &&
4551  srcend.z() < src.getZRes())
4552  {
4553  bool allconst = true, firsttile = true;
4554  T constval, cval;
4555  // Check if we are all constant...
4556  for (tz = srctileidx.z(); allconst && tz <= srcendtile.z(); tz++)
4557  for (ty = srctileidx.y(); allconst && ty <= srcendtile.y(); ty++)
4558  for (tx = srctileidx.x(); tx <= srcendtile.x(); tx++)
4559  {
4560  srctile = src.getTile(tx, ty, tz);
4561  if (!srctile->isConstant())
4562  {
4563  allconst = false;
4564  break;
4565  }
4566  cval = (*srctile)(0, 0, 0);
4567  if (firsttile)
4568  {
4569  firsttile = false;
4570  constval = cval;
4571  }
4572  if (!SYSisEqual(cval, constval))
4573  {
4574  allconst = false;
4575  break;
4576  }
4577  }
4578  if (allconst)
4579  {
4580  // Should have found at least one tile!
4581  UT_ASSERT(!firsttile);
4582  tile->makeConstant(constval);
4583 
4584  // Onto the next tile!
4585  continue;
4586  }
4587  }
4588 
4589  // All of the fragments aren't constant, or aren't all inside
4590  // our range. So we have to work fragment at a time..
4591  for (tz = srctileidx.z(); tz <= srcendtile.z(); tz++)
4592  for (ty = srctileidx.y(); ty <= srcendtile.y(); ty++)
4593  for (tx = srctileidx.x(); tx <= srcendtile.x(); tx++)
4594  {
4595  int destx, desty, destz;
4596  int srcx, srcy, srcz;
4597 
4598  destx = (tx == srctileidx.x()) ? 0 : (TILESIZE-srctileoff.x());
4599  desty = (ty == srctileidx.y()) ? 0 : (TILESIZE-srctileoff.y());
4600  destz = (tz == srctileidx.z()) ? 0 : (TILESIZE-srctileoff.z());
4601  srcx = (tx == srctileidx.x()) ? srctileoff.x() : 0;
4602  srcy = (ty == srctileidx.y()) ? srctileoff.y() : 0;
4603  srcz = (tz == srctileidx.z()) ? srctileoff.z() : 0;
4604 #if 1
4605  if (tx >= 0 &&
4606  ty >= 0 &&
4607  tz >= 0 &&
4608  tx < src.getTileRes(0) &&
4609  ty < src.getTileRes(1) &&
4610  tz < src.getTileRes(2) &&
4611  ((tx != src.getTileRes(0)-1) ||
4612  srcend.x() < src.getXRes()) &&
4613  ((ty != src.getTileRes(1)-1) ||
4614  srcend.y() < src.getYRes()) &&
4615  ((tz != src.getTileRes(2)-1) ||
4616  srcend.z() < src.getZRes())
4617  )
4618  {
4619  srctile = src.getTile(tx, ty, tz);
4620  // In bounds
4621  tile->copyFragment(
4622  destx, desty, destz,
4623  *srctile,
4624  srcx, srcy, srcz
4625  );
4626  }
4627  else
4628 #endif
4629  {
4630  // Out of bounds!
4631  int maxd = SYSmin(tile->zres(),
4632  TILESIZE-srcz);
4633  int maxh = SYSmin(tile->yres(),
4634  TILESIZE-srcy);
4635  int maxw = SYSmin(tile->xres(),
4636  TILESIZE-srcx);
4637  for (int z = destz; z < maxd; z++)
4638  {
4639  for (int y = desty; y < maxh; y++)
4640  {
4641  for (int x = destx; x < maxw; x++)
4642  {
4643  T val;
4644 
4645  val = src.getValue(x + vit.x() + offx,
4646  y + vit.y() + offy,
4647  z + vit.z() + offz);
4648  tile->setValue(x, y, z, val);
4649  }
4650  }
4651  }
4652  }
4653  }
4654  }
4655 }
4656 
4657 template <typename T>
4658 void
4660  const UT_Filter *filter, float radius,
4661  int clampaxis,
4662  const UT_JobInfo &info)
4663 {
4665  UT_Vector3 pos;
4666  UT_Vector3 ratio;
4667  UT_Interrupt *boss = UTgetInterrupt();
4668 
4669  vit.setArray(this);
4670  vit.splitByTile(info);
4671  vit.setCompressOnExit(true);
4672  vit.setInterrupt(boss);
4673 
4674  ratio.x() = 1.0f / getXRes();
4675  ratio.y() = 1.0f / getYRes();
4676  ratio.z() = 1.0f / getZRes();
4677 
4678  for (vit.rewind(); !vit.atEnd(); vit.advance())
4679  {
4680  pos.x() = vit.x()+0.5f;
4681  pos.y() = vit.y()+0.5f;
4682  pos.z() = vit.z()+0.5f;
4683  pos *= ratio;
4684 
4685  vit.setValue(src.evaluate(pos, *filter, radius, clampaxis));
4686  }
4687 }
4688 
4689 template <typename T>
4690 bool
4691 UT_VoxelArray<T>::posToIndex(UT_Vector3 pos, int &x, int &y, int &z) const
4692 {
4693  // We go from the position in the unit cube into the index.
4694  // The center of cells must map to the exact integer indices.
4695  pos.x() *= myRes[0];
4696  pos.y() *= myRes[1];
4697  pos.z() *= myRes[2];
4698 
4699  // The centers of cells are now mapped .5 too high. Ie, the true
4700  // center of cell index (0,0,0) would be (0.5,0.5,0.5) This, however,
4701  // is exactly what we want for rounding.
4702  x = (int) SYSfloor(pos.x());
4703  y = (int) SYSfloor(pos.y());
4704  z = (int) SYSfloor(pos.z());
4705 
4706  // Determine if out of bounds.
4707  return isValidIndex(x, y, z);
4708 }
4709 
4710 template <typename T>
4711 bool
4713 {
4714  // We go from the position in the unit cube into the index.
4715  // The center of cells must map to the exact integer indices.
4716  pos.x() *= myRes[0];
4717  pos.y() *= myRes[1];
4718  pos.z() *= myRes[2];
4719 
4720  // The centers of cells are now mapped .5 too high. Ie, the true
4721  // center of cell index (0,0,0) would be (0.5,0.5,0.5)
4722  pos.x() -= 0.5;
4723  pos.y() -= 0.5;
4724  pos.z() -= 0.5;
4725 
4726  ipos = pos;
4727 
4728  // Determine if out of bounds.
4729  if (pos.x() < 0 || pos.x() >= myRes[0] ||
4730  pos.y() < 0 || pos.y() >= myRes[1] ||
4731  pos.z() < 0 || pos.z() >= myRes[2])
4732  return false;
4733 
4734  return true;
4735 }
4736 
4737 template <typename T>
4738 bool
4739 UT_VoxelArray<T>::indexToPos(int x, int y, int z, UT_Vector3F &pos) const
4740 {
4741  pos.x() = x;
4742  pos.y() = y;
4743  pos.z() = z;
4744 
4745  // Move the indices to the centers of the cells.
4746  pos.x() += 0.5;
4747  pos.y() += 0.5;
4748  pos.z() += 0.5;
4749 
4750  // And scale into the unit cube.
4751  pos *= myInvRes;
4752 
4753  // Return true if the original coordinates were in range.
4754  return isValidIndex(x, y, z);
4755 }
4756 
4757 template <typename T>
4758 bool
4759 UT_VoxelArray<T>::indexToPos(int x, int y, int z, UT_Vector3D &pos) const
4760 {
4761  pos.x() = x;
4762  pos.y() = y;
4763  pos.z() = z;
4764 
4765  // Move the indices to the centers of the cells.
4766  pos.x() += 0.5;
4767  pos.y() += 0.5;
4768  pos.z() += 0.5;
4769 
4770  // And scale into the unit cube.
4771  pos *= myInvRes;
4772 
4773  // Return true if the original coordinates were in range.
4774  return isValidIndex(x, y, z);
4775 }
4776 
4777 template <typename T>
4778 void
4780 {
4781  pos = index;
4782 
4783  // Move the indices to the centers of the cells.
4784  pos.x() += 0.5;
4785  pos.y() += 0.5;
4786  pos.z() += 0.5;
4787 
4788  // And scale into the unit cube.
4789  pos *= myInvRes;
4790 }
4791 
4792 template <typename T>
4793 void
4795 {
4796  pos = index;
4797 
4798  // Move the indices to the centers of the cells.
4799  pos.x() += 0.5;
4800  pos.y() += 0.5;
4801  pos.z() += 0.5;
4802 
4803  // And scale into the unit cube.
4804  pos *= myInvRes;
4805 }
4806 
4807 template <typename T>
4808 void
4810 {
4811  myBorderType = type;
4812  myBorderValue = t;
4813 }
4814 
4815 template <typename T>
4816 void
4818 {
4819  myBorderScale[0] = sx;
4820  myBorderScale[1] = sy;
4821  myBorderScale[2] = sz;
4822 }
4823 
4824 template <typename T>
4825 void
4827 {
4828  UT_VoxelArrayIterator<T> vit(this);
4829  vit.splitByTile(info);
4830  for (vit.rewind(); !vit.atEnd(); vit.advanceTile())
4831  {
4832  int i = vit.getLinearTileNum();
4833  myTiles[i].tryCompress(getCompressionOptions());
4834  }
4835 }
4836 
4837 template <typename T>
4838 void
4840 {
4841  UT_VoxelArrayIterator<T> vit(this);
4842  vit.splitByTile(info);
4843  for (vit.rewind(); !vit.atEnd(); vit.advanceTile())
4844  {
4845  int i = vit.getLinearTileNum();
4846  myTiles[i].uncompress();
4847  }
4848 }
4849 
4850 template <typename T>
4851 void
4853 {
4854  UT_VoxelArrayIterator<T> vit(this);
4855  vit.splitByTile(info);
4856  for (vit.rewind(); !vit.atEnd(); vit.advanceTile())
4857  {
4858  int i = vit.getLinearTileNum();
4859  if (!myTiles[i].isConstant())
4860  myTiles[i].uncompress();
4861  }
4862 }
4863 
4864 template <typename T>
4865 void
4866 UT_VoxelArray<T>::saveData(std::ostream &os) const
4867 {
4868  T cval;
4869  char version;
4870 
4871  // First determine if we are fully constant.
4872  if (isConstant(&cval))
4873  {
4874  // Save a constant array.
4875  version = 0;
4876  UTwrite(os, &version, 1);
4877  UTwrite<T>(os, &cval);
4878  return;
4879  }
4880 
4881  // Compressed tiles.
4882  version = 1;
4883  UTwrite(os, &version, 1);
4884 
4885  // Store list of compression types
4887 
4888  int i, ntiles;
4889 
4890  ntiles = numTiles();
4891  for (i = 0; i < ntiles; i++)
4892  {
4893  myTiles[i].save(os);
4894  }
4895 }
4896 
4897 
4898 template <typename T>
4899 void
4901 {
4902  T cval;
4903  char version;
4904 
4905  is.readChar(version);
4906 
4907  // First determine if we are fully constant.
4908  if (version == 0)
4909  {
4910  // Save a constant array.
4911  is.read<T>(&cval);
4912 
4913  constant(cval);
4914  return;
4915  }
4916 
4917  if (version == 1)
4918  {
4919  UT_IntArray compressions;
4920 
4921  // Store list of compression types
4922  UT_VoxelTile<T>::loadCompressionTypes(is, compressions);
4923 
4924  int i, ntiles;
4925 
4926  ntiles = numTiles();
4927  for (i = 0; i < ntiles; i++)
4928  {
4929  myTiles[i].load(is, compressions);
4930  }
4931  }
4932 }
4933 
4934 template <typename T>
4935 bool
4936 UT_VoxelArray<T>::saveData(UT_JSONWriter &w, const char *shared_mem_owner) const
4937 {
4938  bool ok = true;
4939  T cval;
4940  int8 version;
4941 
4942  // First determine if we are fully constant.
4943  if (isConstant(&cval))
4944  {
4945  ok = ok && w.jsonBeginArray();
4946  // Save a constant array.
4947  ok = ok && w.jsonKey(UT_VoxelArrayJSON::getToken(
4949  ok = ok && w.jsonValue(cval);
4950  ok = ok && w.jsonEndArray();
4951  return ok;
4952  }
4953 
4954  if (shared_mem_owner)
4955  {
4956  SYS_SharedMemory *shm;
4957  shm = copyToSharedMemory(shared_mem_owner);
4958  if (shm)
4959  {
4960  ok = ok && w.jsonBeginArray();
4961  // Save a shared memory array.
4962  ok = ok && w.jsonKey(UT_VoxelArrayJSON::getToken(
4964  ok = ok && w.jsonString(shm->id());
4965  ok = ok && w.jsonEndArray();
4966  return ok;
4967  }
4968  // Fall back on raw write.
4969  }
4970 
4971 
4972  // Compressed tiles.
4973  ok = ok && w.jsonBeginArray();
4974  ok = ok && w.jsonKey(UT_VoxelArrayJSON::getToken(
4976  ok = ok && w.jsonBeginArray();
4977  ok = ok && w.jsonKey(UT_VoxelArrayJSON::getToken(
4979  version = 1;
4980  ok = ok && w.jsonValue(version);
4981 
4982  // Store list of compression types
4983  ok = ok && w.jsonKey(UT_VoxelArrayJSON::getToken(
4986 
4987  ok = ok && w.jsonKey(UT_VoxelArrayJSON::getToken(
4989  ok = ok && w.jsonBeginArray();
4990  int i, ntiles;
4991 
4992  ntiles = numTiles();
4993  for (i = 0; i < ntiles; i++)
4994  {
4995  ok = ok && myTiles[i].save(w);
4996  }
4997  ok = ok && w.jsonEndArray();
4998 
4999  ok = ok && w.jsonEndArray();
5000 
5001  ok = ok && w.jsonEndArray();
5002  return ok;
5003 }
5004 
5005 template <typename T>
5006 bool
5008 {
5009  T cval;
5010  int8 version;
5011  UT_WorkBuffer key;
5012  int key_id;
5013  bool array_error = false;
5014 
5015  delete mySharedMemView;
5016  mySharedMemView = 0;
5017 
5018  delete mySharedMem;
5019  mySharedMem = 0;
5020 
5021  if (!p.parseBeginArray(array_error) || array_error)
5022  return false;
5023 
5024  if (!p.parseString(key))
5025  return false;
5026  key_id = UT_VoxelArrayJSON::getArrayID(key.buffer());
5028  {
5029  if (!p.parseNumber(cval))
5030  return false;
5031 
5032  constant(cval);
5033  }
5034  else if (key_id == UT_VoxelArrayJSON::ARRAY_SHAREDARRAY)
5035  {
5036  UT_WorkBuffer shm_id;
5037  if (!p.parseString(shm_id))
5038  return false;
5039 
5040  if (!populateFromSharedMemory(shm_id.buffer()))
5041  {
5042  // If the shared memory disappears before we get a chance to read
5043  // it, don't abort, just set a constant value. That way Mantra has
5044  // a chance to recover. Not the most elegant solution but works for
5045  // now.
5046  T cval;
5047  cval = 0;
5048  constant(cval);
5049  }
5050  }
5051  else if (key_id == UT_VoxelArrayJSON::ARRAY_TILEDARRAY)
5052  {
5053  UT_JSONParser::traverser it, tile_it;
5054  UT_IntArray compressions;
5055  int i, ntiles;
5056 
5057 
5058  for (it = p.beginArray(); !it.atEnd(); ++it)
5059  {
5060  if (!it.getLowerKey(key))
5061  return false;
5062  switch (UT_VoxelArrayJSON::getArrayID(key.buffer()))
5063  {
5065  if (!p.parseNumber(version))
5066  return false;
5067  break;
5069  // Store list of compression types
5070  if (!UT_VoxelTile<T>::loadCompressionTypes(p, compressions))
5071  return false;
5072  break;
5073 
5075  ntiles = numTiles();
5076  for (tile_it = p.beginArray(), i = 0; !tile_it.atEnd();
5077  ++tile_it, ++i)
5078  {
5079  if (i < ntiles)
5080  {
5081  if (!myTiles[i].load(p, compressions))
5082  return false;
5083  }
5084  else
5085  {
5086  UT_VoxelTile<T> dummy_tile;
5087  dummy_tile.setRes(TILESIZE, TILESIZE, TILESIZE);
5088  if (!dummy_tile.load(p, compressions))
5089  return false;
5090  }
5091  }
5092  if (i != ntiles)
5093  return false;
5094  break;
5095  default:
5096  p.addWarning("Unexpected key for voxel data: %s",
5097  key.buffer());
5098  if (!p.skipNextObject())
5099  return false;
5100  }
5101  }
5102  }
5103  else
5104  {
5105  p.addWarning("Unexpected voxel key: %s\n", key.buffer());
5106  p.skipNextObject();
5107  return false;
5108  }
5109 
5110  if (!p.parseEndArray(array_error) || array_error)
5111  return false;
5112 
5113  return true;
5114 }
5115 
5116 template<typename T>
5118 UT_VoxelArray<T>::copyToSharedMemory(const char *shared_mem_owner) const
5119 {
5121 
5122  int ntiles;
5123 
5124  ntiles = numTiles();
5125 
5126  // Tally up the size we need.
5127  exint total_mem_size;
5128 
5129  total_mem_size = sizeof(exint); // Tile count.
5130 
5131  UT_Array<exint> tile_sizes;
5132  UT_Array<exint> next_block_offsets;
5133 
5134  for (int i = 0; i < ntiles; i++)
5135  {
5136  exint tile_size, block_size;
5137  if (myTiles[i].isConstant())
5138  tile_size = sizeof(T);
5139  else
5140  tile_size = myTiles[i].getDataLength();
5141 
5142  tile_sizes.append(tile_size);
5143 
5144  block_size = SYSroundUpToMultipleOf<exint>(tile_size, sizeof(exint));
5145  block_size += sizeof(exint) * 2; // Offset + compression type
5146 
5147  if (i < (ntiles-1))
5148  next_block_offsets.append(block_size / sizeof(exint));
5149  else
5150  next_block_offsets.append(0);
5151 
5152  total_mem_size += block_size;
5153  }
5154 
5155  // Start with an empty shared memory.
5156  SYS_SharedMemory *shared_mem;
5157  UT_WorkBuffer shared_mem_key;
5158 
5159  shared_mem_key.sprintf("%s:%p", shared_mem_owner, this);
5160 
5161  shared_mem = shmgr.get(shared_mem_key.buffer());
5162 
5163  // Does the existing block have the same number of tiles and same
5164  // offsets? In that case we don't reset the shared memory, just write
5165  // the tiles out again. There's very little adverse effect in changing
5166  // the voxel tile data.
5167  bool same_meta_data = false;
5168  if (shared_mem->size() >= sizeof(exint))
5169  {
5170  same_meta_data = true;
5171  {
5172  SYS_SharedMemoryView sv_tc(*shared_mem, 0, sizeof(exint));
5173 
5174  exint *ptr_ds = (exint *)sv_tc.data();
5175  if (*ptr_ds != ntiles)
5176  same_meta_data = false;
5177  }
5178 
5179  if (same_meta_data)
5180  {
5181  // Check if the offsets + compression type is the same.
5182  exint offset = sizeof(exint);
5183  for (int i = 0; i < ntiles; i++)
5184  {
5185  SYS_SharedMemoryView sv_tile(*shared_mem, offset,
5186  sizeof(exint) * 2);
5187 
5188  exint *ptr_tile = (exint *)sv_tile.data();
5189  if (ptr_tile[0] != next_block_offsets(i) ||
5190  ptr_tile[1] != (exint)myTiles[i].myCompressionType)
5191  {
5192  same_meta_data = false;
5193  break;
5194  }
5195  offset += next_block_offsets(i) * sizeof(exint);
5196  }
5197  }
5198  }
5199 
5200  if (!same_meta_data)
5201  {
5202  if (!shared_mem->reset(total_mem_size) ||
5203  shared_mem->size() != total_mem_size)
5204  {
5205  return NULL;
5206  }
5207  }
5208 
5209  {
5210  SYS_SharedMemoryView sv_tc(*shared_mem, 0, sizeof(exint));
5211 
5212  exint *ptr_ds = (exint *)sv_tc.data();
5213  *ptr_ds = ntiles;
5214  }
5215 
5216  exint offset = sizeof(exint);
5217  for (int i = 0; i < ntiles; i++)
5218  {
5219  SYS_SharedMemoryView sv_tile(*shared_mem, offset,
5220  sizeof(exint) * 2 + tile_sizes(i));
5221 
5222  exint *ptr_tile = (exint *)sv_tile.data();
5223 
5224  // Offset to next tile in exint units.
5225  ptr_tile[0] = next_block_offsets(i);
5226  ptr_tile[1] = myTiles[i].myCompressionType;
5227 
5228  ::memcpy(&ptr_tile[2], myTiles[i].rawData(), tile_sizes(i));
5229 
5230  offset += ptr_tile[0] * sizeof(exint);
5231  }
5232 
5233  return shared_mem;
5234 }
5235 
5236 template<typename T>
5237 bool
5239 {
5240  mySharedMem = new SYS_SharedMemory(id, /*read_only=*/true);
5241  if (!mySharedMem->size())
5242  {
5243  delete mySharedMem;
5244  mySharedMem = 0;
5245  return false;
5246  }
5247 
5248  exint ntiles;
5249  {
5250  SYS_SharedMemoryView sv_tc(*mySharedMem, 0, sizeof(exint));
5251  ntiles = *(const exint *)sv_tc.data();
5252  }
5253  if (ntiles != numTiles())
5254  {
5255  delete mySharedMem;
5256  mySharedMem = 0;
5257  return false;
5258  }
5259 
5260  mySharedMemView = new SYS_SharedMemoryView(*mySharedMem, sizeof(exint));
5261 
5262  exint *data = (exint *)mySharedMemView->data();
5263 
5264  for (int i = 0; i < ntiles; i++)
5265  {
5266  exint offset = data[0];
5267  int8 ctype = (int8)data[1];
5268 
5269  myTiles[i].setForeignData(&data[2], ctype);
5270 
5271  data += offset;
5272  }
5273 
5274  return true;
5275 }
5276 
5277 
5278 
5279 
5280 //
5281 // UT_VoxelMipMap functions
5282 //
5283 template <typename T>
5285 {
5286  initializePrivate();
5287 }
5288 
5289 template <typename T>
5291 {
5292  destroyPrivate();
5293 }
5294 
5295 template <typename T>
5297 {
5298  initializePrivate();
5299 
5300  *this = src;
5301 }
5302 
5303 template <typename T>
5304 const UT_VoxelMipMap<T> &
5306 {
5308  int level, i;
5309 
5310  if (&src == this)
5311  return *this;
5312 
5313  destroyPrivate();
5314 
5315  // We do a deep copy. We have learned from UT_String that
5316  // shallow copies will destroy us in the end.
5317  myOwnBase = true;
5318  myBaseLevel = new UT_VoxelArray<T>;
5319  *myBaseLevel = *src.myBaseLevel;
5320 
5321  myNumLevels = src.myNumLevels;
5322 
5323  for (i = 0; i < src.myLevels.entries(); i++)
5324  {
5325  levels = new UT_VoxelArray<T> *[myNumLevels];
5326  myLevels.append(levels);
5327 
5328  for (level = 0; level < myNumLevels; level++)
5329  {
5330  levels[level] = new UT_VoxelArray<T>;
5331  *levels[level] = *src.myLevels(i)[level];
5332  }
5333  }
5334  return *this;
5335 }
5336 
5337 template <typename T>
5338 void
5340  mipmaptype function)
5341 {
5342  UT_Array<mipmaptype> functions;
5343 
5344  functions.append(function);
5345  build(baselevel, functions);
5346 }
5347 
5348 template <typename T>
5349 void
5351  const UT_Array<mipmaptype> &functions)
5352 {
5353  destroyPrivate();
5354 
5355  // Setup our baselevel.
5356  myBaseLevel = baselevel;
5357  myOwnBase = false;
5358 
5359  // Now build from bottom up. First, calculate the number
5360  // of levels we'll need.
5361  myNumLevels = 0;
5362  int maxres, level;
5363 
5364  maxres = SYSmax(myBaseLevel->getXRes(), myBaseLevel->getYRes());
5365  maxres = SYSmax(maxres, myBaseLevel->getZRes());
5366 
5367  //cerr << "Max res " << maxres << endl;
5368 
5369  while (maxres > 1)
5370  {
5371  myNumLevels++;
5372  maxres++;
5373  maxres >>= 1;
5374  }
5375 
5376  //cerr << "Levels: " << myNumLevels << endl;
5377 
5378  // Create our desired voxel array list.
5379  for (int i = 0; i < functions.entries(); i++)
5380  {
5381  mipmaptype function = functions(i);
5382 
5383  UT_VoxelArray<T> **levels = new UT_VoxelArray<T> *[myNumLevels];
5384  myLevels.append(levels);
5385 
5386  UT_VoxelArray<T> *lastlevel = myBaseLevel;
5387  for (level = myNumLevels-1; level >= 0; level--)
5388  {
5389  // Find our resolutions.
5390  int xres = (lastlevel->getXRes() + 1) >> 1;
5391  int yres = (lastlevel->getYRes() + 1) >> 1;
5392  int zres = (lastlevel->getZRes() + 1) >> 1;
5393 
5394  levels[level] = new UT_VoxelArray<T>;
5395  levels[level]->size(xres, yres, zres);
5396 
5397  downsample(*levels[level], *lastlevel, function);
5398 
5399  lastlevel = levels[level];
5400  }
5401  }
5402 }
5403 
5404 template <typename T>
5405 void
5408  const UT_VoxelArray<T> &src,
5409  mipmaptype function,
5410  const UT_JobInfo &info)
5411 {
5412  UT_VoxelArrayIterator<T> vit(&dst);
5413 
5414  vit.splitByTile(info);
5415  vit.setCompressOnExit(true);
5416 
5417  for (vit.rewind(); !vit.atEnd(); vit.advance())
5418  {
5419  int x = 2*vit.x();
5420  int y = 2*vit.y();
5421  int z = 2*vit.z();
5422  T sam[8];
5423 
5424  if (src.extractSample(x, y, z, sam))
5425  {
5426  vit.setValue(sam[0]);
5427  }
5428  else
5429  {
5430  vit.setValue(
5431  mixValues(
5432  mixValues(
5433  mixValues(sam[0], sam[1], function),
5434  mixValues(sam[2], sam[3], function), function),
5435  mixValues(
5436  mixValues(sam[4], sam[5], function),
5437  mixValues(sam[6], sam[7], function), function),
5438  function));
5439  }
5440  }
5441 }
5442 
5443 
5444 template <typename T>
5445 int64
5447 {
5448  int64 mem = inclusive ? sizeof(*this) : 0;
5449 
5450  mem += myLevels.getMemoryUsage(false);
5451  mem += myLevels.entries() * myNumLevels * sizeof(UT_VoxelArray<T>*);
5452  for (exint j = 0; j < myLevels.entries(); j++)
5453  {
5454  for (int i = 0; i < myNumLevels; i++)
5455  mem += myLevels(j)[i]->getMemoryUsage(true);
5456  }
5457 
5458  return mem;
5459 }
5460 
5461 template <typename T>
5462 void
5463 UT_VoxelMipMap<T>::traverseTopDown(Callback function, void *data) const
5464 {
5465  doTraverse(0, 0, 0, 0, function, data);
5466 }
5467 
5468 template <typename T>
5469 void
5470 UT_VoxelMipMap<T>::doTraverse(int x, int y, int z, int level,
5471  Callback function, void *data) const
5472 {
5473  UT_StackBuffer<T> tval(myLevels.entries());
5474  bool isfinal;
5475  UT_VoxelArray<T> *vox;
5476  int shift;
5477 
5478  if (level == myNumLevels)
5479  {
5480  isfinal = true;
5481  vox = myBaseLevel;
5482  tval[0] = (*vox)(x, y, z);
5483  for (int i = 1; i < myLevels.entries(); i++)
5484  tval[i] = tval[0];
5485  }
5486  else
5487  {
5488  isfinal = false;
5489  for (int i = 0; i < myLevels.entries(); i++)
5490  {
5491  vox = myLevels(i)[level];
5492  tval[i] = (*vox)(x, y, z);
5493  }
5494  }
5495 
5496  shift = myNumLevels - level;
5497 
5498  UT_BoundingBox box;
5499 
5500  box.initBounds(SYSmin(x << shift, myBaseLevel->getXRes()),
5501  SYSmin(y << shift, myBaseLevel->getYRes()),
5502  SYSmin(z << shift, myBaseLevel->getZRes()));
5503  box.enlargeBounds(SYSmin((x+1) << shift, myBaseLevel->getXRes()),
5504  SYSmin((y+1) << shift, myBaseLevel->getYRes()),
5505  SYSmin((z+1) << shift, myBaseLevel->getZRes()));
5506 
5507  if (!function(tval, box, isfinal, data))
5508  {
5509  // Function asked for an early exit.
5510  // Give it.
5511  return;
5512  }
5513 
5514  // We now want to proceed to the next level.
5515  if (isfinal)
5516  return;
5517 
5518  level++;
5519  shift--;
5520  x <<= 1;
5521  y <<= 1;
5522  z <<= 1;
5523 
5524  bool xinc, yinc, zinc;
5525 
5526  // Determine which of our children are valid.
5527  if ( ((x+1) << shift) < myBaseLevel->getXRes() )
5528  xinc = true;
5529  else
5530  xinc = false;
5531  if ( ((y+1) << shift) < myBaseLevel->getYRes() )
5532  yinc = true;
5533  else
5534  yinc = false;
5535  if ( ((z+1) << shift) < myBaseLevel->getZRes() )
5536  zinc = true;
5537  else
5538  zinc = false;
5539 
5540  //
5541  // Traverse all valid children. Note that this is deliberately a gray
5542  // order traversal for full nodes.
5543  //
5544 
5545  doTraverse(x, y, z, level, function, data);
5546 
5547  if (yinc)
5548  {
5549  doTraverse(x, y+1, z, level, function, data);
5550  if (zinc)
5551  doTraverse(x, y+1, z+1, level, function, data);
5552  }
5553  if (zinc)
5554  doTraverse(x, y, z+1, level, function, data);
5555 
5556  if (xinc)
5557  {
5558  if (zinc)
5559  doTraverse(x+1, y, z+1, level, function, data);
5560  doTraverse(x+1, y, z, level, function, data);
5561  if (yinc)
5562  {
5563  doTraverse(x+1, y+1, z, level, function, data);
5564  if (zinc)
5565  doTraverse(x+1, y+1, z+1, level, function, data);
5566  }
5567  }
5568 }
5569 
5570 template <typename T>
5571 template <typename OP>
5572 void
5574 {
5575  doTraverse(0, 0, 0, numLevels()-1, op);
5576 }
5577 
5578 template <typename T>
5579 template <typename OP>
5580 void
5581 UT_VoxelMipMap<T>::doTraverse(int x, int y, int z, int level,
5582  OP &op) const
5583 {
5584  int shift;
5585 
5586  shift = level;
5587 
5588  UT_BoundingBoxI box;
5589 
5590  box.initBounds(SYSmin(x << shift, myBaseLevel->getXRes()),
5591  SYSmin(y << shift, myBaseLevel->getYRes()),
5592  SYSmin(z << shift, myBaseLevel->getZRes()));
5593  box.enlargeBounds(SYSmin((x+1) << shift, myBaseLevel->getXRes()),
5594  SYSmin((y+1) << shift, myBaseLevel->getYRes()),
5595  SYSmin((z+1) << shift, myBaseLevel->getZRes()));
5596 
5597  if (!op(box, level))
5598  {
5599  // Function asked for an early exit.
5600  // Give it.
5601  return;
5602  }
5603 
5604  level--;
5605  // We now want to proceed to the next level.
5606  if (level < 0)
5607  return;
5608 
5609  x <<= 1;
5610  y <<= 1;
5611  z <<= 1;
5612 
5613  bool xinc, yinc, zinc;
5614 
5615  // Determine which of our children are valid.
5616  if ( ((x+1) << level) < myBaseLevel->getXRes() )
5617  xinc = true;
5618  else
5619  xinc = false;
5620  if ( ((y+1) << level) < myBaseLevel->getYRes() )
5621  yinc = true;
5622  else
5623  yinc = false;
5624  if ( ((z+1) << level) < myBaseLevel->getZRes() )
5625  zinc = true;
5626  else
5627  zinc = false;
5628 
5629  //
5630  // Traverse all valid children. Note that this is deliberately a gray
5631  // order traversal for full nodes.
5632  //
5633 
5634  doTraverse(x, y, z, level, op);
5635 
5636  if (yinc)
5637  {
5638  doTraverse(x, y+1, z, level, op);
5639  if (zinc)
5640  doTraverse(x, y+1, z+1, level, op);
5641  }
5642  if (zinc)
5643  doTraverse(x, y, z+1, level, op);
5644 
5645  if (xinc)
5646  {
5647  if (zinc)
5648  doTraverse(x+1, y, z+1, level, op);
5649  doTraverse(x+1, y, z, level, op);
5650  if (yinc)
5651  {
5652  doTraverse(x+1, y+1, z, level, op);
5653  if (zinc)
5654  doTraverse(x+1, y+1, z+1, level, op);
5655  }
5656  }
5657 }
5658 
5659 template <typename T>
5660 template <typename OP>
5661 void
5663 {
5664  doTraverseSorted(0, 0, 0, numLevels()-1, op);
5665 }
5666 
5668 {
5669 public:
5670  float value;
5671  int key;
5672 };
5673 
5675 {
5676  return lhs.value < rhs.value;
5677 }
5678 
5679 template <typename T>
5680 template <typename OP>
5681 void
5682 UT_VoxelMipMap<T>::doTraverseSorted(int x, int y, int z, int level,
5683  OP &op) const
5684 {
5685  UT_BoundingBoxI box;
5686 
5687  box.initBounds(SYSmin(x << level, myBaseLevel->getXRes()),
5688  SYSmin(y << level, myBaseLevel->getYRes()),
5689  SYSmin(z << level, myBaseLevel->getZRes()));
5690  box.enlargeBounds(SYSmin((x+1) << level, myBaseLevel->getXRes()),
5691  SYSmin((y+1) << level, myBaseLevel->getYRes()),
5692  SYSmin((z+1) << level, myBaseLevel->getZRes()));
5693 
5694  if (!op(box, level))
5695  {
5696  // Function asked for an early exit.
5697  // Give it.
5698  return;
5699  }
5700 
5701  level--;
5702  // We now want to proceed to the next level.
5703  if (level < 0)
5704  return;
5705 
5706  x <<= 1;
5707  y <<= 1;
5708  z <<= 1;
5709 
5710  bool xinc, yinc, zinc;
5711 
5712  // Determine which of our children are valid.
5713  if ( ((x+1) << level) < myBaseLevel->getXRes() )
5714  xinc = true;
5715  else
5716  xinc = false;
5717  if ( ((y+1) << level) < myBaseLevel->getYRes() )
5718  yinc = true;
5719  else
5720  yinc = false;
5721  if ( ((z+1) << level) < myBaseLevel->getZRes() )
5722  zinc = true;
5723  else
5724  zinc = false;
5725 
5726  UT_BoundingBoxI boxes[8];
5727  int numboxes = 0;
5728 
5729  // Build all of our bounding boxes.
5730  for (int dz = 0; dz < 2; dz++)
5731  {
5732  if (dz && !zinc)
5733  continue;
5734  for (int dy = 0; dy < 2; dy++)
5735  {
5736  if (dy && !yinc)
5737  continue;
5738  for (int dx = 0; dx < 2; dx++)
5739  {
5740  if (dx && !xinc)
5741  continue;
5742  boxes[numboxes].initBounds(
5743  SYSmin((x+dx) << level, myBaseLevel->getXRes()),
5744  SYSmin((y+dy) << level, myBaseLevel->getYRes()),
5745  SYSmin((z+dz) << level, myBaseLevel->getZRes()));
5746  boxes[numboxes].enlargeBounds(
5747  SYSmin((x+dx+1) << level, myBaseLevel->getXRes()),
5748  SYSmin((y+dy+1) << level, myBaseLevel->getYRes()),
5749  SYSmin((z+dz+1) << level, myBaseLevel->getZRes()));
5750  numboxes++;
5751  }
5752  }
5753  }
5754 
5755  ut_VoxelMipMapSortCompare sortstats[8];
5756  for (int i = 0; i < numboxes; i++)
5757  {
5758  sortstats[i].value = op.sortValue(boxes[i], level);
5759  sortstats[i].key = i;
5760  }
5761  std::stable_sort(sortstats, &sortstats[numboxes]);
5762 
5763  for (int i = 0; i < numboxes; i++)
5764  {
5765  int whichbox = sortstats[i].key;
5766  doTraverseSorted(boxes[whichbox](0,0)>>level, boxes[whichbox](1,0)>>level, boxes[whichbox](2,0)>>level, level, op);
5767  }
5768 }
5769 
5770 template <typename T>
5771 void
5773 {
5774  myNumLevels = 0;
5775  myBaseLevel = 0;
5776  myOwnBase = false;
5777 }
5778 
5779 template <typename T>
5780 void
5782 {
5783  for (exint i = 0; i < myLevels.entries(); i++)
5784  {
5785  for (exint level = 0; level < myNumLevels; level++)
5786  delete myLevels(i)[level];
5787  delete [] myLevels(i);
5788  }
5789  myLevels.entries(0);
5790 
5791  if (myOwnBase)
5792  delete myBaseLevel;
5793 
5794  initializePrivate();
5795 }
5796 
5797 //
5798 // UT_VoxelArrayIterator implementation
5799 //
5800 template <typename T>
5802 {
5803  myArray = 0;
5804  myHandle.resetHandle();
5805  myCurTile = -1;
5806  myShouldCompressOnExit = false;
5807  myUseTileList = false;
5808  myJobInfo = 0;
5809  myInterrupt = 0;
5810 }
5811 
5812 template <typename T>
5814 {
5815  myShouldCompressOnExit = false;
5816  myUseTileList = false;
5817  myJobInfo = 0;
5818  setArray(vox);
5819  myInterrupt = 0;
5820 }
5821 
5822 template <typename T>
5824 {
5825  myShouldCompressOnExit = false;
5826  myUseTileList = false;
5827  myJobInfo = 0;
5828  setHandle(handle);
5829  myInterrupt = 0;
5830 }
5831 
5832 template <typename T>
5834 {
5835 }
5836 
5837 template <typename T>
5838 void
5840 {
5841  int numtiles;
5842  fpreal tileperrange;
5843 
5844  // Must already have been set.
5845  UT_ASSERT(myArray);
5846  if (!myArray)
5847  return;
5848 
5849  // Divide our tiles in to num ranges groups.
5850  numtiles = myArray->numTiles();
5851 
5852  // If we are using a tile list, use it instead
5853  if (myUseTileList)
5854  numtiles = myTileList.entries();
5855 
5856  // Trivial case.
5857  if (!numtiles)
5858  {
5859  myTileStart = 0;
5860  myTileEnd = 0;
5861  return;
5862  }
5863 
5864  // Sanity check idx.
5865  if (idx < 0 || idx >= numranges)
5866  {
5867  UT_ASSERT(!"Idx out of bounds");
5868  myTileStart = -1;
5869  myTileEnd = -1;
5870  return;
5871  }
5872 
5873  // Sanity test numranges
5874  if (numranges < 1)
5875  {
5876  UT_ASSERT(!"Invalid range count!");
5877  numranges = 1;
5878  }
5879 
5880  // Compute tile per range.
5881  // We use floating point so, if you have 15 tiles and
5882  // 16 processors, you don't decide to do all 15
5883  // tiles on one processor.
5884  tileperrange = (fpreal)numtiles / (fpreal)numranges;
5885 
5886  myTileStart = (int) SYSfloor(idx * tileperrange);
5887  myTileEnd = (int) SYSfloor((idx+1) * tileperrange);
5888 
5889  // Special case to ensure we always get the last tile,
5890  // despite the evils of floating point.
5891  if (idx == numranges-1)
5892  myTileEnd = numtiles;
5893 
5894  UT_ASSERT(myTileStart >= 0);
5895  UT_ASSERT(myTileEnd <= numtiles);
5896 }
5897 
5898 template <typename T>
5899 void
5901 {
5902  // Must already have been set.
5903  UT_ASSERT(myArray);
5904  if (!myArray)
5905  return;
5906 
5907  // If there is one thread, don't bother
5908  if (info.numJobs() == 1)
5909  myJobInfo = 0;
5910  else
5911  myJobInfo = &info;
5912 }
5913 
5914 template <typename T>
5915 void
5917 {
5918  UT_Vector3 pmin, pmax, vmin, vmax;
5919 
5920  pmin = bbox.minvec();
5921  pmax = bbox.maxvec();
5922 
5923  pmin.x() = SYSmax(pmin.x(), 0.0f);
5924  pmin.y() = SYSmax(pmin.y(), 0.0f);
5925  pmin.z() = SYSmax(pmin.z(), 0.0f);
5926  pmax.x() = SYSmin(pmax.x(), 1.0f);
5927  pmax.y() = SYSmin(pmax.y(), 1.0f);
5928  pmax.z() = SYSmin(pmax.z(), 1.0f);
5929 
5930  myArray->posToIndex(pmin, vmin);
5931  myArray->posToIndex(pmax, vmax);
5932 
5933  restrictToBBox(SYSfloor(vmin.x()), SYSceil(vmax.x()),
5934  SYSfloor(vmin.y()), SYSceil(vmax.y()),
5935  SYSfloor(vmin.z()), SYSceil(vmax.z()));
5936 }
5937 
5938 template <typename T>
5939 void
5941  int ymin, int ymax,
5942  int zmin, int zmax)
5943 {
5944  int xres, yres, zres, x, y, z;
5945 
5946  xres = myArray->getXRes();
5947  yres = myArray->getYRes();
5948  zres = myArray->getZRes();
5949 
5950  myTileList.entries(0);
5951  myUseTileList = true;
5952  // Make sure we have any tiles that are valid.
5953  if (xmin < xres && xmax >= 0 &&
5954  ymin < yres && ymax >= 0 &&
5955  zmin < zres && zmax >= 0)
5956  {
5957  // Clamp into valid range.
5958  myArray->clampIndex(xmin, ymin, zmin);
5959  myArray->clampIndex(xmax, ymax, zmax);
5960 
5961  // Convert to tiles...
5962  xmin >>= TILEBITS;
5963  ymin >>= TILEBITS;
5964  zmin >>= TILEBITS;
5965  xmax >>= TILEBITS;
5966  ymax >>= TILEBITS;
5967  zmax >>= TILEBITS;
5968 
5969  // Check if we are all accounted for.
5970  if (myArray->numTiles() == (xmax-xmin+1)*(ymax-ymin+1)*(zmax-zmin+1))
5971  {
5972  UT_ASSERT(xmin == 0 && ymin == 0 && zmin == 0);
5973  // No need for a tile list, just run everything!
5974  myUseTileList = false;
5975  }
5976  else
5977  {
5978  // Iterate over all active tiles, adding to our list.
5979  for (z = zmin; z <= zmax; z++)
5980  {
5981  for (y = ymin; y <= ymax; y++)
5982  {
5983  for (x = xmin; x <= xmax; x++)
5984  {
5985  myTileList.append(myArray->xyzTileToLinear(x, y, z));
5986  }
5987  }
5988  }
5989  }
5990  }
5991 
5992  // Recompute our ranges.
5993  myCurTile = -1;
5994  setPartialRange(0, 1);
5995 }
5996 
5997 
5998 template <typename T>
5999 void
6001 {
6002  // Ensure we have at least one tile in each direction
6003  if (!myArray ||
6004  !myArray->getRes(0) || !myArray->getRes(1) || !myArray->getRes(2))
6005  {
6006  myCurTile = -1;
6007  return;
6008  }
6009 
6010  if (myUseTileList)
6011  {
6012  if (myJobInfo)
6013  myCurTileListIdx = myJobInfo->nextTask();
6014  else
6015  myCurTileListIdx = myTileStart;
6016  if (myCurTileListIdx < 0 || myCurTileListIdx >= myTileEnd)
6017  {
6018  myCurTile = -1;
6019  return;
6020  }
6021  myCurTile = myTileList(myCurTileListIdx);
6022  }
6023  else
6024  {
6025  if (myJobInfo)
6026  myCurTile = myJobInfo->nextTask();
6027  else
6028  myCurTile = myTileStart;
6029  // If this is an empty range, we quit right away.
6030  if (myCurTile < 0 || myCurTile >= myTileEnd)
6031  {
6032  myCurTile = -1;
6033  return;
6034  }
6035  }
6036 
6037  UT_VoxelTile<T> *tile;
6038 
6039  tile = myArray->getLinearTile(myCurTile);
6040 
6041  // Get the tile position
6042  if (myCurTile)
6043  {
6044  myArray->linearTileToXYZ(myCurTile,
6045  myTilePos[0], myTilePos[1], myTilePos[2]);
6046  myPos[0] = TILESIZE * myTilePos[0];
6047  myPos[1] = TILESIZE * myTilePos[1];
6048  myPos[2] = TILESIZE * myTilePos[2];
6049  }
6050  else
6051  {
6052  myTilePos[0] = 0;
6053  myTilePos[1] = 0;
6054  myTilePos[2] = 0;
6055  myPos[0] = 0;
6056  myPos[1] = 0;
6057  myPos[2] = 0;
6058  }
6059 
6060  myTileLocalPos[0] = 0;
6061  myTileLocalPos[1] = 0;
6062  myTileLocalPos[2] = 0;
6063 
6064  myTileSize[0] = tile->xres();
6065  myTileSize[1] = tile->yres();
6066  myTileSize[2] = tile->zres();
6067 }
6068 
6069 template <typename T>
6070 void
6072 {
6073  if (myInterrupt && myInterrupt->opInterrupt())
6074  {
6075  myCurTile = -1;
6076  return;
6077  }
6078  if (myUseTileList)
6079  {
6080  if (getCompressOnExit())
6081  {
6082  // Verify our last tile was a legitimate one.
6083  if (myCurTile >= 0 && myCurTileListIdx < myTileEnd)
6084  {
6085  myArray->getLinearTile(myCurTile)->tryCompress(myArray->getCompressionOptions());
6086  }
6087  }
6088 
6089  // Advance our myCurTileListIdx and rebuild from
6090  if (myJobInfo)
6091  myCurTileListIdx = myJobInfo->nextTask();
6092  else
6093  myCurTileListIdx++;
6094  if (myCurTileListIdx >= myTileEnd)
6095  {
6096  myCurTile = -1;
6097  return;
6098  }
6099 
6100  myCurTile = myTileList(myCurTileListIdx);
6101 
6102  myArray->linearTileToXYZ(myCurTile,
6103  myTilePos[0], myTilePos[1], myTilePos[2]);
6104 
6105  UT_VoxelTile<T> *tile;
6106 
6107  tile = myArray->getLinearTile(myCurTile);
6108  myTileLocalPos[0] = 0;
6109  myTileLocalPos[1] = 0;
6110  myTileLocalPos[2] = 0;
6111  myTileSize[0] = tile->xres();
6112  myTileSize[1] = tile->yres();
6113  myTileSize[2] = tile->zres();
6114 
6115  myPos[0] = TILESIZE * myTilePos[0];
6116  myPos[1] = TILESIZE * myTilePos[1];
6117  myPos[2] = TILESIZE * myTilePos[2];
6118  return;
6119  }
6120 
6121  // If requested, compress our last tile.
6122  if (getCompressOnExit())
6123  {
6124  // Verify our last tile was a legitimate one.
6125  if (myCurTile >= 0 && myCurTile < myTileEnd)
6126  {
6127  myArray->getLinearTile(myCurTile)->tryCompress(myArray->getCompressionOptions());
6128  }
6129  }
6130 
6131  if (myJobInfo)
6132  {
6133  myCurTile = myJobInfo->nextTask();
6134  if (myCurTile >= myTileEnd)
6135  {
6136  myCurTile = -1;
6137  return;
6138  }
6139  myArray->linearTileToXYZ(myCurTile,
6140  myTilePos[0], myTilePos[1], myTilePos[2]);
6141  }
6142  else
6143  {
6144  // Advance our myTilePos, then rebuild everything from there.
6145  myTilePos[0]++;
6146  if (myTilePos[0] >= myArray->getTileRes(0))
6147  {
6148  myTilePos[0] = 0;
6149  myTilePos[1]++;
6150  if (myTilePos[1] >= myArray->getTileRes(1))
6151  {
6152  myTilePos[1] = 0;
6153  myTilePos[2]++;
6154  if (myTilePos[2] >= myArray->getTileRes(2))
6155  {
6156  // All done!
6157  myCurTile = -1;
6158  return;
6159  }
6160  }
6161  }
6162  myCurTile = myArray->xyzTileToLinear(myTilePos[0], myTilePos[1], myTilePos[2]);
6163  }
6164 
6165  UT_VoxelTile<T> *tile;
6166 
6167  // Rebuild our settings from this tile.
6168 
6169  // See if we hit our end.
6170  if (myCurTile >= myTileEnd)
6171  {
6172  myCurTile = -1;
6173  return;
6174  }
6175 
6176  tile = myArray->getLinearTile(myCurTile);
6177  myTileLocalPos[0] = 0;
6178  myTileLocalPos[1] = 0;
6179  myTileLocalPos[2] = 0;
6180  myTileSize[0] = tile->xres();
6181  myTileSize[1] = tile->yres();
6182  myTileSize[2] = tile->zres();
6183 
6184  myPos[0] = TILESIZE * myTilePos[0];
6185  myPos[1] = TILESIZE * myTilePos[1];
6186  myPos[2] = TILESIZE * myTilePos[2];
6187 }
6188 
6189 
6190 template <typename T>
6191 void
6193 {
6194  myTileLocalPos[0] = myTileSize[0];
6195  myTileLocalPos[1] = myTileSize[1];
6196  myTileLocalPos[2] = myTileSize[2];
6197 }
6198 
6199 template <typename T>
6200 template <typename OP>
6201 void
6203 {
6204  rewind();
6205 
6206  while (!atEnd())
6207  {
6208  UT_VoxelTile<T> *tile;
6209 
6210  tile = myArray->getLinearTile(myCurTile);
6211 
6212  if (!tile->isSimpleCompression())
6213  tile->uncompress();
6214 
6215  if (tile->isConstant())
6216  {
6217  // Woohoo! Too simple!
6218  T val;
6219 
6220  val = tile->rawData()[0];
6221 
6222  val = op(val);
6223 
6224  tile->rawData()[0] = val;
6225  }
6226  else
6227  {
6228  T *val;
6229 
6230  val = tile->rawData();
6231 
6232  int n = tile->numVoxels();
6233  for (int i = 0; i < n; i++)
6234  {
6235  val[i] = op(val[i]);
6236  }
6237  }
6238 
6239  advanceTile();
6240  }
6241 }
6242 
6243 template <typename T>
6244 template <typename OP, typename S>
6245 void
6247  const UT_VoxelArray<S> &a)
6248 {
6249  UT_ASSERT(myArray->isMatching(a));
6250 
6251  rewind();
6252 
6253  while (!atEnd())
6254  {
6255  UT_VoxelTile<S> *atile;
6256  UT_VoxelTile<T> *tile;
6257 
6258  tile = myArray->getLinearTile(myCurTile);
6259  atile = a.getLinearTile(myCurTile);
6260 
6261  if (!atile->isSimpleCompression())
6262  {
6263  // Have to use getValue...
6264  for (int z = 0; z < tile->zres(); z++)
6265  for (int y = 0; y < tile->yres(); y++)
6266  for (int x = 0; x < tile->xres(); x++)
6267  {
6268  S aval;
6269  T val;
6270 
6271  val = tile->operator()(x, y, z);
6272  aval = atile->operator()(x, y, z);
6273 
6274  val = op(val, aval);
6275 
6276  tile->setValue(x, y, z, val);
6277  }
6278  }
6279  else
6280  {
6281  if (!tile->isSimpleCompression())
6282  tile->uncompress();
6283 
6284  int ainc = atile->isConstant() ? 0 : 1;
6285 
6286  if (tile->isConstant())
6287  {
6288  if (!ainc)
6289  {
6290  // Woohoo! Too simple!
6291  S aval;
6292  T val;
6293 
6294  val = tile->rawData()[0];
6295  aval = atile->rawData()[0];
6296 
6297  val = op(val, aval);
6298 
6299  tile->rawData()[0] = val;
6300  }
6301  else
6302  {
6303  tile->uncompress();
6304  }
6305  }
6306 
6307  // In case we uncompressed ourself above...
6308  // we have to test again
6309  if (!tile->isConstant())
6310  {
6311  const S *aval;
6312  T *val;
6313 
6314  val = tile->rawData();
6315  aval = atile->rawData();
6316  ainc = atile->isConstant() ? 0 : 1;
6317 
6318  int n = tile->numVoxels();
6319  for (int i = 0; i < n; i++)
6320  {
6321  val[i] = op(val[i], *aval);
6322  aval += ainc;
6323  }
6324  }
6325  }
6326 
6327  advanceTile();
6328  }
6329 }
6330 
6331 template <typename T>
6332 template <typename OP>
6333 void
6335 {
6336  rewind();
6337 
6338  while (!atEnd())
6339  {
6340  UT_VoxelTile<T> *tile;
6341 
6342  tile = myArray->getLinearTile(myCurTile);
6343 
6344  if (!tile->isSimpleCompression())
6345  tile->uncompress();
6346 
6347  if (tile->isConstant())
6348  {
6349  T val;
6350 
6351  val = tile->rawData()[0];
6352 
6353  val = op(val, a);
6354 
6355  tile->rawData()[0] = val;
6356  }
6357  else
6358  {
6359  T *val;
6360 
6361  val = tile->rawData();
6362 
6363  int n = tile->numVoxels();
6364  for (int i = 0; i < n; i++)
6365  {
6366  val[i] = op(val[i], a);
6367  }
6368  }
6369 
6370  advanceTile();
6371  }
6372 }
6373 
6374 template <typename T>
6375 template <typename OP, typename S, typename R>
6376 void
6378  const UT_VoxelArray<S> &a,
6379  const UT_VoxelArray<R> &b)
6380 {
6381  UT_ASSERT(myArray->isMatching(a));
6382  UT_ASSERT(myArray->isMatching(b));
6383 
6384  rewind();
6385 
6386  while (!atEnd())
6387  {
6388  UT_VoxelTile<S> *atile;
6389  UT_VoxelTile<R> *btile;
6390  UT_VoxelTile<T> *tile;
6391 
6392  tile = myArray->getLinearTile(myCurTile);
6393  atile = a.getLinearTile(myCurTile);
6394  btile = b.getLinearTile(myCurTile);
6395 
6396  if (!atile->isSimpleCompression() || !btile->isSimpleCompression())
6397  {
6398  // Have to use getValue...
6399  for (int z = 0; z < tile->zres(); z++)
6400  for (int y = 0; y < tile->yres(); y++)
6401  for (int x = 0; x < tile->xres(); x++)
6402  {
6403  S aval;
6404  R bval;
6405  T val;
6406 
6407  val = tile->operator()(x, y, z);
6408  aval = atile->operator()(x, y, z);
6409  bval = btile->operator()(x, y, z);
6410 
6411  val = op(val, aval, bval);
6412 
6413  tile->setValue(x, y, z, val);
6414  }
6415  }
6416  else
6417  {
6418  if (!tile->isSimpleCompression())
6419  tile->uncompress();
6420 
6421  int ainc = atile->isConstant() ? 0 : 1;
6422  int binc = btile->isConstant() ? 0 : 1;
6423 
6424  if (tile->isConstant())
6425  {
6426  if (!ainc && !binc)
6427  {
6428  // Woohoo! Too simple!
6429  S aval;
6430