volumerender.h
説明を見る。
1 //
2 // Copyright (c) 2003-2011, MIST Project, Nagoya University
3 // All rights reserved.
4 //
5 // Redistribution and use in source and binary forms, with or without modification,
6 // are permitted provided that the following conditions are met:
7 //
8 // 1. Redistributions of source code must retain the above copyright notice,
9 // this list of conditions and the following disclaimer.
10 //
11 // 2. Redistributions in binary form must reproduce the above copyright notice,
12 // this list of conditions and the following disclaimer in the documentation
13 // and/or other materials provided with the distribution.
14 //
15 // 3. Neither the name of the Nagoya University nor the names of its contributors
16 // may be used to endorse or promote products derived from this software
17 // without specific prior written permission.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
20 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
21 // FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
22 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
24 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
25 // IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
26 // THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 //
28 
38 #ifndef __INCLUDE_MIST_VOLUMERENDER__
39 #define __INCLUDE_MIST_VOLUMERENDER__
40 
41 
42 #ifndef __INCLUDE_MIST_CONF_H__
43 #include "config/mist_conf.h"
44 #endif
45 
46 #ifndef __INCLUDE_MIST_H__
47 #include "mist.h"
48 #endif
49 
50 // カラー画像の設定を読み込む
51 #ifndef __INCLUDE_MIST_COLOR_H__
52 #include "config/color.h"
53 #endif
54 
55 #ifndef __INCLUDE_MIST_VECTOR__
56 #include "vector.h"
57 #endif
58 
59 #ifndef __INCLUDE_MIST_LIMITS__
60 #include "limits.h"
61 #endif
62 
63 
64 #ifndef __INCLUDE_MIST_THREAD__
65 #include "thread.h"
66 #endif
67 
68 #ifndef __INCLUDE_MIST_DISTANCE_TRANSFORM__
69 #include "filter/distance.h"
70 #endif
71 
72 
73 #include <vector>
74 
75 
76 // mist名前空間の始まり
77 _MIST_BEGIN
78 
79 
80 
81 
85 
86 
94 
95 
96 namespace volumerender
97 {
98  template < class T >
99  struct attribute
100  {
101  typedef T value_type;
102 
103  value_type pixel;
104  double alpha;
105  bool has_alpha;
106 
107  attribute( ) : pixel( ), alpha( 0 ), has_alpha( false ){ }
108 
109  //template < class TT >
110  //attribute( const attribute &a ) : pixel( a.pixel ), alpha( a.alpha ), has_alpha( a.has_alpha ){ }
111 
112  attribute( const attribute &a ) : pixel( a.pixel ), alpha( a.alpha ), has_alpha( a.has_alpha ){ }
113 
114  attribute( const value_type &pix, double a ) : pixel( pix ), alpha( a ), has_alpha( a > 0.0 ){ }
115 
116  };
117 
118  template < class T >
119  struct attribute_element
120  {
121  typedef T value_type;
122  typedef ptrdiff_t difference_type;
123 
124  value_type pixel;
125  difference_type value_lower;
126  difference_type value_upper;
127  double alpha_lower;
128  double alpha_upper;
129 
130  attribute_element( ) : pixel( ), value_lower( 0 ), value_upper( 0 ), alpha_lower( 0 ), alpha_upper( 0 ){ }
131  attribute_element( const value_type &pixel, difference_type vs, difference_type ve, double as, double ae )
132  : pixel( pixel ), value_lower( vs ), value_upper( ve ), alpha_lower( as ), alpha_upper( ae ){ }
133 
134  attribute_element( const attribute_element &a ) : pixel( a.pixel ), value_lower( a.value_lower ), value_upper( a.value_upper ), alpha_lower( a.alpha_lower ), alpha_upper( a.alpha_upper ){ }
135 
136  };
137 
138  template < class T >
139  class attribute_table : protected std::vector< attribute< T > >
140  {
141  public:
142  typedef std::vector< attribute< T > > base;
143  typedef attribute< T > attribute_type;
144  typedef attribute_element< T > attribute_element_type;
145  typedef typename base::allocator_type allocator_type;
146  typedef typename base::reference reference;
147  typedef typename base::const_reference const_reference;
148  typedef typename base::value_type value_type;
149  typedef typename base::size_type size_type;
150  typedef typename base::difference_type difference_type;
151  typedef typename base::pointer pointer;
152  typedef typename base::const_pointer const_pointer;
153 
154  typedef T pixel_type;
155  typedef std::vector< attribute_element_type > attribute_element_list;
156 
157  protected:
158  difference_type sindex_;
159  difference_type eindex_;
160  difference_type zero_index_;
161 
162  attribute_element_list attribute_element_table;
163 
164  public:
165  void append( const pixel_type &pixel, difference_type vs, difference_type ve, double as, double ae )
166  {
167  if( base::empty( ) )
168  {
169  return;
170  }
171  else if( vs > ve )
172  {
173  difference_type tmp = vs;
174  vs = ve;
175  ve = tmp;
176  }
177  vs = vs < sindex_ ? sindex_: vs;
178  vs = vs > eindex_ ? eindex_: vs;
179  ve = ve < sindex_ ? sindex_: ve;
180  ve = ve > eindex_ ? eindex_: ve;
181 
182  if( as > ae )
183  {
184  double tmp = as;
185  as = ae;
186  ae = tmp;
187  }
188 
189  attribute_element_table.push_back( attribute_element_type( pixel, vs, ve, as, ae ) );
190 
191  double step = ( ae - as ) / static_cast< double >( ve - vs + 1 );
192  for( difference_type i = vs ; i <= ve ; i++ )
193  {
194  value_type &p = at( i );
195  p.pixel = pixel;
196  p.alpha = as + static_cast< double >( i - vs ) * step;
197  p.has_alpha = p.alpha > 0.0;
198  }
199  }
200 
201  public:
202  const difference_type minimum( ) const
203  {
204  return( sindex_ );
205  }
206 
207  const difference_type maximum( ) const
208  {
209  return( eindex_ );
210  }
211 
212  const bool empty( ) const
213  {
214  return( eindex_ < sindex_ );
215  }
216 
217  value_type &operator []( difference_type index )
218  {
219  return( base::operator []( zero_index_ + index ) );
220  }
221 
222  value_type &at( difference_type index )
223  {
224  return( base::at( zero_index_ + index ) );
225  }
226 
227  const value_type &operator []( difference_type index ) const
228  {
229  return( base::operator []( zero_index_ + index ) );
230  }
231 
232  const value_type &at( difference_type index ) const
233  {
234  return( base::at( zero_index_ + index ) );
235  }
236 
237  bool has_alpha( difference_type index ) const
238  {
239  return( operator []( index ).has_alpha );
240  }
241 
242  const attribute_table &operator =( const attribute_table &a )
243  {
244  if( &a != this )
245  {
246  base::operator =( a );
247  sindex_ = a.sindex_;
248  eindex_ = a.eindex_;
249  zero_index_ = a.zero_index_;
250  attribute_element_table = a.attribute_element_table;
251  }
252  return( *this );
253  }
254 
255  const attribute_element_list &attribute_elements( ) const { return( attribute_element_table ); }
256 
257 
258  public:
259  void create( difference_type si, difference_type ei )
260  {
261  sindex_ = si;
262  eindex_ = ei;
263 
264  base::resize( ei - si + 1 );
265  zero_index_ = - si;
266 
267  fill( );
268  attribute_element_table.clear( );
269  }
270 
271  void clear( )
272  {
273  base::clear( );
274  sindex_ = 0;
275  eindex_ = -1;
276  zero_index_ = 0;
277  attribute_element_table.clear( );
278  }
279 
280  void fill( )
281  {
282  for( size_type i = 0 ; i < base::size( ) ; i++ )
283  {
284  base::operator []( i ) = attribute_type( );
285  }
286  attribute_element_table.clear( );
287  }
288 
289  attribute_table( ) : sindex_( 0 ), eindex_( -1 ), zero_index_( 0 ){ }
290 
291  attribute_table( difference_type si, difference_type ei ) : base( ei - si + 1 ), sindex_( si ), eindex_( ei ), zero_index_( 0 )
292  {
293  zero_index_ = - si;
294  }
295 
296  attribute_table( const attribute_table &a ) : base( a ), sindex_( a.sindex_ ), eindex_( a.eindex_ ), zero_index_( a.zero_index_ ), attribute_element_table( a.attribute_element_table )
297  {
298  }
299  };
300 
301  struct boundingbox
302  {
303  public:
304  double a;
305  double b;
306  double c;
307  double d;
308 
309  void clear( )
310  {
311  a = 0;
312  b = 0;
313  c = 0;
314  d = 0;
315  }
316 
317  boundingbox( ) : a( 0 ), b( 0 ), c( 0 ), d( 0 ){ }
318  boundingbox( double aa, double bb, double cc, double dd ) : a( aa ), b( bb ), c( cc ), d( dd ){ }
319  };
320 
321 
322  template < class T >
323  inline bool check_intersection( vector3< T > &p1, vector3< T > &p2, const boundingbox &box )
324  {
325  typedef vector3< T > vector_type;
326  typedef typename vector3< T >::value_type value_type;
327 
328  vector_type n( box.a, box.b, box.c );
329 
330  value_type d1 = p1.inner( n );
331  value_type d2 = p2.inner( n );
332  bool c1 = box.d + d1 >= 0;
333  bool c2 = box.d + d2 >= 0;
334 
335  if( c1 && c2 )
336  {
337  return( true );
338  }
339  else if( !c1 && !c2 )
340  {
341  return( false );
342  }
343  vector_type v1 = p1;
344  vector_type v2 = p2;
345  if( !c1 )
346  {
347  p1 = v2 + ( v1 - v2 ) * ( ( box.d + d2 ) / ( d2 - d1 ) );
348  }
349  if( !c2 )
350  {
351  p2 = v1 + ( v2 - v1 ) * ( ( box.d + d1 ) / ( d1 - d2 ) );
352  }
353  return( true );
354  }
355 
356  template < class T >
357  inline bool check_intersection( vector3< T > &p1, vector3< T > &p2, const boundingbox &box, vector3< T > &normal )
358  {
359  typedef vector3< T > vector_type;
360  typedef typename vector3< T >::value_type value_type;
361 
362  vector_type n( box.a, box.b, box.c );
363 
364  value_type d1 = p1.inner( n );
365  value_type d2 = p2.inner( n );
366  bool c1 = box.d + d1 >= 0;
367  bool c2 = box.d + d2 >= 0;
368 
369  if( c1 && c2 )
370  {
371  return( true );
372  }
373  else if( !c1 && !c2 )
374  {
375  return( false );
376  }
377  vector_type v1 = p1;
378  vector_type v2 = p2;
379  if( !c1 )
380  {
381  p1 = v2 + ( v1 - v2 ) * ( ( box.d + d2 ) / ( d2 - d1 ) );
382  normal = n;
383  }
384  if( !c2 )
385  {
386  p2 = v1 + ( v2 - v1 ) * ( ( box.d + d1 ) / ( d1 - d2 ) );
387  }
388  return( true );
389  }
390 
391  struct parameter
392  {
393  typedef vector3< double > vector_type;
394 
395  vector_type pos; // 画像座標系でのカメラ位置
396  vector_type dir; // 画像座標系でのカメラ視線方向
397  vector_type up; // 画像座標系でのカメラ上向き方向
398 
399  bool perspective_view; // 透視投影で描画するかどうか
400  bool mirror_view; // 鏡に映った画像として描画するかどうか
401  bool value_interpolation; // 値補間を利用して描画するかどうか（false の場合は色補間を利用）
402  double fovy; // 投影面のY軸方向の視野角
403  double ambient_ratio; // 環境光の強さ（最大1）
404  double diffuse_ratio; // 拡散反射光の強さ（最大1）
405  double light_attenuation; // 光の減衰係数
406  double sampling_step; // レイキャスティング時のサンプリング間隔
407  double termination; // レイキャスティングの終了条件（積算不透明度が個の値未満になるとキャスティングを打ち切り）
408  double specular; // 鏡面反射光の強さ（最大1）
409  double background_R; // 背景色のR成分
410  double background_G; // 背景色のG成分
411  double background_B; // 背景色のB成分
412 
413  vector_type offset; // バウンディングボックスの中心座標
414  boundingbox box[ 6 ]; // バウンディングボックス
415 
416  parameter( ) : perspective_view( true ), mirror_view( false ), value_interpolation( true ), fovy( 80.0 ), ambient_ratio( 0.4 ), diffuse_ratio( 0.6 ),
417  light_attenuation( 0.0 ), sampling_step( 1.0 ), termination( 0.01 ), specular( 1.0 ), background_R( 0.0 ), background_G( 0.0 ), background_B( 0.0 )
418  {
419  }
420  };
421 
422  inline std::ostream &operator <<( std::ostream &out, const parameter &p )
423  {
424  out << "Pos = ( " << p.pos << " )" << std::endl;
425  out << "Dir = ( " << p.dir << " )" << std::endl;
426  out << "Up = ( " << p.up << " )" << std::endl;
427 
428  out << "BOX Center = ( " << p.offset << " )" << std::endl;
429 
430  out << "Perspective = " << p.perspective_view << std::endl;
431  out << "Mirror View = " << p.mirror_view << std::endl;
432  out << "ValueInterpolation = " << p.value_interpolation << std::endl;
433  out << "Fovy = " << p.fovy << std::endl;
434  out << "Ambient = " << p.ambient_ratio << std::endl;
435  out << "Diffuse = " << p.diffuse_ratio << std::endl;
436  out << "Specular = " << p.specular << std::endl;
437  out << "LightAtten = " << p.light_attenuation << std::endl;
438  out << "SamplingStep = " << p.sampling_step << std::endl;
439  out << "Termination = " << p.termination << std::endl;
440 
441  return( out );
442  }
443 
444  struct no_depth_map
445  {
446  typedef ptrdiff_t difference_type;
447  double operator()( difference_type i, difference_type j, difference_type k ) const
448  {
449  return( 2.0 );
450  }
451  };
452 
453  template < class DepthMap >
454  struct depth_map
455  {
456  typedef DepthMap depth_map_type;
457  typedef typename depth_map_type::difference_type difference_type;
458 
459  const depth_map_type &depth_map_;
460 
461  depth_map( const depth_map_type &dmap ) : depth_map_( dmap )
462  {
463  }
464 
465  double operator()( difference_type i, difference_type j, difference_type k ) const
466  {
467  double l = depth_map_( i >> 2, j >> 2, k >> 2 );
468  return( l < 1.0 ? 2.0 : l * 4.0 - 2.0 );
469  }
470  };
471 
472 
473  // 倍精度浮動小数を整数に丸める
474  // 正の数の場合は，通常の int キャストと同じ動作をする
475  // 負の数のときは，値の小さいほうに丸められる
476  inline long to_integer( double val )
477  {
478  return( static_cast< long >( val ) );
479  //val += 68719476736.0 * 1.5;
480  //return( ( ( long * )&val )[ 0 ] >> 16 );
481  }
482 
483  template < class T >
484  inline void normalize_vector( vector3< T > &v )
485  {
486  double len = std::sqrt( v.x * v.x + v.y * v.y + v.z * v.z ) + 1.0e-10;
487  v.x /= len;
488  v.y /= len;
489  v.z /= len;
490  }
491 
492  inline double _3( double v )
493  {
494  if( v < 0.0 )
495  {
496  return( - std::pow( -v, 1.0 / 3.0 ) );
497  }
498  else if( v > 0.0 )
499  {
500  return( std::pow( v, 1.0 / 3.0 ) );
501  }
502  else
503  {
504  return( 0.0 );
505  }
506  }
507 
508  inline double __maximum__( double v1, double v2 )
509  {
510  return( v1 > v2 ? v1 : v2 );
511  }
512 
513  inline double __maximum__( double v1, double v2, double v3 )
514  {
515  if( v1 > v2 )
516  {
517  if( v1 > v3 )
518  {
519  return( v1 );
520  }
521  else
522  {
523  return( v3 );
524  }
525  }
526  else if( v2 > v3 )
527  {
528  return( v2 );
529  }
530  else
531  {
532  return( v3 );
533  }
534  }
535 
536  inline double __minimum__( double v1, double v2 )
537  {
538  return( v1 < v2 ? v1 : v2 );
539  }
540 
541  inline double __minimum__( double v1, double v2, double v3 )
542  {
543  if( v1 < v2 )
544  {
545  if( v1 < v3 )
546  {
547  return( v1 );
548  }
549  else
550  {
551  return( v3 );
552  }
553  }
554  else if( v2 < v3 )
555  {
556  return( v2 );
557  }
558  else
559  {
560  return( v3 );
561  }
562  }
563 
564  inline size_t solve_equation( double p, double q, double &v1, double &v2, double &v3 )
565  {
566  const double pai = 3.1415926535897932384626433832795;
567  double p_3 = p / 3.0;
568  double q_2 = q / 2.0;
569  double D = - ( q_2 * q_2 + p_3 * p_3 * p_3 );
570 
571  if( D > 0.0 )
572  {
573  double r = std::sqrt( - 4.0 * p / 3.0 );
574  double t1 = 1.0 / ( 3.0 * std::cos( 3.0 * q / ( p * r ) ) );
575  double t2 = 2.0 * pai / 3.0 + t1;
576  double t3 = 2.0 * pai / 3.0 - t1;
577  v1 = r * std::cos( t1 );
578  v2 = r * std::cos( t2 );
579  v3 = r * std::cos( t3 );
580  return( 3 );
581  }
582  else if( D == 0.0 )
583  {
584  if( p == 0.0 )
585  {
586  v1 = v2 = v3 = 0.0;
587  return( 1 );
588  }
589  else
590  {
591  double r = std::sqrt( - 4.0 * p / 3.0 );
592  v1 = r;
593  v2 = r * std::cos( 2.0 * pai / 3.0 );
594  v3 = 0.0;
595  return( 2 );
596  }
597  }
598  //else if( D < 0.0 && p < 0.0 )
599  //{
600  // double _D = std::sqrt( -D );
601  // v1 = _3( - q / 2.0 + _D ) + _3( - q / 2.0 - _D );
602  // v2 = v3 = 0.0;
603  // return( 1 );
604  //}
605  else
606  {
607  double _D = std::sqrt( -D );
608  v1 = _3( - q / 2.0 + _D ) + _3( - q / 2.0 - _D );
609  v2 = v3 = 0.0;
610  return( 1 );
611  }
612  }
613 
614  inline size_t solve_equation( double a, double b, double c, double d, double &v1, double &v2, double &v3 )
615  {
616  if( a == 0.0 )
617  {
618  return( 0 );
619  }
620  else
621  {
622  return( solve_equation( ( - b * b / ( 3.0 * a ) + c ) / a, ( 2.0 * b * b * b / ( 27.0 * a * a ) - b * c / ( 3.0 * a ) + d ) / a, v1, v2, v3 ) );
623  }
624  }
625 }
626 
627 
629 // ボリュームレンダリンググループの終わり
630 
632 // 可視化グループの終わり
633 
634 namespace rendering_helper
635 {
636  template < class Array, class T >
637  class value_interpolation
638  {
639  protected:
640  typedef typename volumerender::parameter::vector_type vector_type;
641  typedef typename Array::size_type size_type;
642  typedef typename Array::const_pointer const_pointer;
643  typedef typename Array::difference_type difference_type;
644  typedef typename volumerender::attribute_table< T >::attribute_type attribute_type;
645 
646  public:
647  const Array &in;
648  const volumerender::parameter &p;
649  const volumerender::attribute_table< T > &table;
650  difference_type d0, d1, d2, d3, d4, d5, d6, d7, _1, _2, _3;
651 
652  value_interpolation( const Array &image, const volumerender::parameter &param, const volumerender::attribute_table< T > &tbl ) : in( image ), p( param ), table( tbl )
653  {
654  difference_type cx = in.width( ) / 2;
655  difference_type cy = in.height( ) / 2;
656  difference_type cz = in.depth( ) / 2;
657  const_pointer ppp = &in( cx, cy, cz );
658  d0 = 0;
659  d1 = &in( cx , cy + 1, cz ) - ppp;
660  d2 = &in( cx + 1, cy + 1, cz ) - ppp;
661  d3 = &in( cx + 1, cy , cz ) - ppp;
662  d4 = &in( cx , cy , cz + 1 ) - ppp;
663  d5 = &in( cx , cy + 1, cz + 1 ) - ppp;
664  d6 = &in( cx + 1, cy + 1, cz + 1 ) - ppp;
665  d7 = &in( cx + 1, cy , cz + 1 ) - ppp;
666  _1 = &in( cx + 1, cy , cz ) - ppp;
667  _2 = &in( cx , cy + 1, cz ) - ppp;
668  _3 = &in( cx , cy , cz + 1 ) - ppp;
669  }
670 
671  bool check( difference_type i, difference_type j, difference_type k ) const
672  {
673  const_pointer p = &in( i, j, k );
674 
675  // この位置における物体が不透明の場合は次のステップへ移行する
676  return( table.has_alpha( p[ d0 ] ) || table.has_alpha( p[ d1 ] ) ||
677  table.has_alpha( p[ d2 ] ) || table.has_alpha( p[ d3 ] ) ||
678  table.has_alpha( p[ d4 ] ) || table.has_alpha( p[ d5 ] ) ||
679  table.has_alpha( p[ d6 ] ) || table.has_alpha( p[ d7 ] ) );
680  }
681 
682  bool render( difference_type i, difference_type j, difference_type k, double xx, double yy, double zz, attribute_type &oc ) const
683  {
684  const_pointer p = &in( i, j, k );
685 
686  // CT値に対応する色と不透明度を取得
687  double ct = ( p[ d0 ] + ( p[ d3 ] - p[ d0 ] ) * xx ) + ( p[ d1 ] - p[ d0 ] + ( p[ d0 ] - p[ d1 ] + p[ d2 ] - p[ d3 ] ) * xx ) * yy;
688  ct += ( ( p[ d4 ] + ( p[ d7 ] - p[ d4 ] ) * xx ) + ( p[ d5 ] - p[ d4 ] + ( p[ d4 ] - p[ d5 ] + p[ d6 ] - p[ d7 ] ) * xx ) * yy - ct ) * zz;
689 
690  oc = table[ volumerender::to_integer( ct ) ];
691 
692  return( oc.has_alpha );
693  }
694 
695  vector_type normal( difference_type i, difference_type j, difference_type k, double xx, double yy, double zz ) const
696  {
697  const_pointer p0 = &in( i, j, k );
698  const_pointer p1 = p0 + d1;
699  const_pointer p2 = p0 + d2;
700  const_pointer p3 = p0 + d3;
701  const_pointer p4 = p0 + d4;
702  const_pointer p5 = p0 + d5;
703  const_pointer p6 = p0 + d6;
704  const_pointer p7 = p0 + d7;
705 
706  double n0x = p3[ 0 ] - p0[ -_1 ];
707  double n0y = p1[ 0 ] - p0[ -_2 ];
708  double n0z = p4[ 0 ] - p0[ -_3 ];
709  double n1x = p2[ 0 ] - p1[ -_1 ];
710  double n1y = p1[ _2 ] - p0[ 0 ];
711  double n1z = p5[ 0 ] - p1[ -_3 ];
712  double n2x = p2[ _1 ] - p1[ 0 ];
713  double n2y = p2[ _2 ] - p3[ 0 ];
714  double n2z = p6[ 0 ] - p2[ -_3 ];
715  double n3x = p3[ _1 ] - p0[ 0 ];
716  double n3y = p2[ 0 ] - p3[ -_2 ];
717  double n3z = p7[ 0 ] - p3[ -_3 ];
718  double n4x = p7[ 0 ] - p4[ -_1 ];
719  double n4y = p5[ 0 ] - p4[ -_2 ];
720  double n4z = p4[ _3 ] - p0[ 0 ];
721  double n5x = p6[ 0 ] - p5[ -_1 ];
722  double n5y = p5[ _2 ] - p4[ 0 ];
723  double n5z = p5[ _3 ] - p1[ 0 ];
724  double n6x = p6[ _1 ] - p5[ 0 ];
725  double n6y = p6[ _2 ] - p7[ 0 ];
726  double n6z = p6[ _3 ] - p2[ 0 ];
727  double n7x = p7[ _1 ] - p4[ 0 ];
728  double n7y = p6[ 0 ] - p7[ -_2 ];
729  double n7z = p7[ _3 ] - p3[ 0 ];
730 
731  double nx = ( n0x + ( n3x - n0x ) * xx ) + ( n1x - n0x + ( n0x - n1x + n2x - n3x ) * xx ) * yy;
732  nx += ( ( n4x + ( n7x - n4x ) * xx ) + ( n5x - n4x + ( n4x - n5x + n6x - n7x ) * xx ) * yy - nx ) * zz;
733  double ny = ( n0y + ( n3y - n0y ) * xx ) + ( n1y - n0y + ( n0y - n1y + n2y - n3y ) * xx ) * yy;
734  ny += ( ( n4y + ( n7y - n4y ) * xx ) + ( n5y - n4y + ( n4y - n5y + n6y - n7y ) * xx ) * yy - ny ) * zz;
735  double nz = ( n0z + ( n3z - n0z ) * xx ) + ( n1z - n0z + ( n0z - n1z + n2z - n3z ) * xx ) * yy;
736  nz += ( ( n4z + ( n7z - n4z ) * xx ) + ( n5z - n4z + ( n4z - n5z + n6z - n7z ) * xx ) * yy - nz ) * zz;
737 
738  nx /= in.reso1( );
739  ny /= in.reso2( );
740  nz /= in.reso3( );
741  double len = std::sqrt( nx * nx + ny * ny + nz * nz ) + type_limits< double >::tiny( );
742  nx /= len;
743  ny /= len;
744  nz /= len;
745 
746  return( vector_type( nx, ny, nz ) );
747  }
748  };
749 
750 
751  template < class Array, class T >
752  class color_interpolation
753  {
754  protected:
755  typedef typename volumerender::parameter::vector_type vector_type;
756  typedef typename Array::size_type size_type;
757  typedef typename Array::const_pointer const_pointer;
758  typedef typename Array::difference_type difference_type;
759  typedef typename volumerender::attribute_table< T >::attribute_type attribute_type;
760 
761  public:
762  const Array &in;
763  const volumerender::parameter &p;
764  const volumerender::attribute_table< T > &table;
765  difference_type d0, d1, d2, d3, d4, d5, d6, d7, _1, _2, _3;
766 
767  color_interpolation( const Array &image, const volumerender::parameter &param, const volumerender::attribute_table< T > &tbl ) : in( image ), p( param ), table( tbl )
768  {
769  difference_type cx = in.width( ) / 2;
770  difference_type cy = in.height( ) / 2;
771  difference_type cz = in.depth( ) / 2;
772  const_pointer ppp = &in( cx, cy, cz );
773  d0 = 0;
774  d1 = &in( cx , cy + 1, cz ) - ppp;
775  d2 = &in( cx + 1, cy + 1, cz ) - ppp;
776  d3 = &in( cx + 1, cy , cz ) - ppp;
777  d4 = &in( cx , cy , cz + 1 ) - ppp;
778  d5 = &in( cx , cy + 1, cz + 1 ) - ppp;
779  d6 = &in( cx + 1, cy + 1, cz + 1 ) - ppp;
780  d7 = &in( cx + 1, cy , cz + 1 ) - ppp;
781  _1 = &in( cx + 1, cy , cz ) - ppp;
782  _2 = &in( cx , cy + 1, cz ) - ppp;
783  _3 = &in( cx , cy , cz + 1 ) - ppp;
784  }
785 
786  bool check( difference_type i, difference_type j, difference_type k ) const
787  {
788  const_pointer p = &in( i, j, k );
789 
790  // この位置における物体が不透明の場合は次のステップへ移行する
791  return( table.has_alpha( p[ d0 ] ) || table.has_alpha( p[ d1 ] ) ||
792  table.has_alpha( p[ d2 ] ) || table.has_alpha( p[ d3 ] ) ||
793  table.has_alpha( p[ d4 ] ) || table.has_alpha( p[ d5 ] ) ||
794  table.has_alpha( p[ d6 ] ) || table.has_alpha( p[ d7 ] ) );
795  }
796 
797  bool render( difference_type i, difference_type j, difference_type k, double xx, double yy, double zz, attribute_type &oc ) const
798  {
799  const_pointer p = &in( i, j, k );
800 
801  // 近傍の8点に対応する色と不透明度を取得
802  attribute_type oc0 = table[ volumerender::to_integer( p[ d0 ] ) ];
803  attribute_type oc1 = table[ volumerender::to_integer( p[ d1 ] ) ];
804  attribute_type oc2 = table[ volumerender::to_integer( p[ d2 ] ) ];
805  attribute_type oc3 = table[ volumerender::to_integer( p[ d3 ] ) ];
806  attribute_type oc4 = table[ volumerender::to_integer( p[ d4 ] ) ];
807  attribute_type oc5 = table[ volumerender::to_integer( p[ d5 ] ) ];
808  attribute_type oc6 = table[ volumerender::to_integer( p[ d6 ] ) ];
809  attribute_type oc7 = table[ volumerender::to_integer( p[ d7 ] ) ];
810 
811  int number_of_visible_voxels = oc0.has_alpha;
812  number_of_visible_voxels += oc1.has_alpha;
813  number_of_visible_voxels += oc2.has_alpha;
814  number_of_visible_voxels += oc3.has_alpha;
815  number_of_visible_voxels += oc4.has_alpha;
816  number_of_visible_voxels += oc5.has_alpha;
817  number_of_visible_voxels += oc6.has_alpha;
818  number_of_visible_voxels += oc7.has_alpha;
819 
820  if( number_of_visible_voxels == 0 )
821  {
822  return( false );
823  }
824  else
825  {
826  // まず平均的な色を決定する
827  oc.pixel = ( oc0.pixel * oc0.has_alpha + oc1.pixel * oc1.has_alpha
828  + oc2.pixel * oc2.has_alpha + oc3.pixel * oc3.has_alpha
829  + oc4.pixel * oc4.has_alpha + oc5.pixel * oc5.has_alpha
830  + oc6.pixel * oc6.has_alpha + oc7.pixel * oc7.has_alpha ) / static_cast< double >( number_of_visible_voxels );
831 
832  // 透明物体があった場合は，周りの不透明物体の色で置き換えることでもあれを回避
833  if( !oc0.has_alpha ){ oc0.pixel = oc.pixel; }
834  if( !oc1.has_alpha ){ oc1.pixel = oc.pixel; }
835  if( !oc2.has_alpha ){ oc2.pixel = oc.pixel; }
836  if( !oc3.has_alpha ){ oc3.pixel = oc.pixel; }
837  if( !oc4.has_alpha ){ oc4.pixel = oc.pixel; }
838  if( !oc5.has_alpha ){ oc5.pixel = oc.pixel; }
839  if( !oc6.has_alpha ){ oc6.pixel = oc.pixel; }
840  if( !oc7.has_alpha ){ oc7.pixel = oc.pixel; }
841 
842  oc.pixel = ( oc0.pixel + ( oc3.pixel - oc0.pixel ) * xx ) + ( oc1.pixel - oc0.pixel + ( oc0.pixel - oc1.pixel + oc2.pixel - oc3.pixel ) * xx ) * yy;
843  oc.pixel = oc.pixel + ( ( oc4.pixel + ( oc7.pixel - oc4.pixel ) * xx ) + ( oc5.pixel - oc4.pixel + ( oc4.pixel - oc5.pixel + oc6.pixel - oc7.pixel ) * xx ) * yy - oc.pixel ) * zz;
844  oc.alpha = ( oc0.alpha + ( oc3.alpha - oc0.alpha ) * xx ) + ( oc1.alpha - oc0.alpha + ( oc0.alpha - oc1.alpha + oc2.alpha - oc3.alpha ) * xx ) * yy;
845  oc.alpha = oc.alpha + ( ( oc4.alpha + ( oc7.alpha - oc4.alpha ) * xx ) + ( oc5.alpha - oc4.alpha + ( oc4.alpha - oc5.alpha + oc6.alpha - oc7.alpha ) * xx ) * yy - oc.alpha ) * zz;
846 
847  return( true );
848  }
849  }
850 
851  vector_type normal( difference_type i, difference_type j, difference_type k, double xx, double yy, double zz ) const
852  {
853  const_pointer p0 = &in( i, j, k );
854  const_pointer p1 = p0 + d1;
855  const_pointer p2 = p0 + d2;
856  const_pointer p3 = p0 + d3;
857  const_pointer p4 = p0 + d4;
858  const_pointer p5 = p0 + d5;
859  const_pointer p6 = p0 + d6;
860  const_pointer p7 = p0 + d7;
861 
862  double n0x = p3[ 0 ] - p0[ -_1 ];
863  double n0y = p1[ 0 ] - p0[ -_2 ];
864  double n0z = p4[ 0 ] - p0[ -_3 ];
865  double n1x = p2[ 0 ] - p1[ -_1 ];
866  double n1y = p1[ _2 ] - p0[ 0 ];
867  double n1z = p5[ 0 ] - p1[ -_3 ];
868  double n2x = p2[ _1 ] - p1[ 0 ];
869  double n2y = p2[ _2 ] - p3[ 0 ];
870  double n2z = p6[ 0 ] - p2[ -_3 ];
871  double n3x = p3[ _1 ] - p0[ 0 ];
872  double n3y = p2[ 0 ] - p3[ -_2 ];
873  double n3z = p7[ 0 ] - p3[ -_3 ];
874  double n4x = p7[ 0 ] - p4[ -_1 ];
875  double n4y = p5[ 0 ] - p4[ -_2 ];
876  double n4z = p4[ _3 ] - p0[ 0 ];
877  double n5x = p6[ 0 ] - p5[ -_1 ];
878  double n5y = p5[ _2 ] - p4[ 0 ];
879  double n5z = p5[ _3 ] - p1[ 0 ];
880  double n6x = p6[ _1 ] - p5[ 0 ];
881  double n6y = p6[ _2 ] - p7[ 0 ];
882  double n6z = p6[ _3 ] - p2[ 0 ];
883  double n7x = p7[ _1 ] - p4[ 0 ];
884  double n7y = p6[ 0 ] - p7[ -_2 ];
885  double n7z = p7[ _3 ] - p3[ 0 ];
886 
887  double nx = ( n0x + ( n3x - n0x ) * xx ) + ( n1x - n0x + ( n0x - n1x + n2x - n3x ) * xx ) * yy;
888  nx += ( ( n4x + ( n7x - n4x ) * xx ) + ( n5x - n4x + ( n4x - n5x + n6x - n7x ) * xx ) * yy - nx ) * zz;
889  double ny = ( n0y + ( n3y - n0y ) * xx ) + ( n1y - n0y + ( n0y - n1y + n2y - n3y ) * xx ) * yy;
890  ny += ( ( n4y + ( n7y - n4y ) * xx ) + ( n5y - n4y + ( n4y - n5y + n6y - n7y ) * xx ) * yy - ny ) * zz;
891  double nz = ( n0z + ( n3z - n0z ) * xx ) + ( n1z - n0z + ( n0z - n1z + n2z - n3z ) * xx ) * yy;
892  nz += ( ( n4z + ( n7z - n4z ) * xx ) + ( n5z - n4z + ( n4z - n5z + n6z - n7z ) * xx ) * yy - nz ) * zz;
893 
894  nx /= in.reso1( );
895  ny /= in.reso2( );
896  nz /= in.reso3( );
897  double len = std::sqrt( nx * nx + ny * ny + nz * nz ) + type_limits< double >::tiny( );
898  nx /= len;
899  ny /= len;
900  nz /= len;
901 
902  return( vector_type( nx, ny, nz ) );
903  }
904  };
905 
906 
907  template < class Array1, class Array2, class T >
908  class value_interpolation_with_mark
909  {
910  protected:
911  typedef typename volumerender::parameter::vector_type vector_type;
912  typedef typename Array1::size_type size_type;
913  typedef typename Array1::const_pointer const_pointer;
914  typedef typename Array2::const_pointer const_mk_pointer;
915  typedef typename Array1::difference_type difference_type;
916  typedef typename volumerender::attribute_table< T >::attribute_type attribute_type;
917 
918  public:
919  const Array1 &in;
920  const Array2 &mk;
921  const volumerender::parameter &p;
922  const volumerender::attribute_table< T > &table;
923  const volumerender::attribute_table< T > &mktable;
924  difference_type d0, d1, d2, d3, d4, d5, d6, d7, _1, _2, _3;
925 
926  value_interpolation_with_mark( const Array1 &image, const Array2 &mark, const volumerender::parameter &param, const volumerender::attribute_table< T > &tbl, const volumerender::attribute_table< T > &mktbl )
927  : in( image ), mk( mark ), p( param ), table( tbl ), mktable( mktbl )
928  {
929  difference_type cx = in.width( ) / 2;
930  difference_type cy = in.height( ) / 2;
931  difference_type cz = in.depth( ) / 2;
932  const_pointer ppp = &in( cx, cy, cz );
933  d0 = 0;
934  d1 = &in( cx , cy + 1, cz ) - ppp;
935  d2 = &in( cx + 1, cy + 1, cz ) - ppp;
936  d3 = &in( cx + 1, cy , cz ) - ppp;
937  d4 = &in( cx , cy , cz + 1 ) - ppp;
938  d5 = &in( cx , cy + 1, cz + 1 ) - ppp;
939  d6 = &in( cx + 1, cy + 1, cz + 1 ) - ppp;
940  d7 = &in( cx + 1, cy , cz + 1 ) - ppp;
941  _1 = &in( cx + 1, cy , cz ) - ppp;
942  _2 = &in( cx , cy + 1, cz ) - ppp;
943  _3 = &in( cx , cy , cz + 1 ) - ppp;
944  }
945 
946  bool check( difference_type i, difference_type j, difference_type k ) const
947  {
948  const_pointer p = &in( i, j, k );
949 
950  // この位置における物体が不透明の場合は次のステップへ移行する
951  return( table.has_alpha( p[ d0 ] ) || table.has_alpha( p[ d1 ] ) ||
952  table.has_alpha( p[ d2 ] ) || table.has_alpha( p[ d3 ] ) ||
953  table.has_alpha( p[ d4 ] ) || table.has_alpha( p[ d5 ] ) ||
954  table.has_alpha( p[ d6 ] ) || table.has_alpha( p[ d7 ] ) );
955  }
956 
957  bool render( difference_type i, difference_type j, difference_type k, double xx, double yy, double zz, attribute_type &oc ) const
958  {
959  const_pointer p = &in( i, j, k );
960 
961  // CT値に対応する色と不透明度を取得
962  double ct = ( p[ d0 ] + ( p[ d3 ] - p[ d0 ] ) * xx ) + ( p[ d1 ] - p[ d0 ] + ( p[ d0 ] - p[ d1 ] + p[ d2 ] - p[ d3 ] ) * xx ) * yy;
963  ct += ( ( p[ d4 ] + ( p[ d7 ] - p[ d4 ] ) * xx ) + ( p[ d5 ] - p[ d4 ] + ( p[ d4 ] - p[ d5 ] + p[ d6 ] - p[ d7 ] ) * xx ) * yy - ct ) * zz;
964 
965  oc = table[ volumerender::to_integer( ct ) ];
966 
967  if( !oc.has_alpha )
968  {
969  return( false );
970  }
971  else
972  {
973  // マーキング結果を重ねる
974  const_mk_pointer pm = &mk( i, j, k );
975 
976  // 近傍の8点に対応する色と不透明度を取得
977  attribute_type mc0 = mktable[ volumerender::to_integer( pm[ d0 ] ) ];
978  attribute_type mc1 = mktable[ volumerender::to_integer( pm[ d1 ] ) ];
979  attribute_type mc2 = mktable[ volumerender::to_integer( pm[ d2 ] ) ];
980  attribute_type mc3 = mktable[ volumerender::to_integer( pm[ d3 ] ) ];
981  attribute_type mc4 = mktable[ volumerender::to_integer( pm[ d4 ] ) ];
982  attribute_type mc5 = mktable[ volumerender::to_integer( pm[ d5 ] ) ];
983  attribute_type mc6 = mktable[ volumerender::to_integer( pm[ d6 ] ) ];
984  attribute_type mc7 = mktable[ volumerender::to_integer( pm[ d7 ] ) ];
985 
986  int number_of_visible_voxels = mc0.has_alpha;
987  number_of_visible_voxels += mc1.has_alpha;
988  number_of_visible_voxels += mc2.has_alpha;
989  number_of_visible_voxels += mc3.has_alpha;
990  number_of_visible_voxels += mc4.has_alpha;
991  number_of_visible_voxels += mc5.has_alpha;
992  number_of_visible_voxels += mc6.has_alpha;
993  number_of_visible_voxels += mc7.has_alpha;
994 
995  if( number_of_visible_voxels > 0 )
996  {
997  attribute_type mc;
998  // まず平均的な色を決定する
999  mc.pixel = ( mc0.pixel * mc0.has_alpha + mc1.pixel * mc1.has_alpha
1000  + mc2.pixel * mc2.has_alpha + mc3.pixel * mc3.has_alpha
1001  + mc4.pixel * mc4.has_alpha + mc5.pixel * mc5.has_alpha
1002  + mc6.pixel * mc6.has_alpha + mc7.pixel * mc7.has_alpha ) / static_cast< double >( number_of_visible_voxels );
1003 
1004  // 透明物体があった場合は，周りの不透明物体の色で置き換えることでもあれを回避
1005  if( !mc0.has_alpha ){ mc0.pixel = mc.pixel; }
1006  if( !mc1.has_alpha ){ mc1.pixel = mc.pixel; }
1007  if( !mc2.has_alpha ){ mc2.pixel = mc.pixel; }
1008  if( !mc3.has_alpha ){ mc3.pixel = mc.pixel; }
1009  if( !mc4.has_alpha ){ mc4.pixel = mc.pixel; }
1010  if( !mc5.has_alpha ){ mc5.pixel = mc.pixel; }
1011  if( !mc6.has_alpha ){ mc6.pixel = mc.pixel; }
1012  if( !mc7.has_alpha ){ mc7.pixel = mc.pixel; }
1013 
1014  mc.pixel = ( mc0.pixel + ( mc3.pixel - mc0.pixel ) * xx ) + ( mc1.pixel - mc0.pixel + ( mc0.pixel - mc1.pixel + mc2.pixel - mc3.pixel ) * xx ) * yy;
1015  mc.pixel = mc.pixel + ( ( mc4.pixel + ( mc7.pixel - mc4.pixel ) * xx ) + ( mc5.pixel - mc4.pixel + ( mc4.pixel - mc5.pixel + mc6.pixel - mc7.pixel ) * xx ) * yy - mc.pixel ) * zz;
1016  mc.alpha = ( mc0.alpha + ( mc3.alpha - mc0.alpha ) * xx ) + ( mc1.alpha - mc0.alpha + ( mc0.alpha - mc1.alpha + mc2.alpha - mc3.alpha ) * xx ) * yy;
1017  mc.alpha = mc.alpha + ( ( mc4.alpha + ( mc7.alpha - mc4.alpha ) * xx ) + ( mc5.alpha - mc4.alpha + ( mc4.alpha - mc5.alpha + mc6.alpha - mc7.alpha ) * xx ) * yy - mc.alpha ) * zz;
1018 
1019  // マーキングの色を合成する
1020  oc.pixel = oc.pixel * ( 1.0 - mc.alpha ) + mc.pixel * mc.alpha;
1021  }
1022 
1023  return( true );
1024  }
1025  }
1026 
1027  vector_type normal( difference_type i, difference_type j, difference_type k, double xx, double yy, double zz ) const
1028  {
1029  const_pointer p0 = &in( i, j, k );
1030  const_pointer p1 = p0 + d1;
1031  const_pointer p2 = p0 + d2;
1032  const_pointer p3 = p0 + d3;
1033  const_pointer p4 = p0 + d4;
1034  const_pointer p5 = p0 + d5;
1035  const_pointer p6 = p0 + d6;
1036  const_pointer p7 = p0 + d7;
1037 
1038  double n0x = p3[ 0 ] - p0[ -_1 ];
1039  double n0y = p1[ 0 ] - p0[ -_2 ];
1040  double n0z = p4[ 0 ] - p0[ -_3 ];
1041  double n1x = p2[ 0 ] - p1[ -_1 ];
1042  double n1y = p1[ _2 ] - p0[ 0 ];
1043  double n1z = p5[ 0 ] - p1[ -_3 ];
1044  double n2x = p2[ _1 ] - p1[ 0 ];
1045  double n2y = p2[ _2 ] - p3[ 0 ];
1046  double n2z = p6[ 0 ] - p2[ -_3 ];
1047  double n3x = p3[ _1 ] - p0[ 0 ];
1048  double n3y = p2[ 0 ] - p3[ -_2 ];
1049  double n3z = p7[ 0 ] - p3[ -_3 ];
1050  double n4x = p7[ 0 ] - p4[ -_1 ];
1051  double n4y = p5[ 0 ] - p4[ -_2 ];
1052  double n4z = p4[ _3 ] - p0[ 0 ];
1053  double n5x = p6[ 0 ] - p5[ -_1 ];
1054  double n5y = p5[ _2 ] - p4[ 0 ];
1055  double n5z = p5[ _3 ] - p1[ 0 ];
1056  double n6x = p6[ _1 ] - p5[ 0 ];
1057  double n6y = p6[ _2 ] - p7[ 0 ];
1058  double n6z = p6[ _3 ] - p2[ 0 ];
1059  double n7x = p7[ _1 ] - p4[ 0 ];
1060  double n7y = p6[ 0 ] - p7[ -_2 ];
1061  double n7z = p7[ _3 ] - p3[ 0 ];
1062 
1063  double nx = ( n0x + ( n3x - n0x ) * xx ) + ( n1x - n0x + ( n0x - n1x + n2x - n3x ) * xx ) * yy;
1064  nx += ( ( n4x + ( n7x - n4x ) * xx ) + ( n5x - n4x + ( n4x - n5x + n6x - n7x ) * xx ) * yy - nx ) * zz;
1065  double ny = ( n0y + ( n3y - n0y ) * xx ) + ( n1y - n0y + ( n0y - n1y + n2y - n3y ) * xx ) * yy;
1066  ny += ( ( n4y + ( n7y - n4y ) * xx ) + ( n5y - n4y + ( n4y - n5y + n6y - n7y ) * xx ) * yy - ny ) * zz;
1067  double nz = ( n0z + ( n3z - n0z ) * xx ) + ( n1z - n0z + ( n0z - n1z + n2z - n3z ) * xx ) * yy;
1068  nz += ( ( n4z + ( n7z - n4z ) * xx ) + ( n5z - n4z + ( n4z - n5z + n6z - n7z ) * xx ) * yy - nz ) * zz;
1069 
1070  nx /= in.reso1( );
1071  ny /= in.reso2( );
1072  nz /= in.reso3( );
1073  double len = std::sqrt( nx * nx + ny * ny + nz * nz ) + type_limits< double >::tiny( );
1074  nx /= len;
1075  ny /= len;
1076  nz /= len;
1077 
1078  return( vector_type( nx, ny, nz ) );
1079  }
1080  };
1081 
1082 
1083  template < class Array1, class Array2, class T >
1084  class value_interpolation_and_mark
1085  {
1086  protected:
1087  typedef typename volumerender::parameter::vector_type vector_type;
1088  typedef typename Array1::size_type size_type;
1089  typedef typename Array1::const_pointer const_pointer;
1090  typedef typename Array2::const_pointer const_mk_pointer;
1091  typedef typename Array1::difference_type difference_type;
1092  typedef typename volumerender::attribute_table< T >::attribute_type attribute_type;
1093 
1094  public:
1095  const Array1 &in;
1096  const Array2 &mk;
1097  const volumerender::parameter &p;
1098  const volumerender::attribute_table< T > &table;
1099  const volumerender::attribute_table< T > &mktable;
1100  difference_type d0, d1, d2, d3, d4, d5, d6, d7, _1, _2, _3;
1101 
1102  value_interpolation_and_mark( const Array1 &image, const Array2 &mark, const volumerender::parameter &param, const volumerender::attribute_table< T > &tbl, const volumerender::attribute_table< T > &mktbl )
1103  : in( image ), mk( mark ), p( param ), table( tbl ), mktable( mktbl )
1104  {
1105  difference_type cx = in.width( ) / 2;
1106  difference_type cy = in.height( ) / 2;
1107  difference_type cz = in.depth( ) / 2;
1108  const_pointer ppp = &in( cx, cy, cz );
1109  d0 = 0;
1110  d1 = &in( cx , cy + 1, cz ) - ppp;
1111  d2 = &in( cx + 1, cy + 1, cz ) - ppp;
1112  d3 = &in( cx + 1, cy , cz ) - ppp;
1113  d4 = &in( cx , cy , cz + 1 ) - ppp;
1114  d5 = &in( cx , cy + 1, cz + 1 ) - ppp;
1115  d6 = &in( cx + 1, cy + 1, cz + 1 ) - ppp;
1116  d7 = &in( cx + 1, cy , cz + 1 ) - ppp;
1117  _1 = &in( cx + 1, cy , cz ) - ppp;
1118  _2 = &in( cx , cy + 1, cz ) - ppp;
1119  _3 = &in( cx , cy , cz + 1 ) - ppp;
1120  }
1121 
1122  bool check( difference_type i, difference_type j, difference_type k ) const
1123  {
1124  const_pointer p = &in( i, j, k );
1125  const_mk_pointer pm = &mk( i, j, k );
1126 
1127  // この位置における物体が不透明の場合は次のステップへ移行する
1128  return( ( table.has_alpha( p[ d0 ] ) || table.has_alpha( p[ d1 ] ) || table.has_alpha( p[ d2 ] ) ||
1129  table.has_alpha( p[ d3 ] ) || table.has_alpha( p[ d4 ] ) || table.has_alpha( p[ d5 ] ) ||
1130  table.has_alpha( p[ d6 ] ) || table.has_alpha( p[ d7 ] ) ) &&
1131  ( mktable.has_alpha( pm[ d0 ] ) || mktable.has_alpha( pm[ d1 ] ) || mktable.has_alpha( pm[ d2 ] ) ||
1132  mktable.has_alpha( pm[ d3 ] ) || mktable.has_alpha( pm[ d4 ] ) || mktable.has_alpha( pm[ d5 ] ) ||
1133  mktable.has_alpha( pm[ d6 ] ) || mktable.has_alpha( pm[ d7 ] ) ) );
1134  }
1135 
1136  bool render( difference_type i, difference_type j, difference_type k, double xx, double yy, double zz, attribute_type &oc ) const
1137  {
1138  // マーキング結果を重ねる
1139  const_mk_pointer pm = &mk( i, j, k );
1140 
1141  // 近傍の8点に対応する色と不透明度を取得
1142  const attribute_type &mc0 = mktable[ volumerender::to_integer( pm[ d0 ] ) ];
1143  const attribute_type &mc1 = mktable[ volumerender::to_integer( pm[ d1 ] ) ];
1144  const attribute_type &mc2 = mktable[ volumerender::to_integer( pm[ d2 ] ) ];
1145  const attribute_type &mc3 = mktable[ volumerender::to_integer( pm[ d3 ] ) ];
1146  const attribute_type &mc4 = mktable[ volumerender::to_integer( pm[ d4 ] ) ];
1147  const attribute_type &mc5 = mktable[ volumerender::to_integer( pm[ d5 ] ) ];
1148  const attribute_type &mc6 = mktable[ volumerender::to_integer( pm[ d6 ] ) ];
1149  const attribute_type &mc7 = mktable[ volumerender::to_integer( pm[ d7 ] ) ];
1150 
1151  int number_of_visible_voxels = mc0.has_alpha;
1152  number_of_visible_voxels += mc1.has_alpha;
1153  number_of_visible_voxels += mc2.has_alpha;
1154  number_of_visible_voxels += mc3.has_alpha;
1155  number_of_visible_voxels += mc4.has_alpha;
1156  number_of_visible_voxels += mc5.has_alpha;
1157  number_of_visible_voxels += mc6.has_alpha;
1158  number_of_visible_voxels += mc7.has_alpha;
1159 
1160  if( number_of_visible_voxels > 0 )
1161  {
1162  const_pointer p = &in( i, j, k );
1163 
1164  // CT値に対応する色と不透明度を取得
1165  double ct = ( p[ d0 ] + ( p[ d3 ] - p[ d0 ] ) * xx ) + ( p[ d1 ] - p[ d0 ] + ( p[ d0 ] - p[ d1 ] + p[ d2 ] - p[ d3 ] ) * xx ) * yy;
1166  ct += ( ( p[ d4 ] + ( p[ d7 ] - p[ d4 ] ) * xx ) + ( p[ d5 ] - p[ d4 ] + ( p[ d4 ] - p[ d5 ] + p[ d6 ] - p[ d7 ] ) * xx ) * yy - ct ) * zz;
1167 
1168  oc = table[ volumerender::to_integer( ct ) ];
1169 
1170  // マークデータの透明度を線形補間する
1171  double alpha = ( mc0.alpha + ( mc3.alpha - mc0.alpha ) * xx ) + ( mc1.alpha - mc0.alpha + ( mc0.alpha - mc1.alpha + mc2.alpha - mc3.alpha ) * xx ) * yy;
1172  alpha = alpha + ( ( mc4.alpha + ( mc7.alpha - mc4.alpha ) * xx ) + ( mc5.alpha - mc4.alpha + ( mc4.alpha - mc5.alpha + mc6.alpha - mc7.alpha ) * xx ) * yy - alpha ) * zz;
1173 
1174  // マーキングの色を合成する
1175  oc.alpha *= alpha;
1176 
1177  return( true );
1178  }
1179  else
1180  {
1181  return( false );
1182  }
1183  }
1184 
1185  vector_type normal( difference_type i, difference_type j, difference_type k, double xx, double yy, double zz ) const
1186  {
1187  const_pointer p0 = &in( i, j, k );
1188  const_pointer p1 = p0 + d1;
1189  const_pointer p2 = p0 + d2;
1190  const_pointer p3 = p0 + d3;
1191  const_pointer p4 = p0 + d4;
1192  const_pointer p5 = p0 + d5;
1193  const_pointer p6 = p0 + d6;
1194  const_pointer p7 = p0 + d7;
1195 
1196  double n0x = p3[ 0 ] - p0[ -_1 ];
1197  double n0y = p1[ 0 ] - p0[ -_2 ];
1198  double n0z = p4[ 0 ] - p0[ -_3 ];
1199  double n1x = p2[ 0 ] - p1[ -_1 ];
1200  double n1y = p1[ _2 ] - p0[ 0 ];
1201  double n1z = p5[ 0 ] - p1[ -_3 ];
1202  double n2x = p2[ _1 ] - p1[ 0 ];
1203  double n2y = p2[ _2 ] - p3[ 0 ];
1204  double n2z = p6[ 0 ] - p2[ -_3 ];
1205  double n3x = p3[ _1 ] - p0[ 0 ];
1206  double n3y = p2[ 0 ] - p3[ -_2 ];
1207  double n3z = p7[ 0 ] - p3[ -_3 ];
1208  double n4x = p7[ 0 ] - p4[ -_1 ];
1209  double n4y = p5[ 0 ] - p4[ -_2 ];
1210  double n4z = p4[ _3 ] - p0[ 0 ];
1211  double n5x = p6[ 0 ] - p5[ -_1 ];
1212  double n5y = p5[ _2 ] - p4[ 0 ];
1213  double n5z = p5[ _3 ] - p1[ 0 ];
1214  double n6x = p6[ _1 ] - p5[ 0 ];
1215  double n6y = p6[ _2 ] - p7[ 0 ];
1216  double n6z = p6[ _3 ] - p2[ 0 ];
1217  double n7x = p7[ _1 ] - p4[ 0 ];
1218  double n7y = p6[ 0 ] - p7[ -_2 ];
1219  double n7z = p7[ _3 ] - p3[ 0 ];
1220 
1221  double nx = ( n0x + ( n3x - n0x ) * xx ) + ( n1x - n0x + ( n0x - n1x + n2x - n3x ) * xx ) * yy;
1222  nx += ( ( n4x + ( n7x - n4x ) * xx ) + ( n5x - n4x + ( n4x - n5x + n6x - n7x ) * xx ) * yy - nx ) * zz;
1223  double ny = ( n0y + ( n3y - n0y ) * xx ) + ( n1y - n0y + ( n0y - n1y + n2y - n3y ) * xx ) * yy;
1224  ny += ( ( n4y + ( n7y - n4y ) * xx ) + ( n5y - n4y + ( n4y - n5y + n6y - n7y ) * xx ) * yy - ny ) * zz;
1225  double nz = ( n0z + ( n3z - n0z ) * xx ) + ( n1z - n0z + ( n0z - n1z + n2z - n3z ) * xx ) * yy;
1226  nz += ( ( n4z + ( n7z - n4z ) * xx ) + ( n5z - n4z + ( n4z - n5z + n6z - n7z ) * xx ) * yy - nz ) * zz;
1227 
1228  nx /= in.reso1( );
1229  ny /= in.reso2( );
1230  nz /= in.reso3( );
1231  double len = std::sqrt( nx * nx + ny * ny + nz * nz ) + type_limits< double >::tiny( );
1232  nx /= len;
1233  ny /= len;
1234  nz /= len;
1235 
1236  return( vector_type( nx, ny, nz ) );
1237  }
1238  };
1239 
1240 
1241  template < class Array1, class Array2, class T >
1242  class value_interpolation_or_mark
1243  {
1244  protected:
1245  typedef typename volumerender::parameter::vector_type vector_type;
1246  typedef typename Array1::size_type size_type;
1247  typedef typename Array1::const_pointer const_pointer;
1248  typedef typename Array2::const_pointer const_mk_pointer;
1249  typedef typename Array1::difference_type difference_type;
1250  typedef typename volumerender::attribute_table< T >::attribute_type attribute_type;
1251 
1252  public:
1253  const Array1 &in;
1254  const Array2 &mk;
1255  const volumerender::parameter &p;
1256  const volumerender::attribute_table< T > &table;
1257  const volumerender::attribute_table< T > &mktable;
1258  difference_type d0, d1, d2, d3, d4, d5, d6, d7, _1, _2, _3;
1259 
1260  value_interpolation_or_mark( const Array1 &image, const Array2 &mark, const volumerender::parameter &param, const volumerender::attribute_table< T > &tbl, const volumerender::attribute_table< T > &mktbl )
1261  : in( image ), mk( mark ), p( param ), table( tbl ), mktable( mktbl )
1262  {
1263  difference_type cx = in.width( ) / 2;
1264  difference_type cy = in.height( ) / 2;
1265  difference_type cz = in.depth( ) / 2;
1266  const_pointer ppp = &in( cx, cy, cz );
1267  d0 = 0;
1268  d1 = &in( cx , cy + 1, cz ) - ppp;
1269  d2 = &in( cx + 1, cy + 1, cz ) - ppp;
1270  d3 = &in( cx + 1, cy , cz ) - ppp;
1271  d4 = &in( cx , cy , cz + 1 ) - ppp;
1272  d5 = &in( cx , cy + 1, cz + 1 ) - ppp;
1273  d6 = &in( cx + 1, cy + 1, cz + 1 ) - ppp;
1274  d7 = &in( cx + 1, cy , cz + 1 ) - ppp;
1275  _1 = &in( cx + 1, cy , cz ) - ppp;
1276  _2 = &in( cx , cy + 1, cz ) - ppp;
1277  _3 = &in( cx , cy , cz + 1 ) - ppp;
1278  }
1279 
1280  bool check( difference_type i, difference_type j, difference_type k ) const
1281  {
1282  const_pointer p = &in( i, j, k );
1283  const_mk_pointer pm = &mk( i, j, k );
1284 
1285  // この位置における物体が不透明の場合は次のステップへ移行する
1286  return( table.has_alpha( p[ d0 ] ) || table.has_alpha( p[ d1 ] ) ||
1287  table.has_alpha( p[ d2 ] ) || table.has_alpha( p[ d3 ] ) ||
1288  table.has_alpha( p[ d4 ] ) || table.has_alpha( p[ d5 ] ) ||
1289  table.has_alpha( p[ d6 ] ) || table.has_alpha( p[ d7 ] ) ||
1290  mktable.has_alpha( pm[ d0 ] ) || mktable.has_alpha( pm[ d1 ] ) ||
1291  mktable.has_alpha( pm[ d2 ] ) || mktable.has_alpha( pm[ d3 ] ) ||
1292  mktable.has_alpha( pm[ d4 ] ) || mktable.has_alpha( pm[ d5 ] ) ||
1293  mktable.has_alpha( pm[ d6 ] ) || mktable.has_alpha( pm[ d7 ] ) );
1294  }
1295 
1296  bool render( difference_type i, difference_type j, difference_type k, double xx, double yy, double zz, attribute_type &oc ) const
1297  {
1298  const_pointer p = &in( i, j, k );
1299 
1300  // CT値に対応する色と不透明度を取得
1301  double ct = ( p[ d0 ] + ( p[ d3 ] - p[ d0 ] ) * xx ) + ( p[ d1 ] - p[ d0 ] + ( p[ d0 ] - p[ d1 ] + p[ d2 ] - p[ d3 ] ) * xx ) * yy;
1302  ct += ( ( p[ d4 ] + ( p[ d7 ] - p[ d4 ] ) * xx ) + ( p[ d5 ] - p[ d4 ] + ( p[ d4 ] - p[ d5 ] + p[ d6 ] - p[ d7 ] ) * xx ) * yy - ct ) * zz;
1303 
1304  oc = table[ volumerender::to_integer( ct ) ];
1305 
1306  // マーキング結果を重ねる
1307  const_mk_pointer pm = &mk( i, j, k );
1308 
1309  // 近傍の8点に対応する色と不透明度を取得
1310  attribute_type mc0 = mktable[ volumerender::to_integer( pm[ d0 ] ) ];
1311  attribute_type mc1 = mktable[ volumerender::to_integer( pm[ d1 ] ) ];
1312  attribute_type mc2 = mktable[ volumerender::to_integer( pm[ d2 ] ) ];
1313  attribute_type mc3 = mktable[ volumerender::to_integer( pm[ d3 ] ) ];
1314  attribute_type mc4 = mktable[ volumerender::to_integer( pm[ d4 ] ) ];
1315  attribute_type mc5 = mktable[ volumerender::to_integer( pm[ d5 ] ) ];
1316  attribute_type mc6 = mktable[ volumerender::to_integer( pm[ d6 ] ) ];
1317  attribute_type mc7 = mktable[ volumerender::to_integer( pm[ d7 ] ) ];
1318 
1319  int number_of_visible_voxels = mc0.has_alpha;
1320  number_of_visible_voxels += mc1.has_alpha;
1321  number_of_visible_voxels += mc2.has_alpha;
1322  number_of_visible_voxels += mc3.has_alpha;
1323  number_of_visible_voxels += mc4.has_alpha;
1324  number_of_visible_voxels += mc5.has_alpha;
1325  number_of_visible_voxels += mc6.has_alpha;
1326  number_of_visible_voxels += mc7.has_alpha;
1327 
1328  if( number_of_visible_voxels > 0 )
1329  {
1330 
1331  attribute_type mc;
1332  // まず平均的な色を決定する
1333  mc.pixel = ( mc0.pixel * mc0.has_alpha + mc1.pixel * mc1.has_alpha
1334  + mc2.pixel * mc2.has_alpha + mc3.pixel * mc3.has_alpha
1335  + mc4.pixel * mc4.has_alpha + mc5.pixel * mc5.has_alpha
1336  + mc6.pixel * mc6.has_alpha + mc7.pixel * mc7.has_alpha ) / static_cast< double >( number_of_visible_voxels );
1337 
1338  // 透明物体があった場合は，周りの不透明物体の色で置き換えることでもあれを回避
1339  if( !mc0.has_alpha ){ mc0.pixel = mc.pixel; }
1340  if( !mc1.has_alpha ){ mc1.pixel = mc.pixel; }
1341  if( !mc2.has_alpha ){ mc2.pixel = mc.pixel; }
1342  if( !mc3.has_alpha ){ mc3.pixel = mc.pixel; }
1343  if( !mc4.has_alpha ){ mc4.pixel = mc.pixel; }
1344  if( !mc5.has_alpha ){ mc5.pixel = mc.pixel; }
1345  if( !mc6.has_alpha ){ mc6.pixel = mc.pixel; }
1346  if( !mc7.has_alpha ){ mc7.pixel = mc.pixel; }
1347 
1348  mc.pixel = ( mc0.pixel + ( mc3.pixel - mc0.pixel ) * xx ) + ( mc1.pixel - mc0.pixel + ( mc0.pixel - mc1.pixel + mc2.pixel - mc3.pixel ) * xx ) * yy;
1349  mc.pixel = mc.pixel + ( ( mc4.pixel + ( mc7.pixel - mc4.pixel ) * xx ) + ( mc5.pixel - mc4.pixel + ( mc4.pixel - mc5.pixel + mc6.pixel - mc7.pixel ) * xx ) * yy - mc.pixel ) * zz;
1350  mc.alpha = ( mc0.alpha + ( mc3.alpha - mc0.alpha ) * xx ) + ( mc1.alpha - mc0.alpha + ( mc0.alpha - mc1.alpha + mc2.alpha - mc3.alpha ) * xx ) * yy;
1351  mc.alpha = mc.alpha + ( ( mc4.alpha + ( mc7.alpha - mc4.alpha ) * xx ) + ( mc5.alpha - mc4.alpha + ( mc4.alpha - mc5.alpha + mc6.alpha - mc7.alpha ) * xx ) * yy - mc.alpha ) * zz;
1352 
1353  // マーキングの色を合成する
1354  oc.pixel = oc.pixel * ( 1.0 - mc.alpha ) + mc.pixel * mc.alpha;
1355  oc.alpha = ( oc.alpha + mc.alpha ) / 2.0;
1356 
1357  return( true );
1358  }
1359  else
1360  {
1361  return( oc.has_alpha );
1362  }
1363  }
1364 
1365  template < class Pointer >
1366  vector_type normal( Pointer p, double xx, double yy, double zz ) const
1367  {
1368  const_pointer p0 = p;
1369  const_pointer p1 = p0 + d1;
1370  const_pointer p2 = p0 + d2;
1371  const_pointer p3 = p0 + d3;
1372  const_pointer p4 = p0 + d4;
1373  const_pointer p5 = p0 + d5;
1374  const_pointer p6 = p0 + d6;
1375  const_pointer p7 = p0 + d7;
1376 
1377  double n0x = p3[ 0 ] - p0[ -_1 ];
1378  double n0y = p1[ 0 ] - p0[ -_2 ];
1379  double n0z = p4[ 0 ] - p0[ -_3 ];
1380  double n1x = p2[ 0 ] - p1[ -_1 ];
1381  double n1y = p1[ _2 ] - p0[ 0 ];
1382  double n1z = p5[ 0 ] - p1[ -_3 ];
1383  double n2x = p2[ _1 ] - p1[ 0 ];
1384  double n2y = p2[ _2 ] - p3[ 0 ];
1385  double n2z = p6[ 0 ] - p2[ -_3 ];
1386  double n3x = p3[ _1 ] - p0[ 0 ];
1387  double n3y = p2[ 0 ] - p3[ -_2 ];
1388  double n3z = p7[ 0 ] - p3[ -_3 ];
1389  double n4x = p7[ 0 ] - p4[ -_1 ];
1390  double n4y = p5[ 0 ] - p4[ -_2 ];
1391  double n4z = p4[ _3 ] - p0[ 0 ];
1392  double n5x = p6[ 0 ] - p5[ -_1 ];
1393  double n5y = p5[ _2 ] - p4[ 0 ];
1394  double n5z = p5[ _3 ] - p1[ 0 ];
1395  double n6x = p6[ _1 ] - p5[ 0 ];
1396  double n6y = p6[ _2 ] - p7[ 0 ];
1397  double n6z = p6[ _3 ] - p2[ 0 ];
1398  double n7x = p7[ _1 ] - p4[ 0 ];
1399  double n7y = p6[ 0 ] - p7[ -_2 ];
1400  double n7z = p7[ _3 ] - p3[ 0 ];
1401 
1402  double nx = ( n0x + ( n3x - n0x ) * xx ) + ( n1x - n0x + ( n0x - n1x + n2x - n3x ) * xx ) * yy;
1403  nx += ( ( n4x + ( n7x - n4x ) * xx ) + ( n5x - n4x + ( n4x - n5x + n6x - n7x ) * xx ) * yy - nx ) * zz;
1404  double ny = ( n0y + ( n3y - n0y ) * xx ) + ( n1y - n0y + ( n0y - n1y + n2y - n3y ) * xx ) * yy;
1405  ny += ( ( n4y + ( n7y - n4y ) * xx ) + ( n5y - n4y + ( n4y - n5y + n6y - n7y ) * xx ) * yy - ny ) * zz;
1406  double nz = ( n0z + ( n3z - n0z ) * xx ) + ( n1z - n0z + ( n0z - n1z + n2z - n3z ) * xx ) * yy;
1407  nz += ( ( n4z + ( n7z - n4z ) * xx ) + ( n5z - n4z + ( n4z - n5z + n6z - n7z ) * xx ) * yy - nz ) * zz;
1408 
1409  nx /= in.reso1( );
1410  ny /= in.reso2( );
1411  nz /= in.reso3( );
1412  double len = std::sqrt( nx * nx + ny * ny + nz * nz ) + type_limits< double >::tiny( );
1413  nx /= len;
1414  ny /= len;
1415  nz /= len;
1416 
1417  return( vector_type( nx, ny, nz ) );
1418  }
1419 
1420 
1421  vector_type normal( difference_type i, difference_type j, difference_type k, double xx, double yy, double zz ) const
1422  {
1423  const_pointer p = &in( i, j, k );
1424  const_mk_pointer pm = &mk( i, j, k );
1425 
1426  // 近傍の8点に対応する色と不透明度を取得
1427  const attribute_type &oc0 = table[ volumerender::to_integer( p[ d0 ] ) ];
1428  const attribute_type &oc1 = table[ volumerender::to_integer( p[ d1 ] ) ];
1429  const attribute_type &oc2 = table[ volumerender::to_integer( p[ d2 ] ) ];
1430  const attribute_type &oc3 = table[ volumerender::to_integer( p[ d3 ] ) ];
1431  const attribute_type &oc4 = table[ volumerender::to_integer( p[ d4 ] ) ];
1432  const attribute_type &oc5 = table[ volumerender::to_integer( p[ d5 ] ) ];
1433  const attribute_type &oc6 = table[ volumerender::to_integer( p[ d6 ] ) ];
1434  const attribute_type &oc7 = table[ volumerender::to_integer( p[ d7 ] ) ];
1435 
1436  const attribute_type &mc0 = mktable[ volumerender::to_integer( pm[ d0 ] ) ];
1437  const attribute_type &mc1 = mktable[ volumerender::to_integer( pm[ d1 ] ) ];
1438  const attribute_type &mc2 = mktable[ volumerender::to_integer( pm[ d2 ] ) ];
1439  const attribute_type &mc3 = mktable[ volumerender::to_integer( pm[ d3 ] ) ];
1440  const attribute_type &mc4 = mktable[ volumerender::to_integer( pm[ d4 ] ) ];
1441  const attribute_type &mc5 = mktable[ volumerender::to_integer( pm[ d5 ] ) ];
1442  const attribute_type &mc6 = mktable[ volumerender::to_integer( pm[ d6 ] ) ];
1443  const attribute_type &mc7 = mktable[ volumerender::to_integer( pm[ d7 ] ) ];
1444 
1445  double oc_alpha = oc0.alpha + oc1.alpha + oc2.alpha + oc3.alpha + oc4.alpha + oc5.alpha + oc6.alpha + oc7.alpha;
1446  double mc_alpha = mc0.alpha + mc1.alpha + mc2.alpha + mc3.alpha + mc4.alpha + mc5.alpha + mc6.alpha + mc7.alpha;
1447  if( oc_alpha > mc_alpha )
1448  {
1449  return( normal( p, xx, yy, zz ) );
1450  }
1451  else
1452  {
1453  return( normal( pm, xx, yy, zz ) );
1454  }
1455  }
1456  };
1457 }
1458 
1459 
1460 // ボリュームレンダリング
1461 namespace __volumerendering_specialized__
1462 {
1463  // CT値補間タイプのレンダリングに特化したボリュームレンダリングエンジン
1464  template < class Array1, class Array2, class DepthMap, class T >
1465  bool volumerendering( const Array1 &in, Array2 &out, const DepthMap &depth_map, const volumerender::parameter &param, const volumerender::attribute_table< T > &volrtable, typename Array1::size_type thread_id, typename Array1::size_type thread_num )
1466  {
1467  typedef typename volumerender::parameter::vector_type vector_type;
1468  typedef typename volumerender::attribute_table< T >::attribute_type attribute_type;
1469  typedef typename volumerender::attribute_table< T >::pixel_type pixel_type;
1470  typedef typename Array1::size_type size_type;
1471  typedef typename Array1::difference_type difference_type;
1472  typedef typename Array1::value_type value_type;
1473  typedef typename Array1::const_pointer const_pointer;
1474  typedef typename Array2::value_type out_value_type;
1475 
1476  vector_type offset = param.offset;
1477  vector_type pos = param.pos - offset;
1478  const volumerender::boundingbox *box = param.box;
1479  double fovy = param.fovy;
1480  double ambient_ratio = param.ambient_ratio;
1481  double diffuse_ratio = param.diffuse_ratio;
1482  double specular = param.specular;
1483  bool bSpecular = specular > 0.0;
1484  double lightAtten = param.light_attenuation;
1485  bool bLightAtten = lightAtten > 0.0;
1486  double sampling_step = param.sampling_step;
1487  double termination = param.termination;
1488  bool bperspective = param.perspective_view;
1489 
1490  const size_type w = in.width( );
1491  const size_type h = in.height( );
1492  const size_type d = in.depth( );
1493 
1494  const size_type image_width = out.width( );
1495  const size_type image_height = out.height( );
1496 
1497  const attribute_type *table = &volrtable[ 0 ];
1498 
1499  const pixel_type background = _pixel_converter_< pixel_type >::convert_to( param.background_R, param.background_G, param.background_B );
1500 
1501  // 高速にアドレス計算を行うためのポインタの差分
1502  difference_type d0, d1, d2, d3, d4, d5, d6, d7, _1, _2, _3;
1503  {
1504  difference_type cx = in.width( ) / 2;
1505  difference_type cy = in.height( ) / 2;
1506  difference_type cz = in.depth( ) / 2;
1507  const_pointer ppp = &in( cx, cy, cz );
1508  d0 = 0;
1509  d1 = &in( cx , cy + 1, cz ) - ppp;
1510  d2 = &in( cx + 1, cy + 1, cz ) - ppp;
1511  d3 = &in( cx + 1, cy , cz ) - ppp;
1512  d4 = &in( cx , cy , cz + 1 ) - ppp;
1513  d5 = &in( cx , cy + 1, cz + 1 ) - ppp;
1514  d6 = &in( cx + 1, cy + 1, cz + 1 ) - ppp;
1515  d7 = &in( cx + 1, cy , cz + 1 ) - ppp;
1516  _1 = &in( cx + 1, cy , cz ) - ppp;
1517  _2 = &in( cx , cy + 1, cz ) - ppp;
1518  _3 = &in( cx , cy , cz + 1 ) - ppp;
1519  }
1520 
1521  // スライス座標系の実寸をワールドと考える
1522  vector_type casting_start, casting_end;
1523 
1524  const double pai = 3.1415926535897932384626433832795;
1525  double focal = ( static_cast< double >( image_height ) / 2.0 ) / std::tan( fovy * pai / 180.0 / 2.0 );
1526 
1527  double cx = static_cast< double >( image_width ) / 2.0;
1528  double cy = static_cast< double >( image_height ) / 2.0;
1529  double ax = in.reso1( );
1530  double ay = in.reso2( );
1531  double az = in.reso3( );
1532  double _1_ax = 1.0 / ax;
1533  double _1_ay = 1.0 / ay;
1534  double _1_az = 1.0 / az;
1535 
1536  double masp = ax < ay ? ax : ay;
1537  masp = masp < az ? masp : az;
1538 
1539  vector_type eZ = -param.dir.unit( ); // カメラ座標系のZ軸方向ベクトル
1540  vector_type eY = param.up.unit( ); // カメラ座標系のY軸方向ベクトル
1541  vector_type eX = ( eY * eZ ).unit( ); // カメラ座標系のX軸方向ベクトル
1542 
1543  if( param.mirror_view )
1544  {
1545  // 鏡に映ったように描画する
1546  // データの軸等に反転が必要な描画を行う際に利用する
1547  eX = -eX;
1548  }
1549 
1550  if( out.reso1( ) < out.reso2( ) )
1551  {
1552  eX *= out.reso1( ) / out.reso2( );
1553  }
1554  else
1555  {
1556  eY *= out.reso2( ) / out.reso1( );
1557  focal *= out.reso2( ) / out.reso1( );
1558  }
1559 
1560  double max_distance = pos.length( ) + std::sqrt( static_cast< double >( w * w + h * h + d * d ) );
1561 
1562  for( size_type j = thread_id ; j < image_height ; j += thread_num )
1563  {
1564  for( size_type i = 0 ; i < image_width ; i++ )
1565  {
1566  // 投影面上の点をカメラ座標系に変換
1567  vector_type Pos( static_cast< double >( i ) - cx, cy - static_cast< double >( j ), -focal );
1568 
1569  // レイ方向をカメラ座標系からワールド座標系に変換
1570  vector_type light;
1571  if( bperspective )
1572  {
1573  light = ( eX * Pos.x + eY * Pos.y + eZ * Pos.z ).unit( );
1574  }
1575  else
1576  {
1577  pos = param.pos - offset + eX * Pos.x + eY * Pos.y;
1578  light = -eZ;
1579  }
1580 
1581  pixel_type add_intensity( 0 );
1582  double add_opacity = 1;
1583 
1584  casting_start = pos;
1585  casting_end = pos + light * max_distance;
1586  vector_type normal;
1587 
1588  // 物体との衝突判定
1589  if( volumerender::check_intersection( casting_start, casting_end, box[ 0 ], normal )
1590  && volumerender::check_intersection( casting_start, casting_end, box[ 1 ], normal )
1591  && volumerender::check_intersection( casting_start, casting_end, box[ 2 ], normal )
1592  && volumerender::check_intersection( casting_start, casting_end, box[ 3 ], normal )
1593  && volumerender::check_intersection( casting_start, casting_end, box[ 4 ], normal )
1594  && volumerender::check_intersection( casting_start, casting_end, box[ 5 ], normal ) )
1595  {
1596  // 光の減衰を実現するために、カメラからの距離を測る
1597  Pos.x = ( pos.x + offset.x ) * _1_ax;
1598  Pos.y = ( pos.y + offset.y ) * _1_ay;
1599  Pos.z = ( pos.z + offset.z ) * _1_az;
1600 
1601  // ワールド座標系からスライス座標系に変換する
1602  // 以降は、全てスライス座標系で計算する
1603  casting_start.x = ( casting_start.x + offset.x ) * _1_ax;
1604  casting_start.y = ( casting_start.y + offset.y ) * _1_ay;
1605  casting_start.z = ( casting_start.z + offset.z ) * _1_az;
1606  casting_end.x = ( casting_end.x + offset.x ) * _1_ax;
1607  casting_end.y = ( casting_end.y + offset.y ) * _1_ay;
1608  casting_end.z = ( casting_end.z + offset.z ) * _1_az;
1609 
1610  vector_type spos = casting_start;
1611  vector_type ray = ( casting_end - casting_start ).unit( );
1612 
1613  // 光の減衰の距離を実測に直すためのパラメータ
1614  double dlen = vector_type( ray.x * ax, ray.y * ay, ray.z * az ).length( );
1615 
1616  // 直方体画素の画像上では方向によってサンプリング間隔が変わってしまう問題に対応
1617  double ray_sampling_step = sampling_step * masp / dlen;
1618 
1619  vector_type ray_step = ray * ray_sampling_step;
1620 
1621  double n = ( casting_end - casting_start ).length( );
1622  double l = 0, of = ( Pos - casting_start ).length( );
1623 
1624  while( l < n )
1625  {
1626  difference_type si = volumerender::to_integer( spos.x );
1627  difference_type sj = volumerender::to_integer( spos.y );
1628  difference_type sk = volumerender::to_integer( spos.z );
1629 
1630  // 判定を行う範囲の先頭ポインタを取得
1631  const_pointer p = &in( si, sj, sk );
1632 
1633  // この位置における物体が不透明の場合は次のステップへ移行する
1634  if( table[ p[ d0 ] ].has_alpha || table[ p[ d1 ] ].has_alpha ||
1635  table[ p[ d2 ] ].has_alpha || table[ p[ d3 ] ].has_alpha ||
1636  table[ p[ d4 ] ].has_alpha || table[ p[ d5 ] ].has_alpha ||
1637  table[ p[ d6 ] ].has_alpha || table[ p[ d7 ] ].has_alpha )
1638  {
1639  if( l > 0 )
1640  {
1641  spos.x -= ray.x;
1642  spos.y -= ray.y;
1643  spos.z -= ray.z;
1644  l -= 1.0;
1645  }
1646  break;
1647  }
1648 
1649  double current_step = depth_map( si, sj, sk );
1650  l += current_step;
1651  spos.x += ray.x * current_step;
1652  spos.y += ray.y * current_step;
1653  spos.z += ray.z * current_step;
1654 
1655  while( l < n )
1656  {
1657  si = volumerender::to_integer( spos.x );
1658  sj = volumerender::to_integer( spos.y );
1659  sk = volumerender::to_integer( spos.z );
1660 
1661  current_step = depth_map( si, sj, sk );
1662 
1663  if( current_step <= 2.0 )
1664  {
1665  break;
1666  }
1667 
1668  l += current_step;
1669  spos.x += ray.x * current_step;
1670  spos.y += ray.y * current_step;
1671  spos.z += ray.z * current_step;
1672  }
1673  }
1674 
1675  // 端まで到達した場合は何もしない
1676  if( l >= n )
1677  {
1678  out( i, j ) = static_cast< out_value_type >( mist::limits_0_255( background ) );
1679  continue;
1680  }
1681 
1682  double nct[ 8 ] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
1683  double ndx[ 8 ] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
1684  double ndy[ 8 ] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
1685  double ndz[ 8 ] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
1686  const_pointer op = NULL;
1687 
1688  if( l == 0.0 )
1689  {
1690  difference_type si = volumerender::to_integer( spos.x );
1691  difference_type sj = volumerender::to_integer( spos.y );
1692  difference_type sk = volumerender::to_integer( spos.z );
1693 
1694  double xx = spos.x - si;
1695  double yy = spos.y - sj;
1696  double zz = spos.z - sk;
1697 
1698  const_pointer p = &in( si, sj, sk );
1699 
1700  nct[ 0 ] = p[ d0 ];
1701  nct[ 1 ] = p[ d3 ] - p[ d0 ];
1702  nct[ 2 ] = p[ d1 ] - p[ d0 ];
1703  nct[ 3 ] = p[ d2 ] - p[ d3 ] - nct[ 2 ];
1704  nct[ 4 ] = p[ d4 ];
1705  nct[ 5 ] = p[ d7 ] - p[ d4 ];
1706  nct[ 6 ] = p[ d5 ] - p[ d4 ];
1707  nct[ 7 ] = p[ d6 ] - p[ d7 ] - nct[ 6 ];
1708 
1709  // CT値に対応する色と不透明度を取得
1710  double ct = ( nct[ 0 ] + nct[ 1 ] * xx ) + ( nct[ 2 ] + nct[ 3 ] * xx ) * yy;
1711  ct += ( ( nct[ 4 ] + nct[ 5 ] * xx ) + ( nct[ 6 ] + nct[ 7 ] * xx ) * yy - ct ) * zz;
1712 
1713  const attribute_type &oc = table[ volumerender::to_integer( ct ) ];
1714 
1715  if( oc.has_alpha )
1716  {
1717  const_pointer p0 = p;
1718  const_pointer p1 = p0 + d1;
1719  const_pointer p2 = p0 + d2;
1720  const_pointer p3 = p0 + d3;
1721  const_pointer p4 = p0 + d4;
1722  const_pointer p5 = p0 + d5;
1723  const_pointer p6 = p0 + d6;
1724  const_pointer p7 = p0 + d7;
1725 
1726  double n0x = p3[ 0 ] - p0[ -_1 ];
1727  double n0y = p1[ 0 ] - p0[ -_2 ];
1728  double n0z = p4[ 0 ] - p0[ -_3 ];
1729  double n1x = p2[ 0 ] - p1[ -_1 ];
1730  double n1y = p1[ _2 ] - p0[ 0 ];
1731  double n1z = p5[ 0 ] - p1[ -_3 ];
1732  double n2x = p2[ _1 ] - p1[ 0 ];
1733  double n2y = p2[ _2 ] - p3[ 0 ];
1734  double n2z = p6[ 0 ] - p2[ -_3 ];
1735  double n3x = p3[ _1 ] - p0[ 0 ];
1736  double n3y = p2[ 0 ] - p3[ -_2 ];
1737  double n3z = p7[ 0 ] - p3[ -_3 ];
1738  double n4x = p7[ 0 ] - p4[ -_1 ];
1739  double n4y = p5[ 0 ] - p4[ -_2 ];
1740  double n4z = p4[ _3 ] - p0[ 0 ];
1741  double n5x = p6[ 0 ] - p5[ -_1 ];
1742  double n5y = p5[ _2 ] - p4[ 0 ];
1743  double n5z = p5[ _3 ] - p1[ 0 ];
1744  double n6x = p6[ _1 ] - p5[ 0 ];
1745  double n6y = p6[ _2 ] - p7[ 0 ];
1746  double n6z = p6[ _3 ] - p2[ 0 ];
1747  double n7x = p7[ _1 ] - p4[ 0 ];
1748  double n7y = p6[ 0 ] - p7[ -_2 ];
1749  double n7z = p7[ _3 ] - p3[ 0 ];
1750 
1751  ndx[ 0 ] = n0x;
1752  ndx[ 1 ] = n3x - n0x;
1753  ndx[ 2 ] = n1x - n0x;
1754  ndx[ 3 ] = n2x - n3x - ndx[ 2 ];
1755  ndx[ 4 ] = n4x;
1756  ndx[ 5 ] = n7x - n4x;
1757  ndx[ 6 ] = n5x - n4x;
1758  ndx[ 7 ] = n6x - n7x - ndx[ 6 ];
1759 
1760  ndy[ 0 ] = n0y;
1761  ndy[ 1 ] = n3y - n0y;
1762  ndy[ 2 ] = n1y - n0y;
1763  ndy[ 3 ] = n2y - n3y - ndy[ 2 ];
1764  ndy[ 4 ] = n4y;
1765  ndy[ 5 ] = n7y - n4y;
1766  ndy[ 6 ] = n5y - n4y;
1767  ndy[ 7 ] = n6y - n7y - ndy[ 6 ];
1768 
1769  ndz[ 0 ] = n0z;
1770  ndz[ 1 ] = n3z - n0z;
1771  ndz[ 2 ] = n1z - n0z;
1772  ndz[ 3 ] = n2z - n3z - ndz[ 2 ];
1773  ndz[ 4 ] = n4z;
1774  ndz[ 5 ] = n7z - n4z;
1775  ndz[ 6 ] = n5z - n4z;
1776  ndz[ 7 ] = n6z - n7z - ndz[ 6 ];
1777 
1778  op = p;
1779 
1780  double nx = ( ndx[ 0 ] + ndx[ 1 ] * xx ) + ( ndx[ 2 ] + ndx[ 3 ] * xx ) * yy;
1781  nx += ( ( ndx[ 4 ] + ndx[ 5 ] * xx ) + ( ndx[ 6 ] + ndx[ 7 ] * xx ) * yy - nx ) * zz;
1782  double ny = ( ndy[ 0 ] + ndy[ 1 ] * xx ) + ( ndy[ 2 ] + ndy[ 3 ] * xx ) * yy;
1783  ny += ( ( ndy[ 4 ] + ndy[ 5 ] * xx ) + ( ndy[ 6 ] + ndy[ 7 ] * xx ) * yy - ny ) * zz;
1784  double nz = ( ndz[ 0 ] + ndz[ 1 ] * xx ) + ( ndz[ 2 ] + ndz[ 3 ] * xx ) * yy;
1785  nz += ( ( ndz[ 4 ] + ndz[ 5 ] * xx ) + ( ndz[ 6 ] + ndz[ 7 ] * xx ) * yy - nz ) * zz;
1786 
1787  // 切断面の法線を反映させる
1788  double _1_len = 0.5 / ( std::sqrt( nx * nx + ny * ny + nz * nz ) + type_limits< double >::tiny( ) );
1789  nx = ( nx * _1_len + normal.x ) * _1_ax;
1790  ny = ( ny * _1_len + normal.y ) * _1_ay;
1791  nz = ( nz * _1_len + normal.z ) * _1_az;
1792 
1793  // 法線が反転している場合への対応
1794  double c = light.x * nx + light.y * ny + light.z * nz;
1795  c = std::sqrt( ( c * c ) / ( nx * nx + ny * ny + nz * nz + type_limits< double >::tiny( ) ) );
1796 
1797  double spec = 0.0;
1798  if( bSpecular )
1799  {
1800  spec = 2.0 * c * c - 1.0;
1801 
1802  if( spec <= 0.0 )
1803  {
1804  spec = 0;
1805  }
1806  else
1807  {
1808  spec *= spec; // 2乗
1809  spec *= spec; // 4乗
1810  spec *= spec; // 8乗
1811  spec *= spec; // 16乗
1812  spec *= spec; // 32乗
1813  spec *= spec; // 64乗
1814  //spec *= spec; // 128乗
1815  spec *= specular * 255.0;
1816  }
1817  }
1818 
1819  double lAtten = 1.0;
1820  if( bLightAtten )
1821  {
1822  double len = ( l + of ) * dlen;
1823  lAtten /= 1.0 + lightAtten * ( len * len );
1824  }
1825 
1826  double alpha = oc.alpha * sampling_step;
1827  add_intensity += alpha * add_opacity * ( oc.pixel * ( c * diffuse_ratio + ambient_ratio ) + spec ) * lAtten;
1828  add_opacity *= ( 1.0 - alpha );
1829 
1830  // 画素がレンダリング結果に与える影響がしきい値以下になった場合は終了
1831  if( add_opacity < termination )
1832  {
1833  out( i, j ) = static_cast< out_value_type >( mist::limits_0_255( add_intensity * ( 1.0 - add_opacity ) + background * add_opacity ) );
1834  continue;
1835  }
1836 
1837  spos.x += ray_step.x;
1838  spos.y += ray_step.y;
1839  spos.z += ray_step.z;
1840  l += ray_sampling_step;
1841  }
1842  else
1843  {
1844  // この位置における物体が透明の場合は次のステップへ移行する
1845  spos += ray_step;
1846  l += ray_sampling_step;
1847 
1848  double ol = l;
1849  while( l < n )
1850  {
1851  difference_type si = volumerender::to_integer( spos.x );
1852  difference_type sj = volumerender::to_integer( spos.y );
1853  difference_type sk = volumerender::to_integer( spos.z );
1854 
1855  const_pointer p = &in( si, sj, sk );
1856 
1857  // この位置における物体が不透明の場合は次のステップへ移行する
1858  if( table[ p[ d0 ] ].has_alpha || table[ p[ d1 ] ].has_alpha ||
1859  table[ p[ d2 ] ].has_alpha || table[ p[ d3 ] ].has_alpha ||
1860  table[ p[ d4 ] ].has_alpha || table[ p[ d5 ] ].has_alpha ||
1861  table[ p[ d6 ] ].has_alpha || table[ p[ d7 ] ].has_alpha )
1862  {
1863  if( l > ol )
1864  {
1865  spos.x -= ray.x;
1866  spos.y -= ray.y;
1867  spos.z -= ray.z;
1868  l -= 1.0;
1869  }
1870  break;
1871  }
1872 
1873  double current_step = depth_map( si, sj, sk );
1874  l += current_step;
1875  spos.x += ray.x * current_step;
1876  spos.y += ray.y * current_step;
1877  spos.z += ray.z * current_step;
1878 
1879  while( l < n )
1880  {
1881  si = volumerender::to_integer( spos.x );
1882  sj = volumerender::to_integer( spos.y );
1883  sk = volumerender::to_integer( spos.z );
1884 
1885  current_step = depth_map( si, sj, sk );
1886 
1887  if( current_step <= 2.0 )
1888  {
1889  break;
1890  }
1891 
1892  l += current_step;
1893  spos.x += ray.x * current_step;
1894  spos.y += ray.y * current_step;
1895  spos.z += ray.z * current_step;
1896  }
1897  }
1898  }
1899  }
1900 
1901  while( l < n )
1902  {
1903  difference_type si = volumerender::to_integer( spos.x );
1904  difference_type sj = volumerender::to_integer( spos.y );
1905  difference_type sk = volumerender::to_integer( spos.z );
1906 
1907  double xx = spos.x - si;
1908  double yy = spos.y - sj;
1909  double zz = spos.z - sk;
1910 
1911  const_pointer p = &in( si, sj, sk );
1912 
1913  bool need_update = p != op;
1914 
1915  if( need_update )
1916  {
1917  nct[ 0 ] = p[ d0 ];
1918  nct[ 1 ] = p[ d3 ] - p[ d0 ];
1919  nct[ 2 ] = p[ d1 ] - p[ d0 ];
1920  nct[ 3 ] = p[ d2 ] - p[ d3 ] - nct[ 2 ];
1921  nct[ 4 ] = p[ d4 ];
1922  nct[ 5 ] = p[ d7 ] - p[ d4 ];
1923  nct[ 6 ] = p[ d5 ] - p[ d4 ];
1924  nct[ 7 ] = p[ d6 ] - p[ d7 ] - nct[ 6 ];
1925  }
1926 
1927  // CT値に対応する色と不透明度を取得
1928  double ct = ( nct[ 0 ] + nct[ 1 ] * xx ) + ( nct[ 2 ] + nct[ 3 ] * xx ) * yy;
1929  ct += ( ( nct[ 4 ] + nct[ 5 ] * xx ) + ( nct[ 6 ] + nct[ 7 ] * xx ) * yy - ct ) * zz;
1930 
1931  const attribute_type &oc = table[ volumerender::to_integer( ct ) ];
1932 
1933  if( oc.has_alpha )
1934  {
1935  if( need_update )
1936  {
1937  const_pointer p0 = p;
1938  const_pointer p1 = p0 + d1;
1939  const_pointer p2 = p0 + d2;
1940  const_pointer p3 = p0 + d3;
1941  const_pointer p4 = p0 + d4;
1942  const_pointer p5 = p0 + d5;
1943  const_pointer p6 = p0 + d6;
1944  const_pointer p7 = p0 + d7;
1945 
1946  double n0x = p3[ 0 ] - p0[ -_1 ];
1947  double n0y = p1[ 0 ] - p0[ -_2 ];
1948  double n0z = p4[ 0 ] - p0[ -_3 ];
1949  double n1x = p2[ 0 ] - p1[ -_1 ];
1950  double n1y = p1[ _2 ] - p0[ 0 ];
1951  double n1z = p5[ 0 ] - p1[ -_3 ];
1952  double n2x = p2[ _1 ] - p1[ 0 ];
1953  double n2y = p2[ _2 ] - p3[ 0 ];
1954  double n2z = p6[ 0 ] - p2[ -_3 ];
1955  double n3x = p3[ _1 ] - p0[ 0 ];
1956  double n3y = p2[ 0 ] - p3[ -_2 ];
1957  double n3z = p7[ 0 ] - p3[ -_3 ];
1958  double n4x = p7[ 0 ] - p4[ -_1 ];
1959  double n4y = p5[ 0 ] - p4[ -_2 ];
1960  double n4z = p4[ _3 ] - p0[ 0 ];
1961  double n5x = p6[ 0 ] - p5[ -_1 ];
1962  double n5y = p5[ _2 ] - p4[ 0 ];
1963  double n5z = p5[ _3 ] - p1[ 0 ];
1964  double n6x = p6[ _1 ] - p5[ 0 ];
1965  double n6y = p6[ _2 ] - p7[ 0 ];
1966  double n6z = p6[ _3 ] - p2[ 0 ];
1967  double n7x = p7[ _1 ] - p4[ 0 ];
1968  double n7y = p6[ 0 ] - p7[ -_2 ];
1969  double n7z = p7[ _3 ] - p3[ 0 ];
1970 
1971  ndx[ 0 ] = n0x;
1972  ndx[ 1 ] = n3x - n0x;
1973  ndx[ 2 ] = n1x - n0x;
1974  ndx[ 3 ] = n2x - n3x - ndx[ 2 ];
1975  ndx[ 4 ] = n4x;
1976  ndx[ 5 ] = n7x - n4x;
1977  ndx[ 6 ] = n5x - n4x;
1978  ndx[ 7 ] = n6x - n7x - ndx[ 6 ];
1979 
1980  ndy[ 0 ] = n0y;
1981  ndy[ 1 ] = n3y - n0y;
1982  ndy[ 2 ] = n1y - n0y;
1983  ndy[ 3 ] = n2y - n3y - ndy[ 2 ];
1984  ndy[ 4 ] = n4y;
1985  ndy[ 5 ] = n7y - n4y;
1986  ndy[ 6 ] = n5y - n4y;
1987  ndy[ 7 ] = n6y - n7y - ndy[ 6 ];
1988 
1989  ndz[ 0 ] = n0z;
1990  ndz[ 1 ] = n3z - n0z;
1991  ndz[ 2 ] = n1z - n0z;
1992  ndz[ 3 ] = n2z - n3z - ndz[ 2 ];
1993  ndz[ 4 ] = n4z;
1994  ndz[ 5 ] = n7z - n4z;
1995  ndz[ 6 ] = n5z - n4z;
1996  ndz[ 7 ] = n6z - n7z - ndz[ 6 ];
1997 
1998  op = p;
1999  }
2000 
2001  double nx = ( ndx[ 0 ] + ndx[ 1 ] * xx ) + ( ndx[ 2 ] + ndx[ 3 ] * xx ) * yy;
2002  nx = ( nx + ( ( ndx[ 4 ] + ndx[ 5 ] * xx ) + ( ndx[ 6 ] + ndx[ 7 ] * xx ) * yy - nx ) * zz ) * _1_ax;
2003  double ny = ( ndy[ 0 ] + ndy[ 1 ] * xx ) + ( ndy[ 2 ] + ndy[ 3 ] * xx ) * yy;
2004  ny = ( ny + ( ( ndy[ 4 ] + ndy[ 5 ] * xx ) + ( ndy[ 6 ] + ndy[ 7 ] * xx ) * yy - ny ) * zz ) * _1_ay;
2005  double nz = ( ndz[ 0 ] + ndz[ 1 ] * xx ) + ( ndz[ 2 ] + ndz[ 3 ] * xx ) * yy;
2006  nz = ( nz + ( ( ndz[ 4 ] + ndz[ 5 ] * xx ) + ( ndz[ 6 ] + ndz[ 7 ] * xx ) * yy - nz ) * zz ) * _1_az;
2007 
2008  // 法線が反転している場合への対応
2009  double c = light.x * nx + light.y * ny + light.z * nz;
2010  c = std::sqrt( ( c * c ) / ( nx * nx + ny * ny + nz * nz + type_limits< double >::tiny( ) ) );
2011 
2012  double spec = 0.0;
2013  if( bSpecular )
2014  {
2015  spec = 2.0 * c * c - 1.0;
2016 
2017  if( spec <= 0.0 )
2018  {
2019  spec = 0;
2020  }
2021  else
2022  {
2023  spec *= spec; // 2乗
2024  spec *= spec; // 4乗
2025  spec *= spec; // 8乗
2026  spec *= spec; // 16乗
2027  spec *= spec; // 32乗
2028  spec *= spec; // 64乗
2029  //spec *= spec; // 128乗
2030  spec *= specular * 255.0;
2031  }
2032  }
2033 
2034  double lAtten = 1.0;
2035  if( bLightAtten )
2036  {
2037  double len = ( l + of ) * dlen;
2038  lAtten /= 1.0 + lightAtten * ( len * len );
2039  }
2040 
2041  double alpha = oc.alpha * sampling_step;
2042  add_intensity += alpha * add_opacity * ( oc.pixel * ( c * diffuse_ratio + ambient_ratio ) + spec ) * lAtten;
2043  add_opacity *= ( 1.0 - alpha );
2044 
2045  // 画素がレンダリング結果に与える影響がしきい値以下になった場合は終了
2046  if( add_opacity < termination )
2047  {
2048  break;
2049  }
2050 
2051  spos.x += ray_step.x;
2052  spos.y += ray_step.y;
2053  spos.z += ray_step.z;
2054  l += ray_sampling_step;
2055  }
2056  else
2057  {
2058  // この位置における物体が透明の場合は次のステップへ移行する
2059  spos += ray_step;
2060  l += ray_sampling_step;
2061 
2062  double ol = l;
2063  while( l < n )
2064  {
2065  difference_type si = volumerender::to_integer( spos.x );
2066  difference_type sj = volumerender::to_integer( spos.y );
2067  difference_type sk = volumerender::to_integer( spos.z );
2068 
2069  const_pointer p = &in( si, sj, sk );
2070 
2071  // この位置における物体が不透明の場合は次のステップへ移行する
2072  if( table[ p[ d0 ] ].has_alpha || table[ p[ d1 ] ].has_alpha ||
2073  table[ p[ d2 ] ].has_alpha || table[ p[ d3 ] ].has_alpha ||
2074  table[ p[ d4 ] ].has_alpha || table[ p[ d5 ] ].has_alpha ||
2075  table[ p[ d6 ] ].has_alpha || table[ p[ d7 ] ].has_alpha )
2076  {
2077  if( l > ol )
2078  {
2079  spos.x -= ray.x;
2080  spos.y -= ray.y;
2081  spos.z -= ray.z;
2082  l -= 1.0;
2083  }
2084  break;
2085  }
2086 
2087  double current_step = depth_map( si, sj, sk );
2088  l += current_step;
2089  spos.x += ray.x * current_step;
2090  spos.y += ray.y * current_step;
2091  spos.z += ray.z * current_step;
2092 
2093  while( l < n )
2094  {
2095  si = volumerender::to_integer( spos.x );
2096  sj = volumerender::to_integer( spos.y );
2097  sk = volumerender::to_integer( spos.z );
2098 
2099  current_step = depth_map( si, sj, sk );
2100 
2101  if( current_step <= 2.0 )
2102  {
2103  break;
2104  }
2105 
2106  l += current_step;
2107  spos.x += ray.x * current_step;
2108  spos.y += ray.y * current_step;
2109  spos.z += ray.z * current_step;
2110  }
2111  }
2112  }
2113  }
2114 
2115  out( i, j ) = static_cast< out_value_type >( mist::limits_0_255( add_intensity * ( 1.0 - add_opacity ) + background * add_opacity ) );
2116  }
2117  else
2118  {
2119  out( i, j ) = static_cast< out_value_type >( mist::limits_0_255( background ) );
2120  }
2121  }
2122  }
2123  return( true );
2124  }
2125 
2126  template < class Array1, class Array2, class DepthMap, class T >
2127  class volumerendering_thread : public mist::thread< volumerendering_thread< Array1, Array2, DepthMap, T > >
2128  {
2129  public:
2130  typedef mist::thread< volumerendering_thread< Array1, Array2, DepthMap, T > > base;
2131  typedef typename base::thread_exit_type thread_exit_type;
2132  typedef typename Array1::size_type size_type;
2133  typedef typename Array1::value_type value_type;
2134 
2135  private:
2136  size_t thread_id_;
2137  size_t thread_num_;
2138 
2139  // 入出力用の画像へのポインタ
2140  const Array1 *in_;
2141  Array2 *out_;
2142  const DepthMap *depth_map_;
2143  const volumerender::parameter *param_;
2144  const volumerender::attribute_table< T > *table_;
2145 
2146  public:
2147  void setup_parameters( const Array1 &in, Array2 &out, const DepthMap &depth_map, const volumerender::parameter &p, const volumerender::attribute_table< T > &t, size_type thread_id, size_type thread_num )
2148  {
2149  in_ = &in;
2150  out_ = &out;
2151  depth_map_ = &depth_map;
2152  param_ = &p;
2153  table_ = &t;
2154  thread_id_ = thread_id;
2155  thread_num_ = thread_num;
2156  }
2157 
2158  volumerendering_thread( size_type id = 0, size_type num = 1 ) : thread_id_( id ), thread_num_( num ),
2159  in_( NULL ), out_( NULL ), depth_map_( NULL ), param_( NULL ), table_( NULL )
2160  {
2161  }
2162  volumerendering_thread( const volumerendering_thread &p ) : base( p ), thread_id_( p.thread_id_ ), thread_num_( p.thread_num_ ),
2163  in_( p.in_ ), out_( p.out_ ), depth_map_( p.depth_map_ ),
2164  param_( p.param_ ), table_( p.table_ )
2165  {
2166  }
2167 
2168  protected:
2169  // 継承した先で必ず実装されるスレッド関数
2170  virtual thread_exit_type thread_function( )
2171  {
2172  volumerendering( *in_, *out_, *depth_map_, *param_, *table_, thread_id_, thread_num_ );
2173  return( true );
2174  }
2175  };
2176 }
2177 
2178 
2179 // ボリュームレンダリング
2180 namespace __volumerendering_controller__
2181 {
2182  // いろいろなレンダラ（色の決定方法）を組み合わせることが可能なボリュームレンダリングエンジン
2183  template < class Array1, class Array2, class DepthMap, class Renderer, class T >
2184  bool volumerendering( const Array1 &in, Array2 &out, const DepthMap &depth_map, const Renderer &renderer, const volumerender::parameter &param, const volumerender::attribute_table< T > &table, typename Array1::size_type thread_id, typename Array1::size_type thread_num )
2185  {
2186  typedef typename volumerender::parameter::vector_type vector_type;
2187  typedef typename volumerender::attribute_table< T >::attribute_type attribute_type;
2188  typedef typename volumerender::attribute_table< T >::pixel_type pixel_type;
2189  typedef typename Array1::size_type size_type;
2190  typedef typename Array1::difference_type difference_type;
2191  typedef typename Array1::value_type value_type;
2192  typedef typename Array1::const_pointer const_pointer;
2193  typedef typename Array2::value_type out_value_type;
2194 
2195  vector_type offset = param.offset;
2196  vector_type pos = param.pos - offset;
2197  const volumerender::boundingbox *box = param.box;
2198  double fovy = param.fovy;
2199  double ambient_ratio = param.ambient_ratio;
2200  double diffuse_ratio = param.diffuse_ratio;
2201  double specular = param.specular;
2202  bool bSpecular = specular > 0.0;
2203  double lightAtten = param.light_attenuation;
2204  double sampling_step = param.sampling_step;
2205  double termination = param.termination;
2206  bool bperspective = param.perspective_view;
2207 
2208  const size_type w = in.width( );
2209  const size_type h = in.height( );
2210  const size_type d = in.depth( );
2211 
2212  const size_type image_width = out.width( );
2213  const size_type image_height = out.height( );
2214 
2215  const pixel_type background = _pixel_converter_< pixel_type >::convert_to( param.background_R, param.background_G, param.background_B );
2216 
2217  // スライス座標系の実寸をワールドと考える
2218  vector_type casting_start, casting_end;
2219 
2220  const double pai = 3.1415926535897932384626433832795;
2221  double focal = ( static_cast< double >( image_height ) / 2.0 ) / std::tan( fovy * pai / 180.0 / 2.0 );
2222 
2223  double cx = static_cast< double >( image_width ) / 2.0;
2224  double cy = static_cast< double >( image_height ) / 2.0;
2225  double ax = in.reso1( );
2226  double ay = in.reso2( );
2227  double az = in.reso3( );
2228  double _1_ax = 1.0 / ax;
2229  double _1_ay = 1.0 / ay;
2230  double _1_az = 1.0 / az;
2231 
2232  double masp = ax < ay ? ax : ay;
2233  masp = masp < az ? masp : az;
2234 
2235  vector_type eZ = -param.dir.unit( ); // カメラ座標系のZ軸方向ベクトル
2236  vector_type eY = param.up.unit( ); // カメラ座標系のY軸方向ベクトル
2237  vector_type eX = ( eY * eZ ).unit( ); // カメラ座標系のX軸方向ベクトル
2238 
2239  if( param.mirror_view )
2240  {
2241  // 鏡に映ったように描画する
2242  // データの軸等に反転が必要な描画を行う際に利用する
2243  eX = -eX;
2244  }
2245 
2246  if( out.reso1( ) < out.reso2( ) )
2247  {
2248  eX *= out.reso1( ) / out.reso2( );
2249  }
2250  else
2251  {
2252  eY *= out.reso2( ) / out.reso1( );
2253  focal *= out.reso2( ) / out.reso1( );
2254  }
2255 
2256  double max_distance = pos.length( ) + std::sqrt( static_cast< double >( w * w + h * h + d * d ) );
2257 
2258  for( size_type j = thread_id ; j < image_height ; j += thread_num )
2259  {
2260  for( size_type i = 0 ; i < image_width ; i++ )
2261  {
2262  // 投影面上の点をカメラ座標系に変換
2263  vector_type Pos( static_cast< double >( i ) - cx, cy - static_cast< double >( j ), -focal );
2264 
2265  // レイ方向をカメラ座標系からワールド座標系に変換
2266  vector_type light;
2267  if( bperspective )
2268  {
2269  light = ( eX * Pos.x + eY * Pos.y + eZ * Pos.z ).unit( );
2270  }
2271  else
2272  {
2273  pos = param.pos - offset + eX * Pos.x + eY * Pos.y;
2274  light = -eZ;
2275  }
2276 
2277  pixel_type add_intensity( 0 );
2278  double add_opacity = 1;
2279 
2280  casting_start = pos;
2281  casting_end = pos + light * max_distance;
2282  vector_type normal;
2283 
2284  // 物体との衝突判定
2285  if( volumerender::check_intersection( casting_start, casting_end, box[ 0 ], normal )
2286  && volumerender::check_intersection( casting_start, casting_end, box[ 1 ], normal )
2287  && volumerender::check_intersection( casting_start, casting_end, box[ 2 ], normal )
2288  && volumerender::check_intersection( casting_start, casting_end, box[ 3 ], normal )
2289  && volumerender::check_intersection( casting_start, casting_end, box[ 4 ], normal )
2290  && volumerender::check_intersection( casting_start, casting_end, box[ 5 ], normal ) )
2291  {
2292  // 光の減衰を実現するために、カメラからの距離を測る
2293  Pos.x = ( pos.x + offset.x ) * _1_ax;
2294  Pos.y = ( pos.y + offset.y ) * _1_ay;
2295  Pos.z = ( pos.z + offset.z ) * _1_az;
2296 
2297  // ワールド座標系からスライス座標系に変換する
2298  // 以降は、全てスライス座標系で計算する
2299  casting_start.x = ( casting_start.x + offset.x ) * _1_ax;
2300  casting_start.y = ( casting_start.y + offset.y ) * _1_ay;
2301  casting_start.z = ( casting_start.z + offset.z ) * _1_az;
2302  casting_end.x = ( casting_end.x + offset.x ) * _1_ax;
2303  casting_end.y = ( casting_end.y + offset.y ) * _1_ay;
2304  casting_end.z = ( casting_end.z + offset.z ) * _1_az;
2305 
2306  vector_type spos = casting_start;
2307  vector_type ray = ( casting_end - casting_start ).unit( );
2308 
2309  // 光の減衰の距離を実測に直すためのパラメータ
2310  double dlen = vector_type( ray.x * ax, ray.y * ay, ray.z * az ).length( );
2311 
2312  // 直方体画素の画像上では方向によってサンプリング間隔が変わってしまう問題に対応
2313  double ray_sampling_step = sampling_step * masp / dlen;
2314 
2315  vector_type ray_step = ray * ray_sampling_step;
2316 
2317  double n = ( casting_end - casting_start ).length( );
2318  double l = 0, of = ( Pos - casting_start ).length( );
2319 
2320  while( l < n )
2321  {
2322  difference_type si = volumerender::to_integer( spos.x );
2323  difference_type sj = volumerender::to_integer( spos.y );
2324  difference_type sk = volumerender::to_integer( spos.z );
2325 
2326  // この位置における物体が不透明の場合は次のステップへ移行する
2327  if( renderer.check( si, sj, sk ) )
2328  {
2329  if( l > 0 )
2330  {
2331  spos.x -= ray.x;
2332  spos.y -= ray.y;
2333  spos.z -= ray.z;
2334  l -= 1.0;
2335  }
2336  break;
2337  }
2338 
2339  double current_step = depth_map( si, sj, sk );
2340  l += current_step;
2341  spos.x += ray.x * current_step;
2342  spos.y += ray.y * current_step;
2343  spos.z += ray.z * current_step;
2344  }
2345 
2346  // 端まで到達した場合は何もしない
2347  if( l >= n )
2348  {
2349  out( i, j ) = static_cast< out_value_type >( mist::limits_0_255( background ) );
2350  continue;
2351  }
2352 
2353  while( l < n )
2354  {
2355  difference_type si = volumerender::to_integer( spos.x );
2356  difference_type sj = volumerender::to_integer( spos.y );
2357  difference_type sk = volumerender::to_integer( spos.z );
2358 
2359  double xx = spos.x - si;
2360  double yy = spos.y - sj;
2361  double zz = spos.z - sk;
2362 
2363  attribute_type oc;
2364 
2365  if( renderer.render( si, sj, sk, xx, yy, zz, oc ) )
2366  {
2367  double lAtten = 1.0;
2368  if( lightAtten > 0.0 )
2369  {
2370  double len = ( l + of ) * dlen;
2371  lAtten /= 1.0 + lightAtten * ( len * len );
2372  }
2373 
2374  double c = light.inner( renderer.normal( si, sj, sk, xx, yy, zz ) );
2375  c = c < 0.0 ? -c : c;
2376 
2377  double spec = 0.0;
2378  if( bSpecular )
2379  {
2380  spec = 2.0 * c * c - 1.0;
2381 
2382  if( spec <= 0.0 )
2383  {
2384  spec = 0;
2385  }
2386  else
2387  {
2388  spec *= spec; // 2乗
2389  spec *= spec; // 4乗
2390  spec *= spec; // 8乗
2391  spec *= spec; // 16乗
2392  spec *= spec; // 32乗
2393  spec *= spec; // 64乗
2394  //spec *= spec; // 128乗
2395  spec *= specular * 255.0;
2396  }
2397  }
2398 
2399  c = c * diffuse_ratio + ambient_ratio;
2400 
2401  double alpha = oc.alpha * sampling_step;
2402  add_intensity += alpha * add_opacity * ( oc.pixel * c + spec ) * lAtten;
2403  add_opacity *= ( 1.0 - alpha );
2404 
2405  // 画素がレンダリング結果に与える影響がしきい値以下になった場合は終了
2406  if( add_opacity < termination )
2407  {
2408  break;
2409  }
2410 
2411  spos.x += ray_step.x;
2412  spos.y += ray_step.y;
2413  spos.z += ray_step.z;
2414  l += ray_sampling_step;
2415  }
2416  else
2417  {
2418  // この位置における物体が透明の場合は次のステップへ移行する
2419  spos += ray_step;
2420  l += ray_sampling_step;
2421 
2422  double ol = l;
2423  while( l < n )
2424  {
2425  difference_type si = volumerender::to_integer( spos.x );
2426  difference_type sj = volumerender::to_integer( spos.y );
2427  difference_type sk = volumerender::to_integer( spos.z );
2428 
2429  // この位置における物体が不透明の場合は次のステップへ移行する
2430  if( renderer.check( si, sj, sk ) )
2431  {
2432  if( l > ol )
2433  {
2434  spos.x -= ray.x;
2435  spos.y -= ray.y;
2436  spos.z -= ray.z;
2437  l -= 1.0;
2438  }
2439  break;
2440  }
2441 
2442  double current_step = depth_map( si, sj, sk );
2443  l += current_step;
2444  spos.x += ray.x * current_step;
2445  spos.y += ray.y * current_step;
2446  spos.z += ray.z * current_step;
2447  }
2448  }
2449  }
2450 
2451  out( i, j ) = static_cast< out_value_type >( mist::limits_0_255( add_intensity * ( 1.0 - add_opacity ) + background * add_opacity ) );
2452  }
2453  else
2454  {
2455  out( i, j ) = static_cast< out_value_type >( mist::limits_0_255( background ) );
2456  }
2457  }
2458  }
2459 
2460  return( true );
2461  }
2462 
2463  // いろいろなレンダラ（色の決定方法）を組み合わせ，指定した位置のレイが最初にあたる位置を返す関数
2464  template < class Array, class DepthMap, class Renderer, class T >
2465  typename volumerender::parameter::vector_type collision_detection( const Array &in, typename Array::size_type image_width, typename Array::size_type image_height, double resoX, double resoY,
2466  const DepthMap &depth_map, const Renderer &renderer, const volumerender::parameter &param, const volumerender::attribute_table< T > &table, typename Array::size_type i, typename Array::size_type j )
2467  {
2468  typedef typename volumerender::parameter::vector_type vector_type;
2469  typedef typename volumerender::attribute_table< T >::attribute_type attribute_type;
2470  typedef typename volumerender::attribute_table< T >::pixel_type pixel_type;
2471  typedef typename Array::size_type size_type;
2472  typedef typename Array::difference_type difference_type;
2473  typedef typename Array::value_type value_type;
2474  typedef typename Array::const_pointer const_pointer;
2475 
2476  vector_type offset = param.offset;
2477  vector_type pos = param.pos - offset;
2478  double fovy = param.fovy;
2479  double ambient_ratio = param.ambient_ratio;
2480  double diffuse_ratio = param.diffuse_ratio;
2481  double specular = param.specular;
2482  bool bSpecular = specular > 0.0;
2483  const volumerender::boundingbox *box = param.box;
2484  double lightAtten = param.light_attenuation;
2485  double sampling_step = param.sampling_step;
2486  double termination = param.termination;
2487  bool bperspective = param.perspective_view;
2488 
2489  const size_type w = in.width( );
2490  const size_type h = in.height( );
2491  const size_type d = in.depth( );
2492 
2493  // スライス座標系の実寸をワールドと考える
2494  vector_type casting_start, casting_end;
2495 
2496  const double pai = 3.1415926535897932384626433832795;
2497  double focal = ( static_cast< double >( image_height ) / 2.0 ) / std::tan( fovy * pai / 180.0 / 2.0 );
2498 
2499  double cx = static_cast< double >( image_width ) / 2.0;
2500  double cy = static_cast< double >( image_height ) / 2.0;
2501  double ax = in.reso1( );
2502  double ay = in.reso2( );
2503  double az = in.reso3( );
2504  double _1_ax = 1.0 / ax;
2505  double _1_ay = 1.0 / ay;
2506  double _1_az = 1.0 / az;
2507 
2508  double asp = resoY / resoX;
2509 
2510  double masp = ax < ay ? ax : ay;
2511  masp = masp < az ? masp : az;
2512 
2513  vector_type eZ = -param.dir.unit( ); // カメラ座標系のZ軸方向ベクトル
2514  vector_type eY = param.up.unit( ); // カメラ座標系のY軸方向ベクトル
2515  vector_type eX = ( eY * eZ ).unit( ); // カメラ座標系のX軸方向ベクトル
2516 
2517  if( param.mirror_view )
2518  {
2519  // 鏡に映ったように描画する
2520  // データの軸等に反転が必要な描画を行う際に利用する
2521  eX = -eX;
2522  }
2523 
2524  if( resoX < resoY )
2525  {
2526  eX *= resoX / resoY;
2527  }
2528  else
2529  {
2530  eY *= resoY / resoX;
2531  focal *= resoY / resoX;
2532  }
2533 
2534  double max_distance = pos.length( ) + std::sqrt( static_cast< double >( w * w + h * h + d * d ) );
2535 
2536  {
2537  // 投影面上の点をカメラ座標系に変換
2538  vector_type Pos( static_cast< double >( i ) - cx, cy - static_cast< double >( j ), -focal );
2539 
2540  // レイ方向をカメラ座標系からワールド座標系に変換
2541  vector_type light;
2542  if( bperspective )
2543  {
2544  light = ( eX * Pos.x + eY * Pos.y + eZ * Pos.z ).unit( );
2545  }
2546  else
2547  {
2548  pos = param.pos - offset + eX * Pos.x + eY * Pos.y;
2549  light = -eZ;
2550  }
2551 
2552  double add_opacity = 1;
2553  double old_add_opacity = 0.0;
2554 
2555  casting_start = pos;
2556  casting_end = pos + light * max_distance;
2557 
2558  // 物体との衝突判定
2559  if( volumerender::check_intersection( casting_start, casting_end, box[ 0 ] )
2560  && volumerender::check_intersection( casting_start, casting_end, box[ 1 ] )
2561  && volumerender::check_intersection( casting_start, casting_end, box[ 2 ] )
2562  && volumerender::check_intersection( casting_start, casting_end, box[ 3 ] )
2563  && volumerender::check_intersection( casting_start, casting_end, box[ 4 ] )
2564  && volumerender::check_intersection( casting_start, casting_end, box[ 5 ] ) )
2565  {
2566  // ワールド座標系からスライス座標系に変換する
2567  // 以降は、全てスライス座標系で計算する
2568  casting_start.x = ( casting_start.x + offset.x ) * _1_ax;
2569  casting_start.y = ( casting_start.y + offset.y ) * _1_ay;
2570  casting_start.z = ( casting_start.z + offset.z ) * _1_az;
2571  casting_end.x = ( casting_end.x + offset.x ) * _1_ax;
2572  casting_end.y = ( casting_end.y + offset.y ) * _1_ay;
2573  casting_end.z = ( casting_end.z + offset.z ) * _1_az;
2574 
2575  vector_type opos = casting_start;
2576  vector_type spos = casting_start;
2577  vector_type ray = ( casting_end - casting_start ).unit( );
2578 
2579  // 光の減衰の距離を実測に直すためのパラメータ
2580  double dlen = vector_type( ray.x * ax, ray.y * ay, ray.z * az ).length( );
2581 
2582  // 直方体画素の画像上では方向によってサンプリング間隔が変わってしまう問題に対応
2583  double ray_sampling_step = sampling_step * masp / dlen;
2584 
2585  vector_type ray_step = ray * ray_sampling_step;
2586 
2587  double n = ( casting_end - casting_start ).length( );
2588  double l = 0, of = ( Pos - casting_start ).length( );
2589 
2590  while( l < n )
2591  {
2592  difference_type si = volumerender::to_integer( spos.x );
2593  difference_type sj = volumerender::to_integer( spos.y );
2594  difference_type sk = volumerender::to_integer( spos.z );
2595 
2596  // この位置における物体が不透明の場合は次のステップへ移行する
2597  if( renderer.check( si, sj, sk ) )
2598  {
2599  if( l > 0 )
2600  {
2601  spos.x -= ray.x;
2602  spos.y -= ray.y;
2603  spos.z -= ray.z;
2604  l -= 1.0;
2605  }
2606  break;
2607  }
2608 
2609  double current_step = depth_map( si, sj, sk );
2610  l += current_step;
2611  spos.x += ray.x * current_step;
2612  spos.y += ray.y * current_step;
2613  spos.z += ray.z * current_step;
2614  }
2615 
2616  while( l < n )
2617  {
2618  difference_type si = volumerender::to_integer( spos.x );
2619  difference_type sj = volumerender::to_integer( spos.y );
2620  difference_type sk = volumerender::to_integer( spos.z );
2621 
2622  double xx = spos.x - si;
2623  double yy = spos.y - sj;
2624  double zz = spos.z - sk;
2625 
2626  attribute_type oc;
2627 
2628  if( renderer.render( si, sj, sk, xx, yy, zz, oc ) )
2629  {
2630  double alpha = oc.alpha * sampling_step;
2631  double aopacity = alpha * add_opacity;
2632 
2633  // 画素がレンダリング結果に与える影響がしきい値以下になった場合は終了
2634  if( old_add_opacity < aopacity )
2635  {
2636  opos = spos;
2637  old_add_opacity = aopacity;
2638  }
2639 
2640  add_opacity = add_opacity * ( 1.0 - alpha );
2641 
2642  // 画素がレンダリング結果に与える影響がしきい値以下になった場合は終了
2643  if( add_opacity < termination )
2644  {
2645  break;
2646  }
2647 
2648  spos.x += ray_step.x;
2649  spos.y += ray_step.y;
2650  spos.z += ray_step.z;
2651  l += ray_sampling_step;
2652  }
2653  else
2654  {
2655  // この位置における物体が透明の場合は次のステップへ移行する
2656  spos += ray_step;
2657  l += ray_sampling_step;
2658 
2659  double ol = l;
2660  while( l < n )
2661  {
2662  difference_type si = volumerender::to_integer( spos.x );
2663  difference_type sj = volumerender::to_integer( spos.y );
2664  difference_type sk = volumerender::to_integer( spos.z );
2665 
2666  // この位置における物体が不透明の場合は次のステップへ移行する
2667  if( renderer.check( si, sj, sk ) )
2668  {
2669  if( l > ol )
2670  {
2671  spos.x -= ray.x;
2672  spos.y -= ray.y;
2673  spos.z -= ray.z;
2674  l -= 1.0;
2675  }
2676  break;
2677  }
2678 
2679  double current_step = depth_map( si, sj, sk );
2680  l += current_step;
2681  spos.x += ray.x * current_step;
2682  spos.y += ray.y * current_step;
2683  spos.z += ray.z * current_step;
2684  }
2685  }
2686  }
2687 
2688  opos.x *= ax;
2689  opos.y *= ay;
2690  opos.z *= az;
2691  return( opos );
2692  }
2693  else
2694  {
2695  double infinity = type_limits< double >::maximum( );
2696  return( vector_type( infinity, infinity, infinity ) );
2697  }
2698  }
2699  }
2700 
2701 
2702 
2703  template < class Array1, class Array2, class DepthMap, class Renderer, class T >
2704  class volumerendering_thread : public mist::thread< volumerendering_thread< Array1, Array2, DepthMap, Renderer, T > >
2705  {
2706  public:
2707  typedef mist::thread< volumerendering_thread< Array1, Array2, DepthMap, Renderer, T > > base;
2708  typedef typename base::thread_exit_type thread_exit_type;
2709  typedef typename Array1::size_type size_type;
2710  typedef typename Array1::value_type value_type;
2711 
2712  private:
2713  size_t thread_id_;
2714  size_t thread_num_;
2715 
2716  // 入出力用の画像へのポインタ
2717  const Array1 *in_;
2718  Array2 *out_;
2719  const DepthMap *depth_map_;
2720  const Renderer *renderer_;
2721  const volumerender::parameter *param_;
2722  const volumerender::attribute_table< T > *table_;
2723 
2724  public:
2725  void setup_parameters( const Array1 &in, Array2 &out, const DepthMap &depth_map, const Renderer &renderer, const volumerender::parameter &p, const volumerender::attribute_table< T > &t, size_type thread_id, size_type thread_num )
2726  {
2727  in_ = &in;
2728  out_ = &out;
2729  depth_map_ = &depth_map;
2730  renderer_ = &renderer;
2731  param_ = &p;
2732  table_ = &t;
2733  thread_id_ = thread_id;
2734  thread_num_ = thread_num;
2735  }
2736 
2737  volumerendering_thread( size_type id = 0, size_type num = 1 ) : thread_id_( id ), thread_num_( num ),
2738  in_( NULL ), out_( NULL ), depth_map_( NULL ), renderer_( NULL ), param_( NULL ), table_( NULL )
2739  {
2740  }
2741 
2742  protected:
2743  // 継承した先で必ず実装されるスレッド関数
2744  virtual thread_exit_type thread_function( )
2745  {
2746  volumerendering( *in_, *out_, *depth_map_, *renderer_, *param_, *table_, thread_id_, thread_num_ );
2747  return( true );
2748  }
2749  };
2750 }
2751 
2752 
2753 namespace __mip_controller__
2754 {
2755  template < class Array1, class Array2 >
2756  bool mip( const Array1 &in, Array2 &out, const volumerender::parameter &param, typename Array1::size_type thread_id, typename Array1::size_type thread_num )
2757  {
2758  typedef typename volumerender::parameter::vector_type vector_type;
2759  typedef typename Array1::size_type size_type;
2760  typedef typename Array1::difference_type difference_type;
2761  typedef typename Array1::value_type value_type;
2762  typedef typename Array1::const_pointer const_pointer;
2763  typedef typename Array2::value_type out_value_type;
2764 
2765  vector_type offset = param.offset;
2766  vector_type pos = param.pos - offset;
2767  const volumerender::boundingbox *box = param.box;
2768  const double sampling_step = param.sampling_step;
2769 
2770  const size_type w = in.width( );
2771  const size_type h = in.height( );
2772  const size_type d = in.depth( );
2773 
2774  const size_type image_width = out.width( );
2775  const size_type image_height = out.height( );
2776 
2777  // 高速にアドレス計算を行うためのポインタの差分
2778  difference_type d0, d1, d2, d3, d4, d5, d6, d7, _1, _2, _3;
2779  {
2780  difference_type cx = in.width( ) / 2;
2781  difference_type cy = in.height( ) / 2;
2782  difference_type cz = in.depth( ) / 2;
2783  const_pointer ppp = &in( cx, cy, cz );
2784  d0 = 0;
2785  d1 = &in( cx , cy + 1, cz ) - ppp;
2786  d2 = &in( cx + 1, cy + 1, cz ) - ppp;
2787  d3 = &in( cx + 1, cy , cz ) - ppp;
2788  d4 = &in( cx , cy , cz + 1 ) - ppp;
2789  d5 = &in( cx , cy + 1, cz + 1 ) - ppp;
2790  d6 = &in( cx + 1, cy + 1, cz + 1 ) - ppp;
2791  d7 = &in( cx + 1, cy , cz + 1 ) - ppp;
2792  _1 = &in( cx + 1, cy , cz ) - ppp;
2793  _2 = &in( cx , cy + 1, cz ) - ppp;
2794  _3 = &in( cx , cy , cz + 1 ) - ppp;
2795  }
2796 
2797  // スライス座標系の実寸をワールドと考える
2798  vector_type casting_start, casting_end;
2799 
2800  double iasp = 1.0;
2801  double cx = static_cast< double >( image_width ) / 2.0;
2802  double cy = static_cast< double >( image_height ) / 2.0;
2803  double ax = in.reso1( );
2804  double ay = in.reso2( );
2805  double az = in.reso3( );
2806  double _1_ax = 1.0 / ax;
2807  double _1_ay = 1.0 / ay;
2808  double _1_az = 1.0 / az;
2809 
2810  double masp = ax < ay ? ax : ay;
2811  masp = masp < az ? masp : az;
2812 
2813  vector_type eZ = -param.dir.unit( ); // カメラ座標系のZ軸方向ベクトル
2814  vector_type eY = param.up.unit( ); // カメラ座標系のY軸方向ベクトル
2815  vector_type eX = ( eY * eZ ).unit( ); // カメラ座標系のX軸方向ベクトル
2816 
2817  if( param.mirror_view )
2818  {
2819  // 鏡に映ったように描画する
2820  // データの軸等に反転が必要な描画を行う際に利用する
2821  eX = -eX;
2822  }
2823 
2824  if( out.reso1( ) < out.reso2( ) )
2825  {
2826  eX *= out.reso1( ) / out.reso2( );
2827  }
2828  else
2829  {
2830  eY *= out.reso2( ) / out.reso1( );
2831  }
2832 
2833  if( out.reso1( ) * image_width >= out.reso2( ) * image_height )
2834  {
2835  iasp = out.reso1( );
2836  }
2837  else
2838  {
2839  iasp = out.reso2( );
2840  }
2841 
2842  double max_distance = pos.length( ) + std::sqrt( static_cast< double >( w * w + h * h + d * d ) );
2843 
2844  for( size_type j = thread_id ; j < image_height ; j += thread_num )
2845  {
2846  for( size_type i = 0 ; i < image_width ; i++ )
2847  {
2848  casting_start = pos + eZ * max_distance + ( static_cast< double >( i ) - cx ) * iasp * eX + ( cy - static_cast< double >( j ) ) * iasp * eY;
2849  casting_end = casting_start - eZ * max_distance * 2.0;
2850 
2851  // 物体との衝突判定
2852  if( volumerender::check_intersection( casting_start, casting_end, box[ 0 ] )
2853  && volumerender::check_intersection( casting_start, casting_end, box[ 1 ] )
2854  && volumerender::check_intersection( casting_start, casting_end, box[ 2 ] )
2855  && volumerender::check_intersection( casting_start, casting_end, box[ 3 ] )
2856  && volumerender::check_intersection( casting_start, casting_end, box[ 4 ] )
2857  && volumerender::check_intersection( casting_start, casting_end, box[ 5 ] ) )
2858  {
2859  // ワールド座標系からスライス座標系に変換する
2860  // 以降は、全てスライス座標系で計算する
2861  casting_start.x = ( casting_start.x + offset.x ) * _1_ax;
2862  casting_start.y = ( casting_start.y + offset.y ) * _1_ay;
2863  casting_start.z = ( casting_start.z + offset.z ) * _1_az;
2864  casting_end.x = ( casting_end.x + offset.x ) * _1_ax;
2865  casting_end.y = ( casting_end.y + offset.y ) * _1_ay;
2866  casting_end.z = ( casting_end.z + offset.z ) * _1_az;
2867 
2868  vector_type spos = casting_start;
2869  vector_type ray = ( casting_end - casting_start ).unit( );
2870 
2871  vector_type ray_step = ray * sampling_step;
2872 
2873  double n = ( casting_end - casting_start ).length( );
2874  double l = 0, maxct = mist::type_limits< double >::minimum( );
2875 
2876  while( l < n )
2877  {
2878  difference_type si = volumerender::to_integer( spos.x );
2879  difference_type sj = volumerender::to_integer( spos.y );
2880  difference_type sk = volumerender::to_integer( spos.z );
2881 
2882  double xx = spos.x - si;
2883  double yy = spos.y - sj;
2884  double zz = spos.z - sk;
2885 
2886  const_pointer p = &in( si, sj, sk );
2887 
2888  double nct0 = p[ d0 ];
2889  double nct1 = p[ d3 ] - p[ d0 ];
2890  double nct2 = p[ d1 ] - p[ d0 ];
2891  double nct3 = p[ d2 ] - p[ d3 ] - nct2;
2892  double nct4 = p[ d4 ];
2893  double nct5 = p[ d7 ] - p[ d4 ];
2894  double nct6 = p[ d5 ] - p[ d4 ];
2895  double nct7 = p[ d6 ] - p[ d7 ] - nct6;
2896 
2897  // CT値に線形補間を用いて取得
2898  double ct = ( nct0 + nct1 * xx ) + ( nct2 + nct3 * xx ) * yy;
2899  ct += ( ( nct4 + nct5 * xx ) + ( nct6 + nct7 * xx ) * yy - ct ) * zz;
2900 
2901  if( maxct < ct )
2902  {
2903  maxct = ct;
2904  }
2905 
2906  spos.x += ray_step.x;
2907  spos.y += ray_step.y;
2908  spos.z += ray_step.z;
2909  l += sampling_step;
2910  }
2911 
2912  out( i, j ) = static_cast< out_value_type >( maxct );
2913  }
2914  else
2915  {
2916  out( i, j ) = 0;
2917  }
2918  }
2919  }
2920 
2921  return( true );
2922  }
2923 
2924 
2925  template < class Array1, class Array2 >
2926  class mip_thread : public mist::thread< mip_thread< Array1, Array2 > >
2927  {
2928  public:
2929  typedef mist::thread< mip_thread< Array1, Array2 > > base;
2930  typedef typename base::thread_exit_type thread_exit_type;
2931  typedef typename Array1::size_type size_type;
2932  typedef typename Array1::value_type value_type;
2933 
2934  private:
2935  size_t thread_id_;
2936  size_t thread_num_;
2937 
2938  // 入出力用の画像へのポインタ
2939  const Array1 *in_;
2940  Array2 *out_;
2941  const volumerender::parameter *param_;
2942 
2943  public:
2944  void setup_parameters( const Array1 &in, Array2 &out, const volumerender::parameter &p, size_type thread_id, size_type thread_num )
2945  {
2946  in_ = &in;
2947  out_ = &out;
2948  param_ = &p;
2949  thread_id_ = thread_id;
2950  thread_num_ = thread_num;
2951  }
2952 
2953  mip_thread( size_type id = 0, size_type num = 1 ) : thread_id_( id ), thread_num_( num ), in_( NULL ), out_( NULL ), param_( NULL )
2954  {
2955  }
2956 
2957  protected:
2958  // 継承した先で必ず実装されるスレッド関数
2959  virtual thread_exit_type thread_function( )
2960  {
2961  mip( *in_, *out_, *param_, thread_id_, thread_num_ );
2962  return( true );
2963  }
2964  };
2965 }
2966 
2967 
2968 
2972 
2973 
2977 
2978 
2995 template < class Array1, class Array2, class DepthMap, class Renderer, class ATTRIBUTETYPE >
2996 bool volumerendering( const Array1 &in, Array2 &out, const DepthMap &dmap, const Renderer &renderer, const volumerender::parameter &param, const volumerender::attribute_table< ATTRIBUTETYPE > &table, typename Array1::size_type thread_num = 0 )
2997 {
2998  if( is_same_object( in, out ) || in.empty( ) )
2999  {
3000  return( false );
3001  }
3002 
3003  typedef typename Array1::size_type size_type;
3004  typedef __volumerendering_controller__::volumerendering_thread< Array1, Array2, DepthMap, Renderer, ATTRIBUTETYPE > volumerendering_thread;
3005 
3006  if( thread_num == 0 )
3007  {
3008  thread_num = static_cast< size_type >( get_cpu_num( ) );
3009  }
3010 
3011  volumerendering_thread *thread = new volumerendering_thread[ thread_num ];
3012 
3013  size_type i;
3014  for( i = 0 ; i < thread_num ; i++ )
3015  {
3016  thread[ i ].setup_parameters( in, out, dmap, renderer, param, table, i, thread_num );
3017  }
3018 
3019  // スレッドを実行して，終了まで待機する
3020  do_threads_( thread, thread_num );
3021 
3022  delete [] thread;
3023 
3024  return( true );
3025 }
3026 
3027 
3043 template < class Array1, class Array2, class DepthMap, class ATTRIBUTETYPE >
3044 bool volumerendering( const Array1 &in, Array2 &out, const DepthMap &dmap, const volumerender::parameter &param,
3045  const volumerender::attribute_table< ATTRIBUTETYPE > &table, typename Array1::size_type thread_num = 0 )
3046 {
3047  if( param.value_interpolation )
3048  {
3049  typedef rendering_helper::value_interpolation< Array1, ATTRIBUTETYPE > Renderer;
3050  return( volumerendering( in, out, dmap, Renderer( in, param, table ), param, table, thread_num ) );
3051  }
3052  else
3053  {
3054  typedef rendering_helper::color_interpolation< Array1, ATTRIBUTETYPE > Renderer;
3055  return( volumerendering( in, out, dmap, Renderer( in, param, table ), param, table, thread_num ) );
3056  }
3057 
3058  return( true );
3059 }
3060 
3061 
3076 template < class Array1, class Array2, class ATTRIBUTETYPE >
3077 bool volumerendering( const Array1 &in, Array2 &out, const volumerender::parameter &param, const volumerender::attribute_table< ATTRIBUTETYPE > &table, typename Array1::size_type thread_num = 0 )
3078 {
3079  return( volumerendering( in, out, volumerender::no_depth_map( ), param, table, thread_num ) );
3080 }
3081 
3082 
3083 
3101 template < class Array1, class Array2, class Array3, class DepthMap, class ATTRIBUTETYPE >
3102 bool volumerendering( const Array1 &in, const Array2 &mk, Array3 &out, const DepthMap &dmap, const volumerender::parameter &param, const volumerender::attribute_table< ATTRIBUTETYPE > &table, const volumerender::attribute_table< ATTRIBUTETYPE > &mktable, typename Array1::size_type thread_num = 0 )
3103 {
3104  if( in.width( ) != mk.width( ) || in.height( ) != mk.height( ) || in.depth( ) != mk.depth( ) )
3105  {
3106  return( false );
3107  }
3108 
3109  typedef rendering_helper::value_interpolation_with_mark< Array1, Array2, ATTRIBUTETYPE > Renderer;
3110  return( volumerendering( in, out, dmap, Renderer( in, mk, param, table, mktable ), param, table, thread_num ) );
3111 }
3112 
3113 
3114 
3133 template < class Array1, class Array2, class Array3, class DepthMap, class ATTRIBUTETYPE >
3134 bool volumerendering( const Array1 &in, const Array2 &mk, Array3 &out, const DepthMap &dmap, const volumerender::parameter &param, const volumerender::attribute_table< ATTRIBUTETYPE > &table, const volumerender::attribute_table< ATTRIBUTETYPE > &mktable, bool apply_and_operation, typename Array1::size_type thread_num = 0 )
3135 {
3136  if( in.width( ) != mk.width( ) || in.height( ) != mk.height( ) || in.depth( ) != mk.depth( ) )
3137  {
3138  return( false );
3139  }
3140 
3141  if( apply_and_operation )
3142  {
3143  typedef rendering_helper::value_interpolation_and_mark< Array1, Array2, ATTRIBUTETYPE > Renderer;
3144  return( volumerendering( in, out, dmap, Renderer( in, mk, param, table, mktable ), param, table, thread_num ) );
3145  }
3146  else
3147  {
3148  typedef rendering_helper::value_interpolation_or_mark< Array1, Array2, ATTRIBUTETYPE > Renderer;
3149  return( volumerendering( in, out, dmap, Renderer( in, mk, param, table, mktable ), param, table, thread_num ) );
3150  }
3151 }
3152 
3153 
3170 template < class Array1, class Array2, class Array3, class ATTRIBUTETYPE >
3171 bool volumerendering( const Array1 &in, const Array2 &mk, Array3 &out, const volumerender::parameter &param, const volumerender::attribute_table< ATTRIBUTETYPE > &table, const volumerender::attribute_table< ATTRIBUTETYPE > &mktable, typename Array1::size_type thread_num = 0 )
3172 {
3173  return( volumerendering( in, mk, out, volumerender::no_depth_map( ), param, table, mktable, thread_num ) );
3174 }
3175 
3176 
3194 template < class Array1, class Array2, class Array3, class ATTRIBUTETYPE >
3195 bool volumerendering( const Array1 &in, const Array2 &mk, Array3 &out, const volumerender::parameter &param, const volumerender::attribute_table< ATTRIBUTETYPE > &table, const volumerender::attribute_table< ATTRIBUTETYPE > &mktable, bool apply_and_operation, typename Array1::size_type thread_num = 0 )
3196 {
3197  return( volumerendering( in, mk, out, volumerender::no_depth_map( ), param, table, mktable, apply_and_operation, thread_num ) );
3198 }
3199 
3200 
3201 namespace specialized
3202 {
3218  template < class Array1, class Array2, class DepthMap, class ATTRIBUTETYPE >
3219  bool volumerendering( const Array1 &in, Array2 &out, const DepthMap &dmap, const volumerender::parameter &param, const volumerender::attribute_table< ATTRIBUTETYPE > &table, typename Array1::size_type thread_num = 0 )
3220  {
3221  if( is_same_object( in, out ) || in.empty( ) )
3222  {
3223  return( false );
3224  }
3225 
3226  typedef typename Array1::size_type size_type;
3227  typedef __volumerendering_specialized__::volumerendering_thread< Array1, Array2, DepthMap, ATTRIBUTETYPE > volumerendering_thread;
3228 
3229  if( thread_num == 0 )
3230  {
3231  thread_num = static_cast< size_type >( get_cpu_num( ) );
3232  }
3233 
3234  volumerendering_thread *thread = new volumerendering_thread[ thread_num ];
3235 
3236  size_type i;
3237  for( i = 0 ; i < thread_num ; i++ )
3238  {
3239  thread[ i ].setup_parameters( in, out, dmap, param, table, i, thread_num );
3240  }
3241 
3242  // スレッドを実行して，終了まで待機する
3243  do_threads_( thread, thread_num );
3244 
3245  delete [] thread;
3246 
3247  return( true );
3248  }
3249 
3250 
3265  template < class Array1, class Array2, class ATTRIBUTETYPE >
3266  bool volumerendering( const Array1 &in, Array2 &out, const volumerender::parameter &param, const volumerender::attribute_table< ATTRIBUTETYPE > &table, typename Array1::size_type thread_num = 0 )
3267  {
3268  return( specialized::volumerendering( in, out, volumerender::no_depth_map( ), param, table, thread_num ) );
3269  }
3270 }
3271 
3272 
3286 template < class Array1, class Array2 >
3287 bool mip( const Array1 &in, Array2 &out, const volumerender::parameter &p, typename Array1::size_type thread_num = 0 )
3288 {
3289  if( is_same_object( in, out ) || in.empty( ) )
3290  {
3291  return( false );
3292  }
3293 
3294  typedef typename Array1::size_type size_type;
3295  typedef __mip_controller__::mip_thread< Array1, Array2 > mip_thread;
3296 
3297  if( thread_num == 0 )
3298  {
3299  thread_num = static_cast< size_type >( get_cpu_num( ) );
3300  }
3301 
3302  mip_thread *thread = new mip_thread[ thread_num ];
3303 
3304  size_type i;
3305  for( i = 0 ; i < thread_num ; i++ )
3306  {
3307  thread[ i ].setup_parameters( in, out, p, i, thread_num );
3308  }
3309 
3310  // スレッドを実行して，終了まで待機する
3311  do_threads_( thread, thread_num );
3312 
3313  delete [] thread;
3314 
3315  return( true );
3316 }
3317 
3337 template < class Array, class DepthMap, class Renderer, class T >
3338 volumerender::parameter::vector_type collision_detection( const Array &in, typename Array::size_type image_width, typename Array::size_type image_height, double resoX, double resoY,
3339  const DepthMap &dmap, const Renderer &renderer, const volumerender::parameter &param, const volumerender::attribute_table< T > &table, typename Array::size_type i, typename Array::size_type j )
3340 {
3341  return( __volumerendering_controller__::collision_detection( in, image_width, image_height, resoX, resoY, dmap, renderer, param, table, i, j ) );
3342 }
3343 
3360 template < class Array1, class Array2, class DepthMap, class Renderer, class T >
3361 volumerender::parameter::vector_type collision_detection( const Array1 &in, const Array2 &out, const DepthMap &dmap, const Renderer &renderer, const volumerender::parameter &param, const volumerender::attribute_table< T > &table, typename Array1::size_type i, typename Array1::size_type j )
3362 {
3363  return( __volumerendering_controller__::collision_detection( in, out.width( ), out.height( ), out.reso1( ), out.reso2( ), dmap, renderer, param, table, i, j ) );
3364 }
3365 
3366 
3385 template < class Array, class DepthMap, class T >
3386 volumerender::parameter::vector_type collision_detection( const Array &in, typename Array::size_type image_width, typename Array::size_type image_height, double resoX, double resoY, const DepthMap &dmap, const volumerender::parameter &param, const volumerender::attribute_table< T > &table, typename Array::size_type i, typename Array::size_type j )
3387 {
3388  if( param.value_interpolation )
3389  {
3390  typedef rendering_helper::value_interpolation< Array, T > Renderer;
3391  return( collision_detection( in, image_width, image_height, resoX, resoY, dmap, Renderer( in, param, table ), param, table, i, j ) );
3392  }
3393  else
3394  {
3395  typedef rendering_helper::color_interpolation< Array, T > Renderer;
3396  return( collision_detection( in, image_width, image_height, resoX, resoY, dmap, Renderer( in, param, table ), param, table, i, j ) );
3397  }
3398 }
3399 
3400 
3416 template < class Array1, class Array2, class DepthMap, class T >
3417 volumerender::parameter::vector_type collision_detection( const Array1 &in, const Array2 &out, const DepthMap &dmap, const volumerender::parameter &param, const volumerender::attribute_table< T > &table, typename Array1::size_type i, typename Array1::size_type j )
3418 {
3419  if( param.value_interpolation )
3420  {
3421  typedef rendering_helper::value_interpolation< Array1, T > Renderer;
3422  return( collision_detection( in, out, dmap, Renderer( in, param, table ), param, table, i, j ) );
3423  }
3424  else
3425  {
3426  typedef rendering_helper::color_interpolation< Array1, T > Renderer;
3427  return( collision_detection( in, out, dmap, Renderer( in, param, table ), param, table, i, j ) );
3428  }
3429 }
3430 
3431 
3449 template < class Array, class T >
3450 volumerender::parameter::vector_type collision_detection( const Array &in, typename Array::size_type image_width, typename Array::size_type image_height, double resoX, double resoY, const volumerender::parameter &param, const volumerender::attribute_table< T > &table, typename Array::size_type i, typename Array::size_type j )
3451 {
3452  return( collision_detection( in, image_width, image_height, resoX, resoY, volumerender::no_depth_map( ), param, table, i, j ) );
3453 }
3454 
3455 
3470 template < class Array1, class Array2, class T >
3471 volumerender::parameter::vector_type collision_detection( const Array1 &in, const Array2 &out, const volumerender::parameter &param, const volumerender::attribute_table< T > &table, typename Array1::size_type i, typename Array1::size_type j )
3472 {
3473  return( collision_detection( in, out, volumerender::no_depth_map( ), param, table, i, j ) );
3474 }
3475 
3476 
3497 template < class Array1, class Array2, class DepthMap, class T >
3498 volumerender::parameter::vector_type collision_detection( const Array1 &in, const Array2 &mk, typename Array1::size_type image_width, typename Array1::size_type image_height, double resoX, double resoY, const DepthMap &dmap, const volumerender::parameter &param, const volumerender::attribute_table< T > &table, const volumerender::attribute_table< T > &mktable, typename Array1::size_type i, typename Array1::size_type j )
3499 {
3500  typedef rendering_helper::value_interpolation_with_mark< Array1, Array2, T > Renderer;
3501  return( collision_detection( in, image_width, image_height, resoX, resoY, dmap, Renderer( in, mk, param, table, mktable ), param, table, i, j ) );
3502 }
3503 
3504 
3522 template < class Array1, class Array2, class Array3, class DepthMap, class T >
3523 volumerender::parameter::vector_type collision_detection( const Array1 &in, const Array2 &mk, const Array3 &out, const DepthMap &dmap, const volumerender::parameter &param, const volumerender::attribute_table< T > &table, const volumerender::attribute_table< T > &mktable, typename Array1::size_type i, typename Array1::size_type j )
3524 {
3525  typedef rendering_helper::value_interpolation_with_mark< Array1, Array2, T > Renderer;
3526  return( collision_detection( in, out, dmap, Renderer( in, mk, param, table, mktable ), param, table, i, j ) );
3527 }
3528 
3529 
3551 template < class Array1, class Array2, class DepthMap, class T >
3552 volumerender::parameter::vector_type collision_detection( const Array1 &in, const Array2 &mk, typename Array1::size_type image_width, typename Array1::size_type image_height, double resoX, double resoY, const DepthMap &dmap, const volumerender::parameter &param, const volumerender::attribute_table< T > &table, const volumerender::attribute_table< T > &mktable, bool apply_and_operation, typename Array1::size_type i, typename Array1::size_type j )
3553 {
3554  if( apply_and_operation )
3555  {
3556  typedef rendering_helper::value_interpolation_and_mark< Array1, Array2, T > Renderer;
3557  return( collision_detection( in, image_width, image_height, resoX, resoY, dmap, Renderer( in, mk, param, table, mktable ), param, table, i, j ) );
3558  }
3559  else
3560  {
3561  typedef rendering_helper::value_interpolation_or_mark< Array1, Array2, T > Renderer;
3562  return( collision_detection( in, image_width, image_height, resoX, resoY, dmap, Renderer( in, mk, param, table, mktable ), param, table, i, j ) );
3563  }
3564 }
3565 
3566 
3585 template < class Array1, class Array2, class Array3, class DepthMap, class T >
3586 volumerender::parameter::vector_type collision_detection( const Array1 &in, const Array2 &mk, const Array3 &out, const DepthMap &dmap, const volumerender::parameter &param, const volumerender::attribute_table< T > &table, const volumerender::attribute_table< T > &mktable, bool apply_and_operation, typename Array1::size_type i, typename Array1::size_type j )
3587 {
3588  if( apply_and_operation )
3589  {
3590  typedef rendering_helper::value_interpolation_and_mark< Array1, Array2, T > Renderer;
3591  return( collision_detection( in, out, dmap, Renderer( in, mk, param, table, mktable ), param, table, i, j ) );
3592  }
3593  else
3594  {
3595  typedef rendering_helper::value_interpolation_or_mark< Array1, Array2, T > Renderer;
3596  return( collision_detection( in, out, dmap, Renderer( in, mk, param, table, mktable ), param, table, i, j ) );
3597  }
3598 }
3599 
3600 
3620 template < class Array1, class Array2, class T >
3621 volumerender::parameter::vector_type collision_detection( const Array1 &in, const Array2 &mk, typename Array1::size_type image_width, typename Array1::size_type image_height, double resoX, double resoY, const volumerender::parameter &param, const volumerender::attribute_table< T > &table, const volumerender::attribute_table< T > &mktable, typename Array1::size_type i, typename Array1::size_type j )
3622 {
3623  return( collision_detection( in, mk, image_width, image_height, resoX, resoY, volumerender::no_depth_map( ), param, table, mktable, i, j ) );
3624 }
3625 
3626 
3643 template < class Array1, class Array2, class Array3, class T >
3644 volumerender::parameter::vector_type collision_detection( const Array1 &in, const Array2 &mk, const Array3 &out, const volumerender::parameter &param, const volumerender::attribute_table< T > &table, const volumerender::attribute_table< T > &mktable, typename Array1::size_type i, typename Array1::size_type j )
3645 {
3646  return( collision_detection( in, mk, out, volumerender::no_depth_map( ), param, table, mktable, i, j ) );
3647 }
3648 
3649 
3670 template < class Array1, class Array2, class T >
3671 volumerender::parameter::vector_type collision_detection( const Array1 &in, const Array2 &mk, typename Array1::size_type image_width, typename Array1::size_type image_height, double resoX, double resoY, const volumerender::parameter &param, const volumerender::attribute_table< T > &table, const volumerender::attribute_table< T > &mktable, bool apply_and_operation, typename Array1::size_type i, typename Array1::size_type j )
3672 {
3673  return( collision_detection( in, mk, image_width, image_height, resoX, resoY, volumerender::no_depth_map( ), param, table, mktable, apply_and_operation, i, j ) );
3674 }
3675 
3676 
3694 template < class Array1, class Array2, class Array3, class T >
3695 volumerender::parameter::vector_type collision_detection( const Array1 &in, const Array2 &mk, const Array3 &out, const volumerender::parameter &param, const volumerender::attribute_table< T > &table, const volumerender::attribute_table< T > &mktable, bool apply_and_operation, typename Array1::size_type i, typename Array1::size_type j )
3696 {
3697  return( collision_detection( in, mk, out, volumerender::no_depth_map( ), param, table, mktable, apply_and_operation, i, j ) );
3698 }
3699 
3700 
3701 
3715 template < class Array, class DepthMap, class ATTRIBUTETYPE >
3716 bool generate_depth_map( const Array &in, DepthMap &dmap, const volumerender::attribute_table< ATTRIBUTETYPE > &table, typename Array::size_type thread_num = 0 )
3717 {
3718  if( is_same_object( in, dmap ) || in.empty( ) || in.depth( ) < 2 )
3719  {
3720  return( false );
3721  }
3722 
3723  typedef typename DepthMap::value_type value_type;
3724  typedef typename Array::size_type size_type;
3725  typedef typename Array::difference_type difference_type;
3726  typedef typename Array::const_pointer const_pointer;
3727 
3728  size_type num = 4;
3729  size_type w = in.width( ) / num;
3730  size_type h = in.height( ) / num;
3731  size_type d = in.depth( ) / num;
3732  size_type rw = in.width( ) > w * num ? 1 : 0;
3733  size_type rh = in.height( ) > h * num ? 1 : 0;
3734  size_type rd = in.depth( ) > d * num ? 1 : 0;
3735 
3736  dmap.resize( w + rw + 1, h + rh + 1, d + rd + 1 );
3737  dmap.reso( 1.0, 1.0, 1.0 );
3738  dmap.fill( 1 );
3739 
3740  for( size_type k = 0 ; k < in.depth( ) ; k++ )
3741  {
3742  size_type _3 = k / num;
3743  for( size_type j = 0 ; j < in.height( ) ; j++ )
3744  {
3745  size_type _2 = j / num;
3746  for( size_type i = 0 ; i < in.width( ) ; i++ )
3747  {
3748  size_type _1 = i / num;
3749 
3750  if( table.has_alpha( in( i, j, k ) ) )
3751  {
3752  dmap( _1, _2, _3 ) = 0;
3753  }
3754  }
3755  }
3756  }
3757 
3758  euclidean::distance_transform( dmap, dmap, thread_num );
3759 
3760  for( size_type i = 0 ; i < dmap.size( ) ; i++ )
3761  {
3762  dmap[ i ] = static_cast< value_type >( std::sqrt( static_cast< double >( dmap[ i ] ) ) );
3763  }
3764 
3765  return( true );
3766 }
3767 
3768 
3770 // ボリュームレンダリンググループの終わり
3771 
3773 // 可視化グループの終わり
3774 
3775 
3776 // mist名前空間の終わり
3777 _MIST_END
3778 
3779 
3780 #endif // __INCLUDE_MIST_VOLUMERENDER__
3781 
Generated on Wed Nov 12 2014 19:44:27 for MIST by doxygen 1.8.1.2