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 
Generated on Wed Nov 12 2014 19:44:22 for MIST by doxygen 1.8.1.2