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

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

registration.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
33 #ifndef __INCLUDE_MIST_REGISTRATION__
34 #define __INCLUDE_MIST_REGISTRATION__
35
36
37 #ifndef __INCLUDE_MIST_CONF_H__
38 #include "config/mist_conf.h"
39 #endif
40
41
42 #ifndef __INCLUDE_MIST_H__
43 #include "mist.h"
44 #endif
45
46 #ifndef __INCLUDE_MIST_VECTOR__
47 #include "vector.h"
48 #endif
49
50 #ifndef __INCLUDE_MIST_MINIMIZATION__
51 #include "minimization.h"
52 #endif
53
54 #ifndef __INCLUDE_MIST_DRAWING__
55 #include "drawing.h"
56 #endif
57
58 #ifndef __INCLUDE_MIST_LIMITS__
59 #include "limits.h"
60 #endif
61
62 #ifndef __INCLUDE_MIST_FFT__
63 #include "fft/fft.h"
64 #endif
65
66 #ifndef __INCLUDE_CONVERTER__
67 #include "converter.h"
68 #endif
69
70 #ifndef __INCLUDE_INTERPOLATE__
71 #include "interpolate.h"
72 #endif
73
74 #ifndef __INCLUDE_GEOMETRY__
75 #include "geometry.h"
76 #endif
77
78 #ifndef M_PI
79 #define M_PI 3.14159265358979323846 /* pi */
80 #endif
81
82 #include <vector>
83 #include <complex>
84
85
86 // mist名前空間の始まり
87 _MIST_BEGIN
88
89 namespace __non_rigid_registration_utility__
90 {
91 inline void FFD( double v, double &f0, double &f1, double &f2, double &f3 )
92 {
93 double v_ = 1.0 - v;
94 double v2 = v * v;
95 double v3 = v2 * v;
96 f0 = v_ * v_ * v_ / 6.0;
97 f1 = 0.5 * v3 - v2 + 4.0 / 6.0;
98 f2 = ( v2 - v3 + v ) * 0.5 + 1.0 / 6.0;
99 f3 = v3 / 6.0;
100 }
101
102
103 // 制御点情報を用いて，ソース画像をターゲット画像に変形する
104 template < class TARGETTYPE, class TARGETTYPEA, class SOURCETYPE, class SOURCETYPEA, class CONTROLMESH >
105 static void non_rigid_transformation_( array2< TARGETTYPE, TARGETTYPEA > &target, const array2< SOURCETYPE, SOURCETYPEA > &source, const CONTROLMESH &control_mesh,
106 typename CONTROLMESH::difference_type mx = -1,
107 typename CONTROLMESH::difference_type my = -1,
108 typename CONTROLMESH::difference_type mz = -1,
109 typename CONTROLMESH::size_type thread_id = 0,
110 typename CONTROLMESH::size_type thread_num = 1
111 )
112 {
113 typedef array2< TARGETTYPE, TARGETTYPEA > target_image_type;
114 typedef array2< SOURCETYPE, SOURCETYPEA > source_image_type;
115 typedef CONTROLMESH control_mesh_type;
116
117 typedef typename target_image_type::size_type size_type;
118 typedef typename target_image_type::difference_type difference_type;
119 typedef typename target_image_type::value_type value_type;
120 value_type minimum = type_limits< value_type >::minimum( );
121
122 difference_type i, j;
123 double u, v;
124 double sx[ 4 ], sy[ 4 ];
125 double stepW = control_mesh.width( ) == 1 ? 1.0 : target.width( ) / static_cast< double >( control_mesh.width( ) - 1 );
126 double stepH = control_mesh.height( ) == 1 ? 1.0 : target.height( ) / static_cast< double >( control_mesh.height( ) - 1 );
127 double _1_stepW = 1.0 / stepW;
128 double _1_stepH = 1.0 / stepH;
129 double _1_ax = 1.0 / source.reso1( );
130 double _1_ay = 1.0 / source.reso2( );
131
132
133 difference_type d0, d1, d2, d3;
134 {
135 difference_type cx = source.width( ) / 2;
136 difference_type cy = source.height( ) / 2;
137 typename source_image_type::const_pointer ppp = &source( cx, cy );
138 d0 = 0;
139 d1 = &source( cx , cy + 1 ) - ppp;
140 d2 = &source( cx + 1, cy + 1 ) - ppp;
141 d3 = &source( cx + 1, cy ) - ppp;
142 }
143
144 difference_type isx, isy, iex, iey;
145 difference_type mw = control_mesh.width( );
146 difference_type mh = control_mesh.height( );
147 if( mx < 0 || mw <= mx || my < 0 || mh <= my )
148 {
149 isx = 0;
150 iex = target.width( ) - 1;
151 isy = 0;
152 iey = target.height( ) - 1;
153 }
154 else
155 {
156 isx = static_cast< difference_type >( ( mx - 2 ) * stepW );
157 iex = static_cast< difference_type >( ( mx + 2 ) * stepW );
158 isy = static_cast< difference_type >( ( my - 2 ) * stepH );
159 iey = static_cast< difference_type >( ( my + 2 ) * stepH );
160
161 difference_type tw = target.width( );
162 difference_type th = target.height( );
163 isx = isx > 0 ? isx : 0;
164 isy = isy > 0 ? isy : 0;
165 iex = iex < tw ? iex : tw - 1;
166 iey = iey < th ? iey : th - 1;
167 }
168
169
170 for( difference_type y = isy + thread_id ; y <= iey ; y += thread_num )
171 {
172 v = y * _1_stepH;
173 j = static_cast< difference_type >( v );
174 v -= j;
175 j--;
176
177 FFD( v, sy[ 0 ], sy[ 1 ], sy[ 2 ], sy[ 3 ] );
178
179 for( difference_type x = isx ; x <= iex ; x++ )
180 {
181 u = x * _1_stepW;
182 i = static_cast< difference_type >( u );
183 u -= i;
184 i--;
185
186 FFD( u, sx[ 0 ], sx[ 1 ], sx[ 2 ], sx[ 3 ] );
187
188 double xx = 0.0, yy = 0.0;
189 for( difference_type m = 0 ; m <= 3 ; m++ )
190 {
191 typename control_mesh_type::const_pointer p = &control_mesh( i, j + m, 0 );
192 for( difference_type l = 0 ; l <= 3 ; l++ )
193 {
194 double val = sx[ l ] * sy[ m ];
195 xx += val * p[ l ].x;
196 yy += val * p[ l ].y;
197 }
198 }
199
200 xx *= _1_ax;
201 yy *= _1_ay;
202
203 if( 0 <= xx && xx < source.width( ) - 1 && 0 <= yy && yy < source.height( ) - 1 )
204 {
205 difference_type ixx = static_cast< size_type >( xx );
206 difference_type iyy = static_cast< size_type >( yy );
207 typename source_image_type::const_pointer p = &source( ixx, iyy, 0 );
208 xx -= ixx;
209 yy -= iyy;
210
211 target( x, y ) = static_cast< typename target_image_type::value_type >( ( p[ d0 ] + ( p[ d3 ] - p[ d0 ] ) * xx ) + ( p[ d1 ] - p[ d0 ] + ( p[ d0 ] - p[ d1 ] + p[ d2 ] - p[ d3 ] ) * xx ) * yy );
212 }
213 else
214 {
215 target( x, y ) = static_cast< typename target_image_type::value_type >( minimum );
216 }
217 }
218 }
219 }
220
221
222 // 制御点情報を用いて，ソース画像をターゲット画像に変形する
223 template < class TARGETTYPE, class TARGETTYPEA, class SOURCETYPE, class SOURCETYPEA, class CONTROLMESH >
224 static void non_rigid_transformation_( array3< TARGETTYPE, TARGETTYPEA > &target, const array3< SOURCETYPE, SOURCETYPEA > &source, const CONTROLMESH &control_mesh,
225 typename CONTROLMESH::difference_type mx = -1,
226 typename CONTROLMESH::difference_type my = -1,
227 typename CONTROLMESH::difference_type mz = -1,
228 typename CONTROLMESH::size_type thread_id = 0,
229 typename CONTROLMESH::size_type thread_num = 1
230 )
231 {
232 typedef array3< TARGETTYPE, TARGETTYPEA > target_image_type;
233 typedef array3< SOURCETYPE, SOURCETYPEA > source_image_type;
234 typedef CONTROLMESH control_mesh_type;
235
236 typedef typename target_image_type::size_type size_type;
237 typedef typename target_image_type::difference_type difference_type;
238
239 typedef typename target_image_type::value_type value_type;
240 value_type minimum = type_limits< value_type >::minimum( );
241
242 difference_type i, j, k;
243 double u, v, w;
244 double sx[ 4 ], sy[ 4 ], sz[ 4 ];
245 double stepW = control_mesh.width( ) == 1 ? 1.0 : target.width( ) / static_cast< double >( control_mesh.width( ) - 1 );
246 double stepH = control_mesh.height( ) == 1 ? 1.0 : target.height( ) / static_cast< double >( control_mesh.height( ) - 1 );
247 double stepD = control_mesh.depth( ) == 1 ? 1.0 : target.depth( ) / static_cast< double >( control_mesh.depth( ) - 1 );
248 double _1_stepW = 1.0 / stepW;
249 double _1_stepH = 1.0 / stepH;
250 double _1_stepD = 1.0 / stepD;
251 double _1_ax = 1.0 / source.reso1( );
252 double _1_ay = 1.0 / source.reso2( );
253 double _1_az = 1.0 / source.reso3( );
254
255 difference_type d0, d1, d2, d3, d4, d5, d6, d7;
256 {
257 difference_type cx = source.width( ) / 2;
258 difference_type cy = source.height( ) / 2;
259 difference_type cz = source.depth( ) / 2;
260 typename source_image_type::const_pointer ppp = &source( cx, cy, cz );
261 d0 = 0;
262 d1 = &source( cx , cy + 1, cz ) - ppp;
263 d2 = &source( cx + 1, cy + 1, cz ) - ppp;
264 d3 = &source( cx + 1, cy , cz ) - ppp;
265 d4 = &source( cx , cy , cz + 1 ) - ppp;
266 d5 = &source( cx , cy + 1, cz + 1 ) - ppp;
267 d6 = &source( cx + 1, cy + 1, cz + 1 ) - ppp;
268 d7 = &source( cx + 1, cy , cz + 1 ) - ppp;
269 }
270
271 difference_type isx, isy, isz, iex, iey, iez;
272 difference_type mw = control_mesh.width( );
273 difference_type mh = control_mesh.height( );
274 difference_type md = control_mesh.depth( );
275 if( mx < 0 || mw <= mx || my < 0 || mh <= my || mz < 0 || md <= mz )
276 {
277 isx = 0;
278 iex = target.width( ) - 1;
279 isy = 0;
280 iey = target.height( ) - 1;
281 isz = 0;
282 iez = target.depth( ) - 1;
283 }
284 else
285 {
286 isx = static_cast< difference_type >( ( mx - 2 ) * stepW );
287 iex = static_cast< difference_type >( ( mx + 2 ) * stepW );
288 isy = static_cast< difference_type >( ( my - 2 ) * stepH );
289 iey = static_cast< difference_type >( ( my + 2 ) * stepH );
290 isz = static_cast< difference_type >( ( mz - 2 ) * stepD );
291 iez = static_cast< difference_type >( ( mz + 2 ) * stepD );
292
293 difference_type tw = target.width( );
294 difference_type th = target.height( );
295 difference_type td = target.depth( );
296 isx = isx > 0 ? isx : 0;
297 isy = isy > 0 ? isy : 0;
298 isz = isz > 0 ? isz : 0;
299 iex = iex < tw ? iex : tw - 1;
300 iey = iey < th ? iey : th - 1;
301 iez = iez < td ? iez : td - 1;
302 }
303
304
305 for( difference_type z = isz + thread_id ; z <= iez ; z += thread_num )
306 {
307 w = z * _1_stepD;
308 k = static_cast< difference_type >( w );
309 w -= k;
310 k--;
311
312 FFD( w, sz[ 0 ], sz[ 1 ], sz[ 2 ], sz[ 3 ] );
313
314 for( difference_type y = isy ; y <= iey ; y++ )
315 {
316 v = y * _1_stepH;
317 j = static_cast< difference_type >( v );
318 v -= j;
319 j--;
320
321 FFD( v, sy[ 0 ], sy[ 1 ], sy[ 2 ], sy[ 3 ] );
322
323 for( difference_type x = isx ; x <= iex ; x++ )
324 {
325 u = x * _1_stepW;
326 i = static_cast< difference_type >( u );
327 u -= i;
328 i--;
329
330 FFD( u, sx[ 0 ], sx[ 1 ], sx[ 2 ], sx[ 3 ] );
331
332 double xx = 0.0, yy = 0.0, zz = 0.0;
333 for( difference_type n = 0 ; n <= 3 ; n++ )
334 {
335 for( difference_type m = 0 ; m <= 3 ; m++ )
336 {
337 typename control_mesh_type::const_pointer p = &control_mesh( i, j + m, k + n );
338 double vvv = sy[ m ] * sz[ n ];
339 for( difference_type l = 0 ; l <= 3 ; l++ )
340 {
341 double vv = sx[ l ] * vvv;
342 xx += vv * p[ l ].x;
343 yy += vv * p[ l ].y;
344 zz += vv * p[ l ].z;
345 }
346 }
347 }
348
349 xx *= _1_ax;
350 yy *= _1_ay;
351 zz *= _1_az;
352
353 if( 0 <= xx && xx < source.width( ) - 1 && 0 <= yy && yy < source.height( ) - 1 && 0 <= zz && zz < source.depth( ) - 1 )
354 {
355 difference_type ixx = static_cast< size_type >( xx );
356 difference_type iyy = static_cast< size_type >( yy );
357 difference_type izz = static_cast< size_type >( zz );
358 typename source_image_type::const_pointer p = &source( ixx, iyy, izz );
359 xx -= ixx;
360 yy -= iyy;
361 zz -= izz;
362
363 double ct = ( p[ d0 ] + ( p[ d3 ] - p[ d0 ] ) * xx ) + ( p[ d1 ] - p[ d0 ] + ( p[ d0 ] - p[ d1 ] + p[ d2 ] - p[ d3 ] ) * xx ) * yy;
364 ct += ( ( p[ d4 ] + ( p[ d7 ] - p[ d4 ] ) * xx ) + ( p[ d5 ] - p[ d4 ] + ( p[ d4 ] - p[ d5 ] + p[ d6 ] - p[ d7 ] ) * xx ) * yy - ct ) * zz;
365 target( x, y, z ) = static_cast< typename target_image_type::value_type >( ct );
366 }
367 else
368 {
369 target( x, y, z ) = static_cast< typename target_image_type::value_type >( minimum );
370 }
371 }
372 }
373 }
374 }
375
376 struct ffd_coefficient
377 {
378 struct coefficient
379 {
380 double value[ 4 ];
381 double &operator []( size_t index ){ return( value[ index ] ); }
382 const double &operator []( size_t index ) const { return( value[ index ] ); }
383 };
384
385 struct coefficients
386 {
387 double value[ 16 ];
388 double &operator []( size_t index ){ return( value[ index ] ); }
389 const double &operator []( size_t index ) const { return( value[ index ] ); }
390 };
391
392 typedef size_t size_type;
393 typedef ptrdiff_t difference_type;
394 typedef mist::array< coefficient > coefficient_type;
395
396 coefficient_type coeff_z;
397 mist::array2< coefficients > coeff_xy;
398
399 template < class CONTROLMESHTYPE >
400 void initialize( size_type width, size_type height, size_type depth, const CONTROLMESHTYPE &control_mesh )
401 {
402 coeff_z.resize( depth );
403 coeff_xy.resize( width, height );
404
405 double stepW = control_mesh.width( ) == 1 ? 1.0 : width / static_cast< double >( control_mesh.width( ) - 1 );
406 double stepH = control_mesh.height( ) == 1 ? 1.0 : height / static_cast< double >( control_mesh.height( ) - 1 );
407 double stepD = control_mesh.depth( ) == 1 ? 1.0 : depth / static_cast< double >( control_mesh.depth( ) - 1 );
408 double _1_stepW = 1.0 / stepW;
409 double _1_stepH = 1.0 / stepH;
410 double _1_stepD = 1.0 / stepD;
411
412 for( size_type z = 0 ; z < coeff_z.size( ) ; z++ )
413 {
414 double w = z * _1_stepD;
415 w -= static_cast< difference_type >( w );
416
417 FFD( w, coeff_z[ z ][ 0 ], coeff_z[ z ][ 1 ], coeff_z[ z ][ 2 ], coeff_z[ z ][ 3 ] );
418 }
419
420 double coeff_x[ 4 ], coeff_y[ 4 ];
421 for( size_type y = 0 ; y < height ; y++ )
422 {
423 double v = y * _1_stepH;
424 v -= static_cast< difference_type >( v );
425
426 FFD( v, coeff_y[ 0 ], coeff_y[ 1 ], coeff_y[ 2 ], coeff_y[ 3 ] );
427
428 for( size_type x = 0 ; x < width ; x++ )
429 {
430 double u = x * _1_stepW;
431 u -= static_cast< difference_type >( u );
432
433 FFD( u, coeff_x[ 0 ], coeff_x[ 1 ], coeff_x[ 2 ], coeff_x[ 3 ] );
434
435 coeff_xy( x, y )[ 0 ] = coeff_x[ 0 ] * coeff_y[ 0 ];
436 coeff_xy( x, y )[ 1 ] = coeff_x[ 1 ] * coeff_y[ 0 ];
437 coeff_xy( x, y )[ 2 ] = coeff_x[ 2 ] * coeff_y[ 0 ];
438 coeff_xy( x, y )[ 3 ] = coeff_x[ 3 ] * coeff_y[ 0 ];
439 coeff_xy( x, y )[ 4 ] = coeff_x[ 0 ] * coeff_y[ 1 ];
440 coeff_xy( x, y )[ 5 ] = coeff_x[ 1 ] * coeff_y[ 1 ];
441 coeff_xy( x, y )[ 6 ] = coeff_x[ 2 ] * coeff_y[ 1 ];
442 coeff_xy( x, y )[ 7 ] = coeff_x[ 3 ] * coeff_y[ 1 ];
443 coeff_xy( x, y )[ 8 ] = coeff_x[ 0 ] * coeff_y[ 2 ];
444 coeff_xy( x, y )[ 9 ] = coeff_x[ 1 ] * coeff_y[ 2 ];
445 coeff_xy( x, y )[ 10 ] = coeff_x[ 2 ] * coeff_y[ 2 ];
446 coeff_xy( x, y )[ 11 ] = coeff_x[ 3 ] * coeff_y[ 2 ];
447 coeff_xy( x, y )[ 12 ] = coeff_x[ 0 ] * coeff_y[ 3 ];
448 coeff_xy( x, y )[ 13 ] = coeff_x[ 1 ] * coeff_y[ 3 ];
449 coeff_xy( x, y )[ 14 ] = coeff_x[ 2 ] * coeff_y[ 3 ];
450 coeff_xy( x, y )[ 15 ] = coeff_x[ 3 ] * coeff_y[ 3 ];
451 }
452 }
453 }
454 };
455
456
457
458 // 制御点情報を用いて，ソース画像をターゲット画像に変形する
459 template < class TARGETTYPE, class TARGETTYPEA, class SOURCETYPE, class SOURCETYPEA, class CONTROLMESH >
460 static void non_rigid_transformation_( array2< TARGETTYPE, TARGETTYPEA > &target, const array2< SOURCETYPE, SOURCETYPEA > &source, const CONTROLMESH &control_mesh,
461 const ffd_coefficient &ffd_coeff,
462 typename CONTROLMESH::difference_type mx = -1,
463 typename CONTROLMESH::difference_type my = -1,
464 typename CONTROLMESH::difference_type mz = -1,
465 typename CONTROLMESH::size_type thread_id = 0,
466 typename CONTROLMESH::size_type thread_num = 1
467 )
468 {
469 typedef array2< TARGETTYPE, TARGETTYPEA > target_image_type;
470 typedef array2< SOURCETYPE, SOURCETYPEA > source_image_type;
471 typedef CONTROLMESH control_mesh_type;
472
473 typedef typename target_image_type::size_type size_type;
474 typedef typename target_image_type::difference_type difference_type;
475 typedef typename target_image_type::value_type value_type;
476 value_type minimum = type_limits< value_type >::minimum( );
477
478 difference_type i, j;
479 double stepW = control_mesh.width( ) == 1 ? 1.0 : target.width( ) / static_cast< double >( control_mesh.width( ) - 1 );
480 double stepH = control_mesh.height( ) == 1 ? 1.0 : target.height( ) / static_cast< double >( control_mesh.height( ) - 1 );
481 double _1_stepW = 1.0 / stepW;
482 double _1_stepH = 1.0 / stepH;
483 double _1_ax = 1.0 / source.reso1( );
484 double _1_ay = 1.0 / source.reso2( );
485
486 difference_type d0, d1, d2, d3, dm;
487 {
488 difference_type cx = source.width( ) / 2;
489 difference_type cy = source.height( ) / 2;
490 typename source_image_type::const_pointer ppp = &source( cx, cy );
491 d0 = 0;
492 d1 = &source( cx , cy + 1 ) - ppp;
493 d2 = &source( cx + 1, cy + 1 ) - ppp;
494 d3 = &source( cx + 1, cy ) - ppp;
495 dm = &control_mesh( 0, 1, 0 ) - &control_mesh( 0, 0, 0 );
496 }
497
498 difference_type isx, isy, iex, iey;
499 difference_type mw = control_mesh.width( );
500 difference_type mh = control_mesh.height( );
501 if( mx < 0 || mw <= mx || my < 0 || mh <= my )
502 {
503 isx = 0;
504 iex = target.width( ) - 1;
505 isy = 0;
506 iey = target.height( ) - 1;
507 }
508 else
509 {
510 isx = static_cast< difference_type >( ( mx - 2 ) * stepW );
511 iex = static_cast< difference_type >( ( mx + 2 ) * stepW );
512 isy = static_cast< difference_type >( ( my - 2 ) * stepH );
513 iey = static_cast< difference_type >( ( my + 2 ) * stepH );
514
515 difference_type tw = target.width( );
516 difference_type th = target.height( );
517 isx = isx > 0 ? isx : 0;
518 isy = isy > 0 ? isy : 0;
519 iex = iex < tw ? iex : tw - 1;
520 iey = iey < th ? iey : th - 1;
521 }
522
523
524 for( difference_type y = isy + thread_id ; y <= iey ; y += thread_num )
525 {
526 j = static_cast< difference_type >( y * _1_stepH ) - 1;
527
528 for( difference_type x = isx ; x <= iex ; x++ )
529 {
530 i = static_cast< difference_type >( x * _1_stepW ) - 1;
531
532 const ffd_coefficient::coefficients &coeff = ffd_coeff.coeff_xy( x, y );
533
534 double xx = 0.0, yy = 0.0;
535 typename control_mesh_type::const_pointer p = &control_mesh( i, j, 0 );
536
537 xx += coeff[ 0 ] * p[ 0 ].x; yy += coeff[ 0 ] * p[ 0 ].y;
538 xx += coeff[ 1 ] * p[ 1 ].x; yy += coeff[ 1 ] * p[ 1 ].y;
539 xx += coeff[ 2 ] * p[ 2 ].x; yy += coeff[ 2 ] * p[ 2 ].y;
540 xx += coeff[ 3 ] * p[ 3 ].x; yy += coeff[ 3 ] * p[ 3 ].y;
541
542 p += dm;
543
544 xx += coeff[ 4 ] * p[ 0 ].x; yy += coeff[ 4 ] * p[ 0 ].y;
545 xx += coeff[ 5 ] * p[ 1 ].x; yy += coeff[ 5 ] * p[ 1 ].y;
546 xx += coeff[ 6 ] * p[ 2 ].x; yy += coeff[ 6 ] * p[ 2 ].y;
547 xx += coeff[ 7 ] * p[ 3 ].x; yy += coeff[ 7 ] * p[ 3 ].y;
548
549 p += dm;
550
551 xx += coeff[ 8 ] * p[ 0 ].x; yy += coeff[ 8 ] * p[ 0 ].y;
552 xx += coeff[ 9 ] * p[ 1 ].x; yy += coeff[ 9 ] * p[ 1 ].y;
553 xx += coeff[ 10 ] * p[ 2 ].x; yy += coeff[ 10 ] * p[ 2 ].y;
554 xx += coeff[ 11 ] * p[ 3 ].x; yy += coeff[ 11 ] * p[ 3 ].y;
555
556 p += dm;
557
558 xx += coeff[ 12 ] * p[ 0 ].x; yy += coeff[ 12 ] * p[ 0 ].y;
559 xx += coeff[ 13 ] * p[ 1 ].x; yy += coeff[ 13 ] * p[ 1 ].y;
560 xx += coeff[ 14 ] * p[ 2 ].x; yy += coeff[ 14 ] * p[ 2 ].y;
561 xx += coeff[ 15 ] * p[ 3 ].x; yy += coeff[ 15 ] * p[ 3 ].y;
562
563 xx *= _1_ax;
564 yy *= _1_ay;
565
566 if( 0 <= xx && xx < source.width( ) - 1 && 0 <= yy && yy < source.height( ) - 1 )
567 {
568 difference_type ixx = static_cast< size_type >( xx );
569 difference_type iyy = static_cast< size_type >( yy );
570 typename source_image_type::const_pointer p = &source( ixx, iyy, 0 );
571 xx -= ixx;
572 yy -= iyy;
573
574 target( x, y ) = static_cast< typename target_image_type::value_type >( ( p[ d0 ] + ( p[ d3 ] - p[ d0 ] ) * xx ) + ( p[ d1 ] - p[ d0 ] + ( p[ d0 ] - p[ d1 ] + p[ d2 ] - p[ d3 ] ) * xx ) * yy );
575 }
576 else
577 {
578 target( x, y ) = static_cast< typename target_image_type::value_type >( minimum );
579 }
580 }
581 }
582 }
583
584
585
586 // 制御点情報を用いて，ソース画像をターゲット画像に変形する
587 template < class TARGETTYPE, class TARGETTYPEA, class SOURCETYPE, class SOURCETYPEA, class CONTROLMESH >
588 static void non_rigid_transformation_( array3< TARGETTYPE, TARGETTYPEA > &target, const array3< SOURCETYPE, SOURCETYPEA > &source, const CONTROLMESH &control_mesh,
589 const ffd_coefficient &ffd_coeff,
590 typename CONTROLMESH::difference_type mx = -1,
591 typename CONTROLMESH::difference_type my = -1,
592 typename CONTROLMESH::difference_type mz = -1,
593 typename CONTROLMESH::size_type thread_id = 0,
594 typename CONTROLMESH::size_type thread_num = 1
595 )
596 {
597 typedef array3< TARGETTYPE, TARGETTYPEA > target_image_type;
598 typedef array3< SOURCETYPE, SOURCETYPEA > source_image_type;
599 typedef CONTROLMESH control_mesh_type;
600
601 typedef typename target_image_type::size_type size_type;
602 typedef typename target_image_type::difference_type difference_type;
603 typedef typename target_image_type::value_type value_type;
604 value_type minimum = type_limits< value_type >::minimum( );
605
606 difference_type i, j, k;
607 double sz[ 4 ];
608 double stepW = control_mesh.width( ) == 1 ? 1.0 : target.width( ) / static_cast< double >( control_mesh.width( ) - 1 );
609 double stepH = control_mesh.height( ) == 1 ? 1.0 : target.height( ) / static_cast< double >( control_mesh.height( ) - 1 );
610 double stepD = control_mesh.depth( ) == 1 ? 1.0 : target.depth( ) / static_cast< double >( control_mesh.depth( ) - 1 );
611 double _1_stepW = 1.0 / stepW;
612 double _1_stepH = 1.0 / stepH;
613 double _1_stepD = 1.0 / stepD;
614 double _1_ax = 1.0 / source.reso1( );
615 double _1_ay = 1.0 / source.reso2( );
616 double _1_az = 1.0 / source.reso3( );
617
618 difference_type d0, d1, d2, d3, d4, d5, d6, d7, dm;
619 {
620 difference_type cx = source.width( ) / 2;
621 difference_type cy = source.height( ) / 2;
622 difference_type cz = source.depth( ) / 2;
623 typename source_image_type::const_pointer ppp = &source( cx, cy, cz );
624 d0 = 0;
625 d1 = &source( cx , cy + 1, cz ) - ppp;
626 d2 = &source( cx + 1, cy + 1, cz ) - ppp;
627 d3 = &source( cx + 1, cy , cz ) - ppp;
628 d4 = &source( cx , cy , cz + 1 ) - ppp;
629 d5 = &source( cx , cy + 1, cz + 1 ) - ppp;
630 d6 = &source( cx + 1, cy + 1, cz + 1 ) - ppp;
631 d7 = &source( cx + 1, cy , cz + 1 ) - ppp;
632 dm = &control_mesh( 0, 1, 0 ) - &control_mesh( 0, 0, 0 );
633 }
634
635 difference_type isx, isy, isz, iex, iey, iez;
636 difference_type mw = control_mesh.width( );
637 difference_type mh = control_mesh.height( );
638 difference_type md = control_mesh.depth( );
639 if( mx < 0 || mw <= mx || my < 0 || mh <= my || mz < 0 || md <= mz )
640 {
641 isx = 0;
642 iex = target.width( ) - 1;
643 isy = 0;
644 iey = target.height( ) - 1;
645 isz = 0;
646 iez = target.depth( ) - 1;
647 }
648 else
649 {
650 isx = static_cast< difference_type >( ( mx - 2 ) * stepW );
651 iex = static_cast< difference_type >( ( mx + 2 ) * stepW );
652 isy = static_cast< difference_type >( ( my - 2 ) * stepH );
653 iey = static_cast< difference_type >( ( my + 2 ) * stepH );
654 isz = static_cast< difference_type >( ( mz - 2 ) * stepD );
655 iez = static_cast< difference_type >( ( mz + 2 ) * stepD );
656
657 difference_type tw = target.width( );
658 difference_type th = target.height( );
659 difference_type td = target.depth( );
660 isx = isx > 0 ? isx : 0;
661 isy = isy > 0 ? isy : 0;
662 isz = isz > 0 ? isz : 0;
663 iex = iex < tw ? iex : tw - 1;
664 iey = iey < th ? iey : th - 1;
665 iez = iez < td ? iez : td - 1;
666 }
667
668 typedef ffd_coefficient::coefficient_type coefficient_type;
669 const coefficient_type &ffd_coeff_z = ffd_coeff.coeff_z;
670
671 for( difference_type z = isz + thread_id ; z <= iez ; z += thread_num )
672 {
673 k = static_cast< difference_type >( z * _1_stepD ) - 1;
674
675 sz[ 0 ] = ffd_coeff_z[ z ][ 0 ];
676 sz[ 1 ] = ffd_coeff_z[ z ][ 1 ];
677 sz[ 2 ] = ffd_coeff_z[ z ][ 2 ];
678 sz[ 3 ] = ffd_coeff_z[ z ][ 3 ];
679
680 for( difference_type y = isy ; y <= iey ; y++ )
681 {
682 j = static_cast< difference_type >( y * _1_stepH ) - 1;
683
684 for( difference_type x = isx ; x <= iex ; x++ )
685 {
686 i = static_cast< difference_type >( x * _1_stepW ) - 1;
687
688 double xx = 0.0, yy = 0.0, zz = 0.0, val;
689 const ffd_coefficient::coefficients &coeff = ffd_coeff.coeff_xy( x, y );
690
691 for( difference_type n = 0 ; n <= 3 ; n++ )
692 {
693 typename control_mesh_type::const_pointer p = &control_mesh( i, j, k + n );
694
695 val = sz[ n ] * coeff[ 0 ]; xx += val * p[ 0 ].x; yy += val * p[ 0 ].y; zz += val * p[ 0 ].z;
696 val = sz[ n ] * coeff[ 1 ]; xx += val * p[ 1 ].x; yy += val * p[ 1 ].y; zz += val * p[ 1 ].z;
697 val = sz[ n ] * coeff[ 2 ]; xx += val * p[ 2 ].x; yy += val * p[ 2 ].y; zz += val * p[ 2 ].z;
698 val = sz[ n ] * coeff[ 3 ]; xx += val * p[ 3 ].x; yy += val * p[ 3 ].y; zz += val * p[ 3 ].z;
699
700 p += dm;
701
702 val = sz[ n ] * coeff[ 4 ]; xx += val * p[ 0 ].x; yy += val * p[ 0 ].y; zz += val * p[ 0 ].z;
703 val = sz[ n ] * coeff[ 5 ]; xx += val * p[ 1 ].x; yy += val * p[ 1 ].y; zz += val * p[ 1 ].z;
704 val = sz[ n ] * coeff[ 6 ]; xx += val * p[ 2 ].x; yy += val * p[ 2 ].y; zz += val * p[ 2 ].z;
705 val = sz[ n ] * coeff[ 7 ]; xx += val * p[ 3 ].x; yy += val * p[ 3 ].y; zz += val * p[ 3 ].z;
706
707 p += dm;
708
709 val = sz[ n ] * coeff[ 8 ]; xx += val * p[ 0 ].x; yy += val * p[ 0 ].y; zz += val * p[ 0 ].z;
710 val = sz[ n ] * coeff[ 9 ]; xx += val * p[ 1 ].x; yy += val * p[ 1 ].y; zz += val * p[ 1 ].z;
711 val = sz[ n ] * coeff[ 10 ]; xx += val * p[ 2 ].x; yy += val * p[ 2 ].y; zz += val * p[ 2 ].z;
712 val = sz[ n ] * coeff[ 11 ]; xx += val * p[ 3 ].x; yy += val * p[ 3 ].y; zz += val * p[ 3 ].z;
713
714 p += dm;
715
716 val = sz[ n ] * coeff[ 12 ]; xx += val * p[ 0 ].x; yy += val * p[ 0 ].y; zz += val * p[ 0 ].z;
717 val = sz[ n ] * coeff[ 13 ]; xx += val * p[ 1 ].x; yy += val * p[ 1 ].y; zz += val * p[ 1 ].z;
718 val = sz[ n ] * coeff[ 14 ]; xx += val * p[ 2 ].x; yy += val * p[ 2 ].y; zz += val * p[ 2 ].z;
719 val = sz[ n ] * coeff[ 15 ]; xx += val * p[ 3 ].x; yy += val * p[ 3 ].y; zz += val * p[ 3 ].z;
720 }
721
722 xx *= _1_ax;
723 yy *= _1_ay;
724 zz *= _1_az;
725
726 if( 0 <= xx && xx < source.width( ) - 1 && 0 <= yy && yy < source.height( ) - 1 && 0 <= zz && zz < source.depth( ) - 1 )
727 {
728 difference_type ixx = static_cast< size_type >( xx );
729 difference_type iyy = static_cast< size_type >( yy );
730 difference_type izz = static_cast< size_type >( zz );
731 typename source_image_type::const_pointer p = &source( ixx, iyy, izz );
732 xx -= ixx;
733 yy -= iyy;
734 zz -= izz;
735
736 double ct = ( p[ d0 ] + ( p[ d3 ] - p[ d0 ] ) * xx ) + ( p[ d1 ] - p[ d0 ] + ( p[ d0 ] - p[ d1 ] + p[ d2 ] - p[ d3 ] ) * xx ) * yy;
737 ct += ( ( p[ d4 ] + ( p[ d7 ] - p[ d4 ] ) * xx ) + ( p[ d5 ] - p[ d4 ] + ( p[ d4 ] - p[ d5 ] + p[ d6 ] - p[ d7 ] ) * xx ) * yy - ct ) * zz;
738 target( x, y, z ) = static_cast< typename target_image_type::value_type >( ct );
739 }
740 else
741 {
742 target( x, y, z ) = static_cast< typename target_image_type::value_type >( minimum );
743 }
744 }
745 }
746 }
747 }
748
749
750 template < class TARGETTYPE, class SOURCETYPE, class CONTROLMESH >
751 class non_rigid_registration_thread : public thread< non_rigid_registration_thread< TARGETTYPE, SOURCETYPE, CONTROLMESH > >
752 {
753 public:
754 typedef TARGETTYPE target_image_type;
755 typedef SOURCETYPE source_image_type;
756 typedef CONTROLMESH control_mesh_type;
757 typedef thread< non_rigid_registration_thread< TARGETTYPE, SOURCETYPE, CONTROLMESH > > base;
758 typedef typename base::thread_exit_type thread_exit_type;
759 typedef typename TARGETTYPE::size_type size_type;
760
761 private:
762 size_t thread_id_;
763 size_t thread_num_;
764
765 // 入出力用の画像へのポインタ
766 target_image_type *target_;
767 const source_image_type *source_;
768 const control_mesh_type *control_mesh_;
769 size_type mx_;
770 size_type my_;
771 size_type mz_;
772 const ffd_coefficient *ffd_coefficient_;
773
774 private:
775 const non_rigid_registration_thread& operator =( const non_rigid_registration_thread &p );
776
777 public:
778 void setup_parameters( target_image_type &target, const source_image_type &source, const control_mesh_type &control_mesh,
779 const ffd_coefficient &ffd_coefficient__, size_t mx, size_type my, size_type mz, size_type thread_id, size_type thread_num )
780 {
781 target_ = &target;
782 source_ = &source;
783 control_mesh_ = &control_mesh;
784 ffd_coefficient_ = &ffd_coefficient__;
785 mx_ = mx;
786 my_ = my;
787 mz_ = mz;
788 thread_id_ = thread_id;
789 thread_num_ = thread_num;
790 }
791
792 void setup_parameters( target_image_type &target, const source_image_type &source, const control_mesh_type &control_mesh,
793 size_t mx, size_type my, size_type mz, size_type thread_id, size_type thread_num )
794 {
795 target_ = &target;
796 source_ = &source;
797 control_mesh_ = &control_mesh;
798 ffd_coefficient_ = NULL;
799 mx_ = mx;
800 my_ = my;
801 mz_ = mz;
802 thread_id_ = thread_id;
803 thread_num_ = thread_num;
804 }
805
806 non_rigid_registration_thread( size_type id = 0, size_type num = 1 ) : thread_id_( id ), thread_num_( num ),
807 target_( NULL ), source_( NULL ), control_mesh_( NULL ), mx_( -1 ), my_( -1 ), mz_( -1 ), ffd_coefficient_( NULL )
808 {
809 }
810
811 non_rigid_registration_thread( const non_rigid_registration_thread &p ) : base( p ), thread_id_( p.thread_id_ ), thread_num_( p.thread_num_ ),
812 target_( p.target_ ), source_( p.source_ ), control_mesh_( p.control_mesh_ ), mx_( p.mx_ ), my_( p.my_ ), mz_( p.mz_ ), ffd_coefficient_( p.ffd_coefficient_ )
813 {
814 }
815
816 protected:
817 // 継承した先で必ず実装されるスレッド関数
818 virtual thread_exit_type thread_function( )
819 {
820 if( ffd_coefficient_ != NULL )
821 {
822 non_rigid_transformation_( *target_, *source_, *control_mesh_, *ffd_coefficient_, mx_, my_, mz_, thread_id_, thread_num_ );
823 }
824 else
825 {
826 non_rigid_transformation_( *target_, *source_, *control_mesh_, mx_, my_, mz_, thread_id_, thread_num_ );
827 }
828 return( true );
829 }
830 };
831
832
833 // 制御点情報を用いて，ソース画像をターゲット画像に変形する
834 template < class TARGETTYPE, class SOURCETYPE, class CONTROLMESH >
835 static bool transformation( TARGETTYPE &target, const SOURCETYPE &source, const CONTROLMESH &control_mesh,
836 typename CONTROLMESH::difference_type mx = -1,
837 typename CONTROLMESH::difference_type my = -1,
838 typename CONTROLMESH::difference_type mz = -1,
839 typename CONTROLMESH::size_type thread_num = 0
840 )
841 {
842 if( is_same_object( target, source ) || target.empty( ) || source.empty( ) )
843 {
844 return( false );
845 }
846
847 typedef TARGETTYPE target_image_type;
848 typedef SOURCETYPE source_image_type;
849 typedef CONTROLMESH control_mesh_type;
850 typedef typename target_image_type::size_type size_type;
851 typedef non_rigid_registration_thread< target_image_type, source_image_type, control_mesh_type > _non_rigid_registration_thread_;
852
853 if( thread_num == 0 )
854 {
855 thread_num = static_cast< size_type >( get_cpu_num( ) );
856 }
857
858 _non_rigid_registration_thread_ *thread = new _non_rigid_registration_thread_[ thread_num ];
859
860 size_type i;
861 for( i = 0 ; i < thread_num ; i++ )
862 {
863 thread[ i ].setup_parameters( target, source, control_mesh, mx, my, mz, i, thread_num );
864 }
865
866 // スレッドを実行して，終了まで待機する
867 do_threads_( thread, thread_num );
868
869 delete [] thread;
870
871 return( true );
872 }
873
874 // 制御点情報を用いて，ソース画像をターゲット画像に変形する
875 template < class TARGETTYPE, class SOURCETYPE, class CONTROLMESH >
876 static bool transformation( TARGETTYPE &target, const SOURCETYPE &source, const CONTROLMESH &control_mesh, const ffd_coefficient &ffd_coeff,
877 typename CONTROLMESH::difference_type mx = -1,
878 typename CONTROLMESH::difference_type my = -1,
879 typename CONTROLMESH::difference_type mz = -1,
880 typename CONTROLMESH::size_type thread_num = 0
881 )
882 {
883 if( is_same_object( target, source ) || target.empty( ) || source.empty( ) )
884 {
885 return( false );
886 }
887
888 typedef TARGETTYPE target_image_type;
889 typedef SOURCETYPE source_image_type;
890 typedef CONTROLMESH control_mesh_type;
891 typedef typename target_image_type::size_type size_type;
892 typedef non_rigid_registration_thread< target_image_type, source_image_type, control_mesh_type > _non_rigid_registration_thread_;
893
894 if( thread_num == 0 )
895 {
896 thread_num = static_cast< size_type >( get_cpu_num( ) );
897 }
898
899 _non_rigid_registration_thread_ *thread = new _non_rigid_registration_thread_[ thread_num ];
900
901 size_type i;
902 for( i = 0 ; i < thread_num ; i++ )
903 {
904 thread[ i ].setup_parameters( target, source, control_mesh, ffd_coeff, mx, my, mz, i, thread_num );
905 }
906
907 // スレッドを実行して，終了まで待機する
908 do_threads_( thread, thread_num );
909
910 delete [] thread;
911
912 return( true );
913 }
914
915
916
917 template < class TARGETTYPE, class SOURCETYPE, class CONTROLMESH >
918 struct registration_functor : public thread< registration_functor< TARGETTYPE, SOURCETYPE, CONTROLMESH > >
919 {
920 typedef TARGETTYPE target_image_type;
921 typedef SOURCETYPE source_image_type;
922 typedef CONTROLMESH control_mesh_type;
923 typedef typename TARGETTYPE::size_type size_type;
924 typedef typename TARGETTYPE::difference_type difference_type;
925 typedef typename CONTROLMESH::value_type vector_type;
926 typedef matrix< double > matrix_type;
927 typedef thread< registration_functor< TARGETTYPE, SOURCETYPE, CONTROLMESH > > base;
928 typedef typename base::thread_exit_type thread_exit_type;
929
930 array< unsigned int * > target;
931 target_image_type transformed_image;
932 const source_image_type &source;
933 const control_mesh_type &control_mesh;
934 control_mesh_type control_mesh_tmp;
935 difference_type x, y, z;
936
937 array2< unsigned int > h;
938 array< unsigned int > hh;
939 array< unsigned int > __no_data_is_associated__;
940
941 difference_type minimum;
942 difference_type maximum;
943 size_type BIN;
944 ffd_coefficient ffd_coeff;
945
946 double H1;
947
948 template < class IMAGETYPE >
949 static difference_type get_minimum( const IMAGETYPE &image )
950 {
951 typedef typename IMAGETYPE::value_type value_type;
952 value_type min = image[ 0 ];
953 for( size_type i = 1 ; i < image.size( ) ; i++ )
954 {
955 if( min > image[ i ] )
956 {
957 min = image[ i ];
958 }
959 }
960 return( static_cast< difference_type >( min ) );
961 }
962
963 template < class T >
964 static const T get_minimum( const T &v1, const T &v2 )
965 {
966 return( v1 < v2 ? v1 : v2 );
967 }
968
969 template < class IMAGETYPE >
970 static difference_type get_maximum( const IMAGETYPE &image )
971 {
972 typedef typename IMAGETYPE::value_type value_type;
973 value_type max = image[ 0 ];
974 for( size_type i = 1 ; i < image.size( ) ; i++ )
975 {
976 if( max < image[ i ] )
977 {
978 max = image[ i ];
979 }
980 }
981 return( static_cast< difference_type >( max ) );
982 }
983
984 template < class T >
985 static const T get_maximum( const T &v1, const T &v2 )
986 {
987 return( v1 > v2 ? v1 : v2 );
988 }
989
990 registration_functor( const target_image_type &tgt, const source_image_type &src, const control_mesh_type &cmesh, size_type bin )
991 : target( tgt.size( ) ), transformed_image( tgt ), source( src ), control_mesh ( cmesh ), control_mesh_tmp( cmesh ), x( 0 ), y( 0 ), z( 0 ), BIN( bin )
992 {
993 // FFD 用の係数テーブルを作成する
994 ffd_coeff.initialize( tgt.width( ), tgt.height( ), tgt.depth( ), control_mesh );
995
996 // 初期画像を作成する
997 transformation( transformed_image, source, control_mesh, ffd_coeff );
998
999 minimum = get_maximum( get_minimum( get_minimum( tgt ), get_minimum( src ) ), static_cast< difference_type >(-1000) );
1000 maximum = get_maximum( get_maximum( tgt ), get_maximum( src ) );
1001 h.resize( ( maximum - minimum + 1 ) / BIN + 1, ( maximum - minimum + 1 ) / BIN + 1 );
1002 hh.resize( ( maximum - minimum + 1 ) / BIN + 1, 1 );
1003 __no_data_is_associated__.resize( ( maximum - minimum + 1 ) / BIN + 1, 1 );
1004
1005 size_type count = 0;
1006 difference_type upper = h.width( );
1007 double _1_bin = 1.0 / BIN;
1008 for( size_t i = 0 ; i < tgt.size( ) ; i++ )
1009 {
1010 difference_type v1 = static_cast< difference_type >( ( tgt[ i ] - minimum ) * _1_bin );
1011 if( v1 < 0 || upper <= v1 )
1012 {
1013 target[ i ] = &__no_data_is_associated__[ 0 ];
1014 continue;
1015 }
1016 target[ i ] = &h( 0, v1 );
1017 hh[ v1 ]++;
1018 count++;
1019 }
1020
1021 H1 = 0.0;
1022
1023 // ヒストグラムの正規化
1024 if( count == 0.0 )
1025 {
1026 count = 1;
1027 }
1028
1029 for( size_type i = 0 ; i < hh.size( ) ; i++ )
1030 {
1031 if( hh[ i ] > 0 )
1032 {
1033 double v = hh[ i ] / static_cast< double >( count );
1034 H1 -= v * std::log( v );
1035 }
1036 }
1037 }
1038
1039 void force_initialize( )
1040 {
1041 ffd_coeff.initialize( transformed_image.width( ), transformed_image.height( ), transformed_image.depth( ), control_mesh );
1042 }
1043
1044 void initialize( difference_type i, difference_type j, difference_type k )
1045 {
1046 x = i;
1047 y = j;
1048 z = k;
1049 control_mesh_tmp = control_mesh;
1050 transformation( transformed_image, source, control_mesh, ffd_coeff );
1051 }
1052
1053 void apply_control_point_to_mesh( control_mesh_type &cmesh ) const
1054 {
1055 cmesh( x, y, z ) = control_mesh_tmp( x, y, z );
1056 }
1057
1058 template < class PARAMETER >
1059 double operator ()( const PARAMETER &p )
1060 {
1061 return( evaluate_error( p, false ) );
1062 }
1063
1064 template < class PARAMETER >
1065 double evaluate_error( const PARAMETER &p, bool use_thread = true )
1066 {
1067 control_mesh_tmp( x, y, z ).x = control_mesh( x, y, z ).x + p[ 0 ];
1068 control_mesh_tmp( x, y, z ).y = control_mesh( x, y, z ).y + p[ 1 ];
1069 control_mesh_tmp( x, y, z ).z = control_mesh( x, y, z ).z + p[ 2 ];
1070
1071 if( !use_thread )
1072 {
1073 non_rigid_transformation_( transformed_image, source, control_mesh_tmp, ffd_coeff, x, y, z );
1074 }
1075 else
1076 {
1077 transformation( transformed_image, source, control_mesh_tmp, ffd_coeff, x, y, z );
1078 }
1079
1080 // データを初期化する
1081 __no_data_is_associated__.fill( );
1082 h.fill( );
1083 hh.fill( );
1084 size_type count = 0;
1085 double _1_bin = 1.0 / BIN;
1086
1087 for( size_t i = 0 ; i < target.size( ) ; i++ )
1088 {
1089 difference_type v2 = static_cast< difference_type >( ( transformed_image[ i ] - minimum ) * _1_bin );
1090 if( v2 < 0 )
1091 {
1092 continue;
1093 }
1094
1095 target[ i ][ v2 ]++;
1096 hh[ v2 ]++;
1097 count++;
1098 }
1099
1100 // ヒストグラムの正規化
1101 if( count == 0 )
1102 {
1103 return ( 5.0 );
1104 }
1105
1106 double H2 = 0.0, H12 = 0.0, _1_count = 1.0 / static_cast< double >( count );
1107
1108 for( size_type i = 0 ; i < hh.size( ) ; i++ )
1109 {
1110 if( hh[ i ] > 0 )
1111 {
1112 double v = hh[ i ] * _1_count;
1113 H2 -= v * std::log( v );
1114
1115 array2< unsigned int >::const_pointer ph = &h( 0, i );
1116 for( size_type j = 0 ; j < hh.size( ) ; j++ )
1117 {
1118 if( ph[ j ] > 0 )
1119 {
1120 double v = ph[ j ] * _1_count;
1121 H12 -= v * std::log( v );
1122 }
1123 }
1124 }
1125 }
1126
1127 return( -( H1 + H2 ) / H12 );
1128 //return( H1 + H2 - H12 );
1129 }
1130
1131 protected:
1132 // 最小のバウンディングボックス用の距離を計算する
1133 double min_range( const double &v1, const double &v2, const double &v3 ) const
1134 {
1135 double v = v1 < v2 ? v1 : v2;
1136 v = v < v3 ? v : v3;
1137 return( v );
1138 }
1139
1140 double max_range( const double &v1, const double &v2, const double &v3 ) const
1141 {
1142 double v = v1 > v2 ? v1 : v2;
1143 v = v > v3 ? v : v3;
1144 return( v );
1145 }
1146
1147 double min_range( const double &v1, const double &v2, const double &v3,
1148 const double &v4, const double &v5, const double &v6,
1149 const double &v7, const double &v8, const double &v9 ) const
1150 {
1151 double vv1 = min_range( v1, v2, v3 );
1152 double vv2 = min_range( v4, v5, v6 );
1153 double vv3 = min_range( v7, v8, v9 );
1154 return( min_range( vv1, vv2, vv3 ) );
1155 }
1156
1157 double max_range( const double &v1, const double &v2, const double &v3,
1158 const double &v4, const double &v5, const double &v6,
1159 const double &v7, const double &v8, const double &v9 ) const
1160 {
1161 double vv1 = max_range( v1, v2, v3 );
1162 double vv2 = max_range( v4, v5, v6 );
1163 double vv3 = max_range( v7, v8, v9 );
1164 return( max_range( vv1, vv2, vv3 ) );
1165 }
1166
1167 // 継承した先で必ず実装されるスレッド関数
1168 virtual thread_exit_type thread_function( )
1169 {
1170 typedef __minimization_utility__::__no_copy_constructor_functor__< registration_functor< TARGETTYPE, SOURCETYPE, CONTROLMESH > > no_constructor_functor_type;
1171
1172 // 最小化を開始
1173 double search_length = 0.95;
1174 matrix_type p( 3, 1 ), dirs = matrix_type::identity( 3, 3 ), bound( 3, 2 );
1175
1176 //const vector_type &p01 = control_mesh( x - 1, y - 1, z - 1 );
1177 //const vector_type &p02 = control_mesh( x , y - 1, z - 1 );
1178 //const vector_type &p03 = control_mesh( x + 1, y - 1, z - 1 );
1179 //const vector_type &p04 = control_mesh( x - 1, y , z - 1 );
1180 //const vector_type &p05 = control_mesh( x , y , z - 1 );
1181 //const vector_type &p06 = control_mesh( x + 1, y , z - 1 );
1182 //const vector_type &p07 = control_mesh( x - 1, y + 1, z - 1 );
1183 //const vector_type &p08 = control_mesh( x , y + 1, z - 1 );
1184 //const vector_type &p09 = control_mesh( x + 1, y + 1, z - 1 );
1185 //const vector_type &p11 = control_mesh( x - 1, y - 1, z );
1186 //const vector_type &p12 = control_mesh( x , y - 1, z );
1187 //const vector_type &p13 = control_mesh( x + 1, y - 1, z );
1188 //const vector_type &p14 = control_mesh( x - 1, y , z );
1189 //const vector_type &p15 = control_mesh( x , y , z );
1190 //const vector_type &p16 = control_mesh( x + 1, y , z );
1191 //const vector_type &p17 = control_mesh( x - 1, y + 1, z );
1192 //const vector_type &p18 = control_mesh( x , y + 1, z );
1193 //const vector_type &p19 = control_mesh( x + 1, y + 1, z );
1194 //const vector_type &p21 = control_mesh( x - 1, y - 1, z + 1 );
1195 //const vector_type &p22 = control_mesh( x , y - 1, z + 1 );
1196 //const vector_type &p23 = control_mesh( x + 1, y - 1, z + 1 );
1197 //const vector_type &p24 = control_mesh( x - 1, y , z + 1 );
1198 //const vector_type &p25 = control_mesh( x , y , z + 1 );
1199 //const vector_type &p26 = control_mesh( x + 1, y , z + 1 );
1200 //const vector_type &p27 = control_mesh( x - 1, y + 1, z + 1 );
1201 //const vector_type &p28 = control_mesh( x , y + 1, z + 1 );
1202 //const vector_type &p29 = control_mesh( x + 1, y + 1, z + 1 );
1203
1204 //double lower_x = max_range( p01.x, p04.x, p07.x, p11.x, p14.x, p17.x, p21.x, p24.x, p27.x );
1205 //double lower_y = max_range( p01.y, p02.y, p03.y, p11.y, p12.y, p13.y, p21.y, p22.y, p23.y );
1206 //double lower_z = max_range( p01.z, p02.z, p03.z, p04.z, p05.z, p06.z, p07.z, p08.z, p09.z );
1207 //double upper_x = min_range( p03.x, p06.x, p09.x, p13.x, p16.x, p19.x, p23.x, p26.x, p29.x );
1208 //double upper_y = min_range( p07.y, p08.y, p09.y, p17.y, p18.y, p19.y, p27.y, p28.y, p29.y );
1209 //double upper_z = min_range( p21.z, p22.z, p23.z, p24.z, p25.z, p26.z, p27.z, p28.z, p29.z );
1210
1211
1212 //bound( 0, 0 ) = ( lower_x - p15.x ) * search_length;
1213 //bound( 0, 1 ) = ( upper_x - p15.x ) * search_length;
1214 //bound( 1, 0 ) = ( lower_y - p15.y ) * search_length;
1215 //bound( 1, 1 ) = ( upper_y - p15.y ) * search_length;
1216 //bound( 2, 0 ) = ( lower_z - p15.z ) * search_length;
1217 //bound( 2, 1 ) = ( upper_z - p15.z ) * search_length;
1218 bound( 0, 0 ) = ( control_mesh( x - 1, y , z ).x - control_mesh( x, y, z ).x ) * search_length;
1219 bound( 0, 1 ) = ( control_mesh( x + 1, y , z ).x - control_mesh( x, y, z ).x ) * search_length;
1220 bound( 1, 0 ) = ( control_mesh( x , y - 1, z ).y - control_mesh( x, y, z ).y ) * search_length;
1221 bound( 1, 1 ) = ( control_mesh( x , y + 1, z ).y - control_mesh( x, y, z ).y ) * search_length;
1222 bound( 2, 0 ) = ( control_mesh( x , y , z - 1 ).z - control_mesh( x, y, z ).z ) * search_length;
1223 bound( 2, 1 ) = ( control_mesh( x , y , z + 1 ).z - control_mesh( x, y, z ).z ) * search_length;
1224 //bound( 0, 0 ) = -( control_mesh( x - 1, y , z ) - control_mesh( x, y, z ) ).length( ) * search_length;
1225 //bound( 0, 1 ) = ( control_mesh( x + 1, y , z ) - control_mesh( x, y, z ) ).length( ) * search_length;
1226 //bound( 1, 0 ) = -( control_mesh( x , y - 1, z ) - control_mesh( x, y, z ) ).length( ) * search_length;
1227 //bound( 1, 1 ) = ( control_mesh( x , y + 1, z ) - control_mesh( x, y, z ) ).length( ) * search_length;
1228 //bound( 2, 0 ) = -( control_mesh( x , y , z - 1 ) - control_mesh( x, y, z ) ).length( ) * search_length;
1229 //bound( 2, 1 ) = ( control_mesh( x , y , z + 1 ) - control_mesh( x, y, z ) ).length( ) * search_length;
1230
1231 matrix_type dir( p.size( ), 1 ), tmp( p.size( ), 1 );
1232 double v1, v2;
1233 size_type i;
1234
1235 // 他変数関数を１変数関数に変換する
1236 no_constructor_functor_type f( *this );
1237 __minimization_utility__::__convert_to_vector_functor__< typename matrix_type::value_type, typename matrix_type::allocator_type, no_constructor_functor_type > functor( p, dir, tmp, f );
1238
1239 // 勾配方向を計算する
1240 double len = 0.0, distance = 1.0;
1241 for( i = 0 ; i < dir.size( ) ; i++ )
1242 {
1243 tmp[ i ] = p[ i ] + distance;
1244 v1 = f( tmp );
1245
1246 tmp[ i ] = p[ i ] - distance;
1247 v2 = f( tmp );
1248
1249 tmp[ i ] = p[ i ];
1250
1251 dir[ i ] = v2 - v1;
1252 len += dir[ i ] * dir[ i ];
1253 }
1254
1255 if( len > 0 )
1256 {
1257 // 勾配方向ベクトルの正規化
1258 len = std::sqrt( len );
1259 for( i = 0 ; i < dir.size( ) ; i++ )
1260 {
1261 dir[ i ] /= len;
1262 }
1263 }
1264 else
1265 {
1266 // 勾配の計算ができなくなったので終了する
1267 return( false );
1268 }
1269
1270 double l1, l2, xx = 0.0;
1271 if( __minimization_utility__::clipping_length( l1, l2, p, dir, bound ) )
1272 {
1273 // Brent の2次収束アルゴリズムを用いて dir 方向への最小化を行う
1274 mist::brent::minimization( l1, l2, xx, functor, 0.1, 100, false );
1275 p = dir * xx;
1276 }
1277 else
1278 {
1279 return( false );
1280 }
1281
1282 // gradient::minimization( p, bound, no_constructor_functor_type( *this ), 1.0 );
1283 // gradient::minimization( p, bound, no_constructor_functor_type( *this ), 0.1, 1.0, 1 );
1284
1285 // 結果を反映
1286 vector_type &v = control_mesh_tmp( x, y, z );
1287 v = control_mesh( x, y, z );
1288 v.x += p[ 0 ];
1289 v.y += p[ 1 ];
1290 v.z += p[ 2 ];
1291
1292 return( true );
1293 }
1294 };
1295 }
1296
1297
1298
1306
1307
1308
1310 namespace non_rigid
1311 {
1324 template < class TARGETTYPE >
1325 class registration
1326 {
1327 public:
1328 typedef TARGETTYPE image_type;
1329 typedef vector3< double > vector_type;
1330 typedef marray< array3< vector_type > > control_mesh_type;
1331 typedef matrix< double > matrix_type;
1332 typedef typename image_type::size_type size_type;
1333 typedef typename image_type::difference_type difference_type;
1334
1335 protected:
1336 const image_type &target;
1337 control_mesh_type control_mesh;
1338 size_type BIN;
1339
1340 public:
1349 registration( const image_type &tgt, size_type divW, size_type divH, size_type divD, size_type bin )
1350 : target( tgt ), control_mesh( divW + 1, divH + 1, divD + 1, 2 ), BIN( bin )
1351 {
1352 }
1353
1354
1363 template < class SOURCETYPE >
1364 void apply( const SOURCETYPE &source, double tolerance, size_type max_loop = 3, size_type coarse_to_fine_step = 1, size_type thread_num = 0 )
1365 {
1366 apply( source, tolerance, max_loop, coarse_to_fine_step, thread_num, __mist_dmy_callback__( ) );
1367 }
1368
1378 template < class SOURCETYPE, class Functor >
1379 void apply( const SOURCETYPE &source, double tolerance, size_type max_loop, size_type coarse_to_fine_step, size_type thread_num, Functor callback )
1380 {
1381 typedef __non_rigid_registration_utility__::registration_functor< TARGETTYPE, SOURCETYPE, control_mesh_type > non_rigid_registration_functor_type;
1382 typedef __minimization_utility__::__no_copy_constructor_functor__< non_rigid_registration_functor_type > no_constructor_functor_type;
1383
1384 size_type divW = control_mesh.width( ) - 1;
1385 size_type divH = control_mesh.height( ) - 1;
1386 size_type divD = control_mesh.depth( ) - 1;
1387
1388 // まず，制御点の初期化
1389 double stepW = divW == 0 ? 0 : source.width( ) / static_cast< double >( divW );
1390 double stepH = divH == 0 ? 0 : source.height( ) / static_cast< double >( divH );
1391 double stepD = divD == 0 ? 0 : source.depth( ) / static_cast< double >( divD );
1392 difference_type w = control_mesh.width( );
1393 difference_type h = control_mesh.height( );
1394 difference_type d = control_mesh.depth( );
1395 for( difference_type k = -2 ; k <= d + 1 ; k++ )
1396 {
1397 for( difference_type j = -2 ; j <= h + 1 ; j++ )
1398 {
1399 for( difference_type i = -2 ; i <= w + 1 ; i++ )
1400 {
1401 vector_type &v = control_mesh( i, j, k );
1402 v.x = i * stepW * source.reso1( );
1403 v.y = j * stepH * source.reso2( );
1404 v.z = k * stepD * source.reso3( );
1405 }
1406 }
1407 }
1408
1409 callback( 0.0 );
1410
1411 if( thread_num == 0 )
1412 {
1413 thread_num = static_cast< size_type >( get_cpu_num( ) );
1414 }
1415
1416 // ノンリジッドレジストレーションを行うファンクタ
1417 non_rigid_registration_functor_type **f = new non_rigid_registration_functor_type*[ thread_num ];
1418
1419 for( size_type i = 0 ; i < thread_num ; i++ )
1420 {
1421 f[ i ] = new non_rigid_registration_functor_type( target, source, control_mesh, BIN );
1422 }
1423
1424 // 一番初期時点での評価値を計算しておく
1425 double err = f[ 0 ]->evaluate_error( matrix_type::zero( 3, 1 ) );
1426 double old_err = err;
1427
1428 typedef mist::vector3< size_type > control_point_index_type;
1429 std::vector< control_point_index_type > control_point_list;
1430
1431 //std::cout << "制御点数1: " << control_point_list.size( ) << std::endl;
1432 //std::cout << "制御点数2: " << w * h * d << std::endl;
1433
1434
1435 size_type fine_step_loop = 0;
1436 while( fine_step_loop++ < coarse_to_fine_step )
1437 {
1438 {
1439 // 1階の探索で，互いに影響しあわない制御点のリストを作成する
1440 difference_type w = control_mesh.width( );
1441 difference_type h = control_mesh.height( );
1442 difference_type d = control_mesh.depth( );
1443 size_type X[] = { 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2 };
1444 size_type Y[] = { 0, 0, 0, 1, 1, 1, 2, 2, 2, 0, 0, 0, 1, 1, 1, 2, 2, 2, 0, 0, 0, 1, 1, 1, 2, 2, 2 };
1445 size_type Z[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2 };
1446
1447 control_point_list.clear( );
1448
1449 for( size_type s = 0 ; s < 27 ; s++ )
1450 {
1451 for( difference_type k = Z[ s ] ; k < d ; k += 3 )
1452 {
1453 for( difference_type j = Y[ s ] ; j < h ; j += 3 )
1454 {
1455 for( difference_type i = X[ s ] ; i < w ; i += 3 )
1456 {
1457 control_point_list.push_back( mist::vector3< size_type >( i, j, k ) );
1458 }
1459 }
1460 }
1461 }
1462 }
1463
1464 __mist_convert_callback__< Functor > callback_( callback, ( fine_step_loop - 1 ) * 100.0 / coarse_to_fine_step, fine_step_loop * 100.0 / coarse_to_fine_step );
1465
1466 size_type loop = 0, count = 0;
1467 double max_iteration_num = control_point_list.size( ) * max_loop / static_cast< double >( thread_num );
1468 while( loop++ < max_loop )
1469 {
1470 for( size_type i = 0 ; i < control_point_list.size( ) ; )
1471 {
1472 size_type numthreads;
1473 for( numthreads = 0 ; numthreads < thread_num && i < control_point_list.size( ) ; numthreads++, i++ )
1474 {
1475 control_point_index_type &v = control_point_list[ i ];
1476 f[ numthreads ]->initialize( v.x, v.y, v.z );
1477 }
1478
1479 // スレッドの生成
1480 for( size_type t = 0 ; t < numthreads ; t++ )
1481 {
1482 f[ t ]->create( );
1483 }
1484
1485 // スレッドの終了待ち
1486 for( size_type t = 0 ; t < numthreads ; t++ )
1487 {
1488 f[ t ]->wait( INFINITE );
1489 }
1490
1491 // リソースの開放
1492 for( size_type t = 0 ; t < numthreads ; t++ )
1493 {
1494 f[ t ]->close( );
1495 }
1496
1497 // 探索の結果を，各制御点に反映する
1498 for( size_type t = 0 ; t < numthreads ; t++ )
1499 {
1500 f[ t ]->apply_control_point_to_mesh( control_mesh );
1501 }
1502
1503 callback_( 100.0 / max_iteration_num * count++ );
1504 }
1505
1506 err = f[ 0 ]->evaluate_error( matrix_type::zero( 3, 1 ) );
1507
1508 if( 2.0 * std::abs( old_err - err ) < tolerance * ( std::abs( old_err ) + std::abs( err ) ) )
1509 {
1510 break;
1511 }
1512
1513 old_err = err;
1514 }
1515
1516
1517 if( fine_step_loop < coarse_to_fine_step )
1518 {
1519 // メッシュの再分割を行う
1520 control_mesh_type cmesh = control_mesh;
1521 difference_type w = control_mesh.width( ) * 2 - 1;
1522 difference_type h = control_mesh.height( ) * 2 - 1;
1523 difference_type d = control_mesh.depth( ) * 2 - 1;
1524 w = w <= 0 ? 1 : w;
1525 h = h <= 0 ? 1 : h;
1526 d = d <= 0 ? 1 : d;
1527 control_mesh.resize( w, h, d );
1528
1529 double stepW = w == 1 ? 0 : source.width( ) / static_cast< double >( w - 1 );
1530 double stepH = h == 1 ? 0 : source.height( ) / static_cast< double >( h - 1 );
1531 double stepD = d == 1 ? 0 : source.depth( ) / static_cast< double >( d - 1 );
1532 for( difference_type k = -2 ; k <= d + 1 ; k++ )
1533 {
1534 for( difference_type j = -2 ; j <= h + 1 ; j++ )
1535 {
1536 for( difference_type i = -2 ; i <= w + 1 ; i++ )
1537 {
1538 vector_type &v = control_mesh( i, j, k );
1539 v.x = i * stepW * source.reso1( );
1540 v.y = j * stepH * source.reso2( );
1541 v.z = k * stepD * source.reso3( );
1542 }
1543 }
1544 }
1545
1546 for( difference_type k = 0 ; k < d ; k++ )
1547 {
1548 bool k_is_odd = ( k % 2 ) == 1;
1549 difference_type kk = k / 2;
1550
1551 for( difference_type j = 0 ; j < h ; j++ )
1552 {
1553 bool j_is_odd = ( j % 2 ) == 1;
1554 difference_type jj = j / 2;
1555
1556 for( difference_type i = 0 ; i < w ; i++ )
1557 {
1558 bool i_is_odd = ( i % 2 ) == 1;
1559 difference_type ii = i / 2;
1560
1561 if( i_is_odd && j_is_odd && k_is_odd )
1562 {
1563 control_mesh( i, j, k ) = ( cmesh( ii , jj , kk ) + cmesh( ii , jj , kk + 1 ) +
1564 cmesh( ii , jj + 1, kk ) + cmesh( ii , jj + 1, kk + 1 ) +
1565 cmesh( ii + 1, jj , kk ) + cmesh( ii + 1, jj , kk + 1 ) +
1566 cmesh( ii + 1, jj + 1, kk ) + cmesh( ii + 1, jj + 1, kk + 1 ) ) / 8.0;
1567 }
1568 else if( !i_is_odd && j_is_odd && k_is_odd )
1569 {
1570 control_mesh( i, j, k ) = (
1571 cmesh( ii - 1, jj , kk ) + cmesh( ii - 1, jj , kk + 1 ) +
1572 cmesh( ii - 1, jj + 1, kk ) + cmesh( ii - 1, jj + 1, kk + 1 ) +
1573 cmesh( ii + 1, jj , kk ) + cmesh( ii + 1, jj , kk + 1 ) +
1574 cmesh( ii + 1, jj + 1, kk ) + cmesh( ii + 1, jj + 1, kk + 1 ) +
1575 6.0 * (
1576 cmesh( ii , jj , kk ) + cmesh( ii , jj , kk + 1 ) +
1577 cmesh( ii , jj + 1, kk ) + cmesh( ii , jj + 1, kk + 1 )
1578 )
1579 ) / 32.0;
1580 }
1581 else if( i_is_odd && !j_is_odd && k_is_odd )
1582 {
1583 control_mesh( i, j, k ) = (
1584 cmesh( ii , jj - 1, kk ) + cmesh( ii , jj - 1, kk + 1 ) +
1585 cmesh( ii + 1, jj - 1, kk ) + cmesh( ii + 1, jj - 1, kk + 1 ) +
1586 cmesh( ii , jj + 1, kk ) + cmesh( ii , jj + 1, kk + 1 ) +
1587 cmesh( ii + 1, jj + 1, kk ) + cmesh( ii + 1, jj + 1, kk + 1 ) +
1588 6.0 * (
1589 cmesh( ii , jj , kk ) + cmesh( ii , jj , kk + 1 ) +
1590 cmesh( ii + 1, jj , kk ) + cmesh( ii + 1, jj , kk + 1 )
1591 )
1592 ) / 32.0;
1593 }
1594 else if( i_is_odd && j_is_odd && !k_is_odd )
1595 {
1596 control_mesh( i, j, k ) = (
1597 cmesh( ii , jj , kk - 1 ) + cmesh( ii + 1, jj , kk - 1 ) +
1598 cmesh( ii , jj + 1, kk - 1 ) + cmesh( ii + 1, jj + 1, kk - 1 ) +
1599 cmesh( ii , jj , kk + 1 ) + cmesh( ii + 1, jj , kk + 1 ) +
1600 cmesh( ii , jj + 1, kk + 1 ) + cmesh( ii + 1, jj + 1, kk + 1 ) +
1601 6.0 * (
1602 cmesh( ii , jj , kk ) + cmesh( ii + 1, jj , kk ) +
1603 cmesh( ii , jj + 1, kk ) + cmesh( ii + 1, jj + 1, kk )
1604 )
1605 ) / 32.0;
1606 }
1607 else if( !i_is_odd && !j_is_odd && k_is_odd )
1608 {
1609 control_mesh( i, j, k ) = (
1610 cmesh( ii - 1, jj - 1, kk ) + cmesh( ii - 1, jj + 1, kk ) +
1611 cmesh( ii + 1, jj - 1, kk ) + cmesh( ii + 1, jj + 1, kk ) +
1612 cmesh( ii - 1, jj - 1, kk + 1 ) + cmesh( ii - 1, jj + 1, kk + 1 ) +
1613 cmesh( ii + 1, jj - 1, kk + 1 ) + cmesh( ii + 1, jj + 1, kk + 1 ) +
1614 6.0 * (
1615 cmesh( ii - 1, jj , kk ) + cmesh( ii , jj - 1, kk ) +
1616 cmesh( ii + 1, jj , kk ) + cmesh( ii , jj + 1, kk ) +
1617 cmesh( ii - 1, jj , kk + 1 ) + cmesh( ii , jj - 1, kk + 1 ) +
1618 cmesh( ii + 1, jj , kk + 1 ) + cmesh( ii , jj + 1, kk + 1 )
1619 ) +
1620 36.0 * (
1621 cmesh( ii , jj , kk ) + cmesh( ii , jj , kk + 1 )
1622 )
1623 ) / 128.0;
1624 }
1625 else if( !i_is_odd && j_is_odd && !k_is_odd )
1626 {
1627 control_mesh( i, j, k ) = (
1628 cmesh( ii - 1, jj , kk - 1 ) + cmesh( ii - 1, jj , kk + 1 ) +
1629 cmesh( ii + 1, jj , kk - 1 ) + cmesh( ii + 1, jj , kk + 1 ) +
1630 cmesh( ii - 1, jj + 1, kk - 1 ) + cmesh( ii - 1, jj + 1, kk + 1 ) +
1631 cmesh( ii + 1, jj + 1, kk - 1 ) + cmesh( ii + 1, jj + 1, kk + 1 ) +
1632 6.0 * (
1633 cmesh( ii - 1, jj , kk ) + cmesh( ii , jj , kk - 1 ) +
1634 cmesh( ii + 1, jj , kk ) + cmesh( ii , jj , kk + 1 ) +
1635 cmesh( ii - 1, jj + 1, kk ) + cmesh( ii , jj + 1, kk - 1 ) +
1636 cmesh( ii + 1, jj + 1, kk ) + cmesh( ii , jj + 1, kk + 1 )
1637 ) +
1638 36.0 * (
1639 cmesh( ii , jj , kk ) + cmesh( ii , jj + 1, kk )
1640 )
1641 ) / 128.0;
1642 }
1643 else if( i_is_odd && !j_is_odd && !k_is_odd )
1644 {
1645 control_mesh( i, j, k ) = (
1646 cmesh( ii , jj - 1, kk - 1 ) + cmesh( ii , jj + 1, kk - 1 ) +
1647 cmesh( ii , jj - 1, kk + 1 ) + cmesh( ii , jj + 1, kk + 1 ) +
1648 cmesh( ii + 1, jj - 1, kk - 1 ) + cmesh( ii + 1, jj + 1, kk - 1 ) +
1649 cmesh( ii + 1, jj - 1, kk + 1 ) + cmesh( ii + 1, jj + 1, kk + 1 ) +
1650 6.0 * (
1651 cmesh( ii , jj , kk - 1 ) + cmesh( ii , jj - 1, kk ) +
1652 cmesh( ii , jj , kk + 1 ) + cmesh( ii , jj + 1, kk ) +
1653 cmesh( ii + 1, jj , kk - 1 ) + cmesh( ii + 1, jj - 1, kk ) +
1654 cmesh( ii + 1, jj , kk + 1 ) + cmesh( ii + 1, jj + 1, kk )
1655 ) +
1656 36.0 * (
1657 cmesh( ii , jj , kk ) + cmesh( ii + 1, jj , kk )
1658 )
1659 ) / 128.0;
1660 }
1661 else
1662 {
1663 control_mesh( i, j, k ) = (
1664 cmesh( ii - 1, jj - 1, kk - 1 ) + cmesh( ii - 1, jj + 1, kk - 1 ) +
1665 cmesh( ii + 1, jj - 1, kk - 1 ) + cmesh( ii + 1, jj + 1, kk - 1 ) +
1666 cmesh( ii - 1, jj - 1, kk + 1 ) + cmesh( ii - 1, jj + 1, kk + 1 ) +
1667 cmesh( ii + 1, jj - 1, kk + 1 ) + cmesh( ii + 1, jj + 1, kk + 1 ) +
1668 6.0 * (
1669 cmesh( ii - 1, jj - 1, kk ) + cmesh( ii - 1, jj + 1, kk ) +
1670 cmesh( ii + 1, jj - 1, kk ) + cmesh( ii + 1, jj + 1, kk ) +
1671 cmesh( ii - 1, jj , kk - 1 ) + cmesh( ii , jj - 1, kk - 1 ) +
1672 cmesh( ii + 1, jj , kk - 1 ) + cmesh( ii , jj + 1, kk - 1 ) +
1673 cmesh( ii - 1, jj , kk + 1 ) + cmesh( ii , jj - 1, kk + 1 ) +
1674 cmesh( ii + 1, jj , kk + 1 ) + cmesh( ii , jj + 1, kk + 1 )
1675 ) +
1676 36.0 * (
1677 cmesh( ii , jj , kk - 1 ) + cmesh( ii , jj , kk + 1 ) +
1678 cmesh( ii - 1, jj , kk ) + cmesh( ii , jj - 1, kk ) +
1679 cmesh( ii + 1, jj , kk ) + cmesh( ii , jj + 1, kk )
1680 ) +
1681 cmesh( ii , jj , kk ) * 216.0
1682 ) / 512.0;
1683 }
1684 }
1685 }
1686 }
1687
1688 // メッシュの再分割を行ったので，FFDの係数の再設定を行う
1689 for( size_type t = 0 ; t < thread_num ; t++ )
1690 {
1691 f[ t ]->force_initialize( );
1692 }
1693 }
1694 }
1695
1696 for( size_type i = 0 ; i < thread_num ; i++ )
1697 {
1698 delete f[ i ];
1699 }
1700 delete [] f;
1701
1702 callback( 100.1 );
1703 }
1704
1705
1715 template < class SOURCETYPE, class SOURCETYPEA, class OUTTYPE, class OUTTYPEA >
1716 void transform( const array2< SOURCETYPE, SOURCETYPEA > &source, array2< OUTTYPE, OUTTYPEA > &out )
1717 {
1718 __non_rigid_registration_utility__::transformation( out, source, control_mesh );
1719 }
1720
1721
1731 template < class SOURCETYPE, class SOURCETYPEA, class OUTTYPE, class OUTTYPEA >
1732 void transform( const array3< SOURCETYPE, SOURCETYPEA > &source, array3< OUTTYPE, OUTTYPEA > &out )
1733 {
1734 __non_rigid_registration_utility__::transformation( out, source, control_mesh );
1735 }
1736
1737
1748 template < class SOURCETYPE, class SOURCETYPEA, class OUTTYPE, class OUTTYPEA >
1749 void transform( const array2< SOURCETYPE, SOURCETYPEA > &source, array2< OUTTYPE, OUTTYPEA > &out, typename array2< OUTTYPE, OUTTYPEA >::value_type value )
1750 {
1751 __non_rigid_registration_utility__::transformation( out, source, control_mesh );
1752
1753 difference_type w = control_mesh.width( );
1754 difference_type h = control_mesh.height( );
1755 difference_type d = control_mesh.depth( );
1756 difference_type i, j;
1757 double ax = out.reso1( );
1758 double ay = out.reso2( );
1759
1760 for( j = 0 ; j < h - 1 ; j++ )
1761 {
1762 for( i = 0 ; i < w - 1 ; i++ )
1763 {
1764 vector_type &vec0 = control_mesh( i , j , 0 );
1765 vector_type &vec1 = control_mesh( i + 1, j , 0 );
1766 vector_type &vec2 = control_mesh( i , j + 1, 0 );
1767
1768 size_type x0 = static_cast< size_type >( vec0.x / ax );
1769 size_type y0 = static_cast< size_type >( vec0.y / ay );
1770 size_type x1 = static_cast< size_type >( vec1.x / ax );
1771 size_type y1 = static_cast< size_type >( vec1.y / ay );
1772 size_type x2 = static_cast< size_type >( vec2.x / ax );
1773 size_type y2 = static_cast< size_type >( vec2.y / ay );
1774 draw_line( out, x0, y0, x1, y1, value );
1775 draw_line( out, x0, y0, x2, y2, value );
1776 }
1777
1778 vector_type &vec0 = control_mesh( i , j , 0 );
1779 vector_type &vec2 = control_mesh( i , j + 1, 0 );
1780 size_type x0 = static_cast< size_type >( vec0.x / ax );
1781 size_type y0 = static_cast< size_type >( vec0.y / ay );
1782 size_type x2 = static_cast< size_type >( vec2.x / ax );
1783 size_type y2 = static_cast< size_type >( vec2.y / ay );
1784 draw_line( out, x0, y0, x2, y2, value );
1785 }
1786 for( i = 0 ; i < w - 1 ; i++ )
1787 {
1788 vector_type &vec0 = control_mesh( i , j , 0 );
1789 vector_type &vec1 = control_mesh( i + 1, j , 0 );
1790
1791 size_type x0 = static_cast< size_type >( vec0.x / ax );
1792 size_type y0 = static_cast< size_type >( vec0.y / ay );
1793 size_type x1 = static_cast< size_type >( vec1.x / ax );
1794 size_type y1 = static_cast< size_type >( vec1.y / ay );
1795 draw_line( out, x0, y0, x1, y1, value );
1796 }
1797 }
1798
1799
1810 template < class SOURCETYPE, class SOURCETYPEA, class OUTTYPE, class OUTTYPEA >
1811 void transform( const array3< SOURCETYPE, SOURCETYPEA > &source, array3< OUTTYPE, OUTTYPEA > &out, typename array2< OUTTYPE, OUTTYPEA >::value_type value )
1812 {
1813 __non_rigid_registration_utility__::transformation( out, source, control_mesh );
1814
1815 difference_type w = control_mesh.width( );
1816 difference_type h = control_mesh.height( );
1817 difference_type d = control_mesh.depth( );
1818 double ax = out.reso1( );
1819 double ay = out.reso2( );
1820 double az = out.reso3( );
1821
1822 for( difference_type k = 0 ; k < d ; k++ )
1823 {
1824 for( difference_type j = 0 ; j < h ; j++ )
1825 {
1826 for( difference_type i = 0 ; i < w ; i++ )
1827 {
1828 const vector_type &vec0 = control_mesh( i, j, k );
1829
1830 size_type x0 = static_cast< size_type >( vec0.x / ax );
1831 size_type y0 = static_cast< size_type >( vec0.y / ay );
1832 size_type z0 = static_cast< size_type >( vec0.z / az );
1833
1834 if( i < w - 1 )
1835 {
1836 const vector_type &vec1 = control_mesh( i + 1, j, k );
1837 size_type x1 = static_cast< size_type >( vec1.x / ax );
1838 size_type y1 = static_cast< size_type >( vec1.y / ay );
1839 size_type z1 = static_cast< size_type >( vec1.z / az );
1840 draw_line( out, x0, y0, z0, x1, y1, z1, value );
1841 }
1842
1843 if( j < h - 1 )
1844 {
1845 const vector_type &vec2 = control_mesh( i, j + 1, k );
1846 size_type x2 = static_cast< size_type >( vec2.x / ax );
1847 size_type y2 = static_cast< size_type >( vec2.y / ay );
1848 size_type z2 = static_cast< size_type >( vec2.z / az );
1849 draw_line( out, x0, y0, z0, x2, y2, z2, value );
1850 }
1851
1852 if( k < d - 1 )
1853 {
1854 const vector_type &vec3 = control_mesh( i, j, k + 1 );
1855 size_type x3 = static_cast< size_type >( vec3.x / ax );
1856 size_type y3 = static_cast< size_type >( vec3.y / ay );
1857 size_type z3 = static_cast< size_type >( vec3.z / az );
1858 draw_line( out, x0, y0, z0, x3, y3, z3, value );
1859 }
1860 }
1861 }
1862 }
1863 }
1864 };
1865
1866 } // 名前空間 non_rigid の終わり
1867
1868 namespace _poc_
1869 {
1870 template < class T> class point
1871 {
1872 public :
1873 T value;
1874 int x;
1875 int y;
1876 };
1877
1878 // 象限の入れ換え(1と3、2と4)
1879 // 中心に直流成分が来る
1880 inline mist::array2< double > shuffle_image( const mist::array2< double > &image )
1881 {
1882 mist::array2< double > tmp;
1883 int width =image.width();
1884 int height = image.height();
1885
1886 tmp.resize( width , height );
1887
1888 #pragma omp parallel for schedule( guided )
1889 for( int h = 0 ; h < height ; h++)
1890 {
1891 for( int w = 0 ; w < width ; w++ )
1892 {
1893 if( w < width / 2 && h < height / 2 )
1894 {
1895 tmp( w + width / 2, h + height / 2 ) = image( w, h );
1896 }
1897 else if( w < width / 2 && h >= height / 2 )
1898 {
1899 tmp( w + width / 2, h - height / 2 ) = image( w, h );
1900
1901 }
1902 else if( w >= width / 2 && h < height / 2 )
1903 {
1904 tmp( w - width / 2, h + height / 2 ) = image( w, h );
1905 }
1906 else if( w >= width / 2 && h >= height / 2 )
1907 {
1908 tmp( w - width / 2, h - height / 2 ) = image( w, h );
1909 }
1910 }
1911 }
1912 return tmp;
1913 }
1914
1915 // ローパスフィルタ
1916 inline mist::array2< std::complex< double > > lowpass_filter( const mist::array2< std::complex< double > > &image, const double lowpass_range )
1917 {
1918 mist::array2< std::complex< double > > tmp;
1919 int width =image.width();
1920 int height = image.height();
1921
1922 tmp.resize( width, height );
1923 std::complex< double > zero( 0 , 0 );
1924
1925 int width1 = (int)(width * lowpass_range / 2);
1926 int width2 = (int)(width * (1-lowpass_range / 2));
1927 int height1 = (int)(height * lowpass_range / 2);
1928 int height2 = (int)(height * (1-lowpass_range / 2));
1929 #pragma omp parallel for schedule( guided )
1930 for( int h = 0 ; h < height ; h++)
1931 {
1932 for( int w = 0 ; w < width ; w++ )
1933 {
1934 if( !( ( w < width1 || w > width2 ) && ( h < height1 || h > height2 ) ) )
1935 {
1936 tmp( w, h ) = zero;
1937 }
1938 else{
1939 tmp( w, h ) = image( w, h );
1940 }
1941 }
1942 }
1943 return tmp;
1944 }
1945
1946 //サブピークを探索
1947 inline void serch_peak( const mist::array2< double > &poc_image, point< double > &peak, point< double > &x_peak, point< double > &y_peak )
1948 {
1949 // ピーク位置の探索
1950 peak.value = poc_image( 0, 0 );
1951 peak.x = 0; peak.y = 0;
1952 for( size_t h = 0 ; h < poc_image.height() ; h++)
1953 {
1954 for( size_t w = 0 ; w < poc_image.width() ; w++ )
1955 {
1956 if( peak.value < poc_image( w, h ) )
1957 {
1958 peak.value = poc_image( w, h );
1959 peak.x = w; peak.y = h;
1960 }
1961 }
1962 }
1963
1964 // サブピーク位置の探索
1966 if( peak.x == 0 )
1967 {
1968 x_peak.value = poc_image( peak.x + 1, peak.y );
1969 x_peak.x = peak.x + 1;
1970 x_peak.y = peak.y;
1971 }
1972 else if( peak.x == poc_image.width() )
1973 {
1974 x_peak.value = poc_image( peak.x - 1, peak.y );
1975 x_peak.x = peak.x - 1;
1976 x_peak.y = peak.y;
1977 }
1978 else if( poc_image( peak.x + 1, peak.y ) >= poc_image( peak.x - 1, peak.y ) )
1979 {
1980 x_peak.value = poc_image( peak.x + 1, peak.y );
1981 x_peak.x = peak.x + 1;
1982 x_peak.y = peak.y;
1983 }
1984 else
1985 {
1986 x_peak.value = poc_image( peak.x - 1, peak.y );
1987 x_peak.x = peak.x - 1;
1988 x_peak.y = peak.y;
1989 }
1990
1992 if( peak.y == 0 )
1993 {
1994 y_peak.value = poc_image( peak.x, peak.y + 1 );
1995 y_peak.x = peak.x;
1996 y_peak.y = peak.y + 1;
1997 }
1998 else if( peak.y == poc_image.height() )
1999 {
2000 y_peak.value = poc_image( peak.x, peak.y - 1 );
2001 y_peak.x = peak.x;
2002 y_peak.y = peak.y - 1;
2003 }
2004 else if( poc_image( peak.x, peak.y + 1 ) >= poc_image( peak.x, peak.y - 1 ) )
2005 {
2006 y_peak.value = poc_image( peak.x, peak.y + 1 );
2007 y_peak.x = peak.x;
2008 y_peak.y = peak.y + 1;
2009 }
2010 else
2011 {
2012 y_peak.value = poc_image( peak.x, peak.y - 1 );
2013 y_peak.x = peak.x;
2014 y_peak.y = peak.y - 1;
2015 }
2016 }
2017 //PEF
2018 inline void PEF( const mist::array2< double > &poc_image, point< double > &peak, point< double > &x_peak, point< double > &y_peak, double &delta_x, double &delta_y, const int distance0, const double lowpass_range )
2019 {
2020 double ux, uy, vx, vy;
2021 double UUx = 0.0 , UVx = 0.0, UUy = 0.0, UVy = 0.0;
2022 int distance = distance0;
2023 int width = poc_image.width();
2024 int height = poc_image.height();
2025
2026 // distanceの例外処理(範囲外のアクセス防止) --------------------------------------
2027 if(distance > peak.x) distance = peak.x;
2028 if(distance > peak.y) distance = peak.y;
2029 if(distance > width - peak.x) distance = width - peak.x;
2030 if(distance > height - peak.y) distance = height - peak.y;
2031
2032 for( int d = 1; d <= distance; d++ )
2033 {
2034 //peak x方向
2035 ux = poc_image( peak.x - d, peak.y ) + poc_image( peak.x + d, peak.y ) - 2 * cos( M_PI * d * lowpass_range ) * poc_image( peak.x, peak.y );
2036 vx = 2 * cos( M_PI * d * lowpass_range ) * poc_image( peak.x, peak.y ) * ( peak.x - width / 2 )
2037 - ( peak.x - width / 2 - d ) * poc_image( peak.x - d, peak.y )
2038 - ( peak.x - width / 2 + d ) * poc_image( peak.x + d, peak.y );
2039 UVx += ux*vx;
2040 UUx += ux*ux;
2041
2042 //peak y方向
2043 uy = poc_image( peak.x, peak.y - d ) + poc_image( peak.x, peak.y + d ) - 2 * cos( M_PI * d * lowpass_range ) * poc_image( peak.x, peak.y );
2044 vy = 2 * cos( M_PI * d * lowpass_range ) * poc_image( peak.x, peak.y ) * ( peak.y - height / 2 )
2045 - ( peak.y - height / 2 - d ) * poc_image( peak.x, peak.y - d )
2046 - ( peak.y - height / 2 + d ) * poc_image( peak.x, peak.y + d );
2047 UVy += uy*vy;
2048 UUy += uy*uy;
2049
2050 //peak2 x方向
2051 ux = poc_image( x_peak.x - d, x_peak.y ) + poc_image( x_peak.x + d, x_peak.y ) - 2 * cos( M_PI * d * lowpass_range ) * poc_image( x_peak.x, x_peak.y );
2052 vx = 2 * cos( M_PI * d * lowpass_range ) * poc_image( x_peak.x, x_peak.y ) * ( x_peak.x - width / 2 )
2053 - ( x_peak.x - width / 2 - d ) * poc_image( x_peak.x - d, x_peak.y )
2054 - ( x_peak.x - width / 2 + d ) * poc_image( x_peak.x + d, x_peak.y );
2055 UVx += ux*vx;
2056 UUx += ux*ux;
2057
2058 //peak2 y方向
2059 uy = poc_image( y_peak.x, y_peak.y - d ) + poc_image( y_peak.x, y_peak.y + d ) - 2 * cos( M_PI * d * lowpass_range ) * poc_image( y_peak.x, y_peak.y );
2060 vy = 2 * cos( M_PI * d * lowpass_range ) * poc_image( y_peak.x, y_peak.y ) * ( y_peak.y - height / 2 )
2061 - ( y_peak.y - height / 2 - d ) * poc_image( y_peak.x, y_peak.y - d )
2062 - ( y_peak.y - height / 2 + d ) * poc_image( y_peak.x, y_peak.y + d );
2063 UVy += uy*vy;
2064 UUy += uy*uy;
2065 }
2066
2067 //std::cout << UVx << "," << UUx << std::endl;
2068 //std::cout << UVy << "," << UUy << std::endl;
2069
2070 delta_x = - UVx / UUx;
2071 delta_y = - UVy / UUy;
2072 }
2073
2084 inline void poc( const mist::array2< double > &input, const mist::array2< double > &reference, double &delta_x, double &delta_y, const double lowpass_range = 0.85, const int distance = 5 )
2085 {
2086 int width = reference.width();
2087 int height = reference.height();
2088 mist::array2< std::complex< double > > freq_ref( width , height ), freq_input( width , height ) ; //周波数空間画像
2089 //mist::array2< std::complex< double > > phase_ref( width , height ), phase_input( width , height ) ; //位相限定画像
2090 mist::array2< std::complex< double > > freq_poc( width , height ); //
2091 mist::array2< double > poc_image;
2092 double norm = 0.0;
2093
2094 // FFT
2095 mist::fft( reference, freq_ref );
2096 mist::fft( input, freq_input );
2097 // 位相限定相関画像作成 --------------------------------------
2098 //for( int h = 0 ; h < height ; h++)
2099 //{
2100 // for( int w = 0 ; w < width ; w++ )
2101 // {
2102 // freq_poc( w, h ) = freq_ref( w , h ) * conj( freq_input( w , h ) );
2103 // norm = abs( freq_poc( w, h) );
2104 // if( norm == 0){
2105 // freq_poc( w, h ) = freq_poc( w, h );
2106 // }
2107 // else
2108 // {
2109 // freq_poc( w, h ) = freq_poc( w, h ) / norm;
2110 // }
2111 // }
2112 //}
2113 for( size_t i = 0 ; i < reference.size() ; i++)
2114 {
2115 freq_poc[i] = freq_ref[i] * conj( freq_input[i] );
2116 norm = abs( freq_poc[i] );
2117 if( norm != 0)
2118 {
2119 freq_poc[i] = freq_poc[i] / norm;
2120 }
2121 }
2122 // end 位相限定相関画像計算 --------------------------------------
2123
2124 // ローパスフィルタ
2125 freq_poc = lowpass_filter( freq_poc, lowpass_range );
2126
2127 // 逆FFT
2128 mist::ifft( freq_poc, poc_image );
2129
2130 // 象限の配置換え(1と3、2と4) ----------------------------
2131 poc_image = shuffle_image( poc_image );
2132
2133 //ピーク、サブピーク位置の探索 ----------------------------
2134 point< double > peak;
2135 point< double > x_peak, y_peak;
2136 serch_peak( poc_image, peak, x_peak, y_peak );
2137
2138 //std::cout << "peak ( " << peak.x - width / 2 << ", " << peak.y - height / 2 << " ) " << peak.value << std::endl;
2139 //std::cout << "x_peak ( " << x_peak.x - width / 2 << ", " << x_peak.y - height / 2 << " ) " << x_peak.value << std::endl;
2140 //std::cout << "y_peak ( " << y_peak.x - width / 2 << ", " << y_peak.y - height / 2 << " ) " << y_peak.value << std::endl;
2141 //PEF
2142 PEF( poc_image, peak, x_peak, y_peak, delta_x, delta_y, distance, lowpass_range );
2143
2144 delta_x *= -1;
2145 delta_y *= -1;
2146 }
2147
2148
2149 } // 名前空間 _poc_ の終わり
2150
2151 namespace poc
2152 {
2153 enum interpolate_type{ LINEAR, BICUBIC };
2166 template < typename Value_type, typename Allocator >
2167 void estimate(
2168 const array2< Value_type, Allocator > &input,
2169 const array2< Value_type, Allocator > &reference,
2170 double &dx,
2171 double &dy,
2172 const double low_pass_range = 0.85,
2173 const int distance = 5)
2174 {
2175 typedef array2< double > image_type;
2176 image_type tmp1, tmp2;
2177 convert( input, tmp1 );
2178 convert( reference, tmp2 );
2179 return _poc_::poc( tmp1, tmp2, dx, dy, low_pass_range, distance );
2180 }
2192 template < typename Value_type, typename Allocator >
2193 void estimate(
2194 const array2< Value_type, Allocator > &input,
2195 const array2< Value_type, Allocator > &reference,
2196 vector2< double > &dvec,
2197 const double low_pass_range = 0.85,
2198 const int distance = 5)
2199 {
2200 typedef array2< double > image_type;
2201 image_type tmp1, tmp2;
2202 convert( input, tmp1 );
2203 convert( reference, tmp2 );
2204 return _poc_::poc( tmp1, tmp2, dvec.x, dvec.y, low_pass_range, distance );
2205 }
2206
2217 template< class Value_type1, class Allocator1, class Value_type2, class Allocator2 >
2218 void transform(
2219 const array2< Value_type1, Allocator1 > &in,
2220 array2< Value_type2, Allocator2 > &out,
2221 const double &dx,
2222 const double &dy,
2223 const double factor = -1,
2224 const int interpolate_type = BICUBIC)
2225 {
2226 if( factor != -1 )
2227 {
2228 out.resize( (size_t)( in.width() * factor ), (size_t )( in.height() * factor ) );
2229 }
2230 switch( interpolate_type )
2231 {
2232 case LINEAR:
2233 return linear::transform( in, out, -dx, -dy );
2234 case BICUBIC:
2235 return cubic::transform( in, out, -dx, -dy );
2236 }
2237 return;
2238 }
2239
2250 template< class Value_type1, class Allocator1, class Value_type2, class Allocator2 >
2251 void transform(
2252 const array2< Value_type1, Allocator1 > &in,
2253 array2< Value_type2, Allocator2 > &out,
2254 const vector2< double > &dvec,
2255 const double factor = -1,
2256 const int interpolate_type = BICUBIC)
2257 {
2258 if( factor != -1 )
2259 {
2260 out.resize( (size_t)( in.width() * factor ), (size_t )( in.height() * factor ) );
2261 }
2262 switch( interpolate_type )
2263 {
2264 case LINEAR:
2265 return linear::transform( in, out, -dvec.x, -dvec.y );
2266 case BICUBIC:
2267 return cubic::transform( in, out, -dvec.x, -dvec.y );
2268 }
2269 return;
2270 }
2271
2272 } // 名前空間 poc の終わり
2274 // 画像のレジストレーショングループの終わり
2275
2276
2277
2278 // mist名前空間の終わり
2279 _MIST_END
2280
2281
2282 #endif // __INCLUDE_MIST_REGISTRATION__
2283