- メインページ
- モジュール
- ネームスペース
- データ構造
- ファイル

- ファイル一覧
- グローバル

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