dct.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 
37 
38 #ifndef __INCLUDE_DCT_H__
39 #define __INCLUDE_DCT_H__
40 
41 
42 #ifndef __INCLUDE_MIST_H__
43 #include "../mist.h"
44 #endif
45 
46 #ifndef __INCLUDE_FFT_UTIL_H__
47 #include "fft_util.h"
48 #endif
49 
50 
51 #if defined( __MIST_WINDOWS__ ) && __MIST_WINDOWS__ > 0
52 #define USE_CDFT_WINTHREADS
53 #define USE_FFT2D_WINTHREADS
54 #define USE_FFT3D_WINTHREADS
55 #else
56 #define USE_CDFT_PTHREADS
57 #define USE_FFT2D_PTHREADS
58 #define USE_FFT3D_PTHREADS
59 #endif
60 
61 #include "fftsg.h"
62 #include "fftsg2d.h"
63 #include "fftsg3d.h"
64 
65 // mist名前空間の始まり
66 _MIST_BEGIN
67 
68 
71 
79 
80 
92 template < class T1, class T2, class Allocator1, class Allocator2 >
93 bool dct( const array1< T1, Allocator1 > &in, array1< T2, Allocator2 > &out )
94 {
95  if( !__fft_util__::size_check( ( unsigned int ) in.size( ) ) )
96  {
97  return( false );
98  }
99 
100  typedef typename Allocator1::size_type size_type;
101  size_type i;
102  __fft_util__::FFT_MEMORY1 mem;
103 
104  if( !__fft_util__::allocate_memory( mem,
105  in.size( ),
106  static_cast< size_t >( std::sqrt( static_cast< double >( in.size( ) / 2 ) ) + 3 ),
107  in.size( ) * 5 / 4 ) )
108  {
109  __fft_util__::deallocate_memory( mem );
110  return( false );
111  }
112 
113  double *data = mem.data;
114  double *w = mem.w;
115  int *ip = mem.ip;
116 
117  for( i = 0 ; i < in.size( ) ; i++ )
118  {
119  std::complex< double > c( __fft_util__::convert_complex< T1 >::convert_to( in[ i ] ) );
120  data[ i ] = c.real( );
121  }
122 
123  ip[ 0 ] = 0;
124 
125  ooura_fft::ddct( static_cast<int>( in.size() ), -1, data, ip, w );
126 
127 
128  out.resize( in.size( ) );
129 
130  for( i = 0 ; i < out.size( ) ; i++ )
131  {
132  out[ i ] = __fft_util__::convert_complex< T2 >::convert_from( data[ i ], 0.0 );
133  }
134 
135  __fft_util__::deallocate_memory( mem );
136 
137  return( true );
138 }
139 
140 
141 
152 template < class T1, class T2, class Allocator1, class Allocator2 >
153 bool idct( const array1< T1, Allocator1 > &in, array1< T2, Allocator2 > &out )
154 {
155  if( !__fft_util__::size_check( ( unsigned int ) in.size( ) ) )
156  {
157  return( false );
158  }
159 
160  typedef typename Allocator1::size_type size_type;
161  size_type i;
162  __fft_util__::FFT_MEMORY1 mem;
163 
164  if( !__fft_util__::allocate_memory( mem,
165  in.size( ),
166  static_cast< size_t >( std::sqrt( static_cast< double >( in.size( ) / 2 ) ) + 3 ),
167  in.size( ) * 5 / 4 ) )
168  {
169  __fft_util__::deallocate_memory( mem );
170  return( false );
171  }
172 
173  double *data = mem.data;
174  double *w = mem.w;
175  int *ip = mem.ip;
176 
177  for( i = 0 ; i < in.size( ) ; i++ )
178  {
179  std::complex< double > c( __fft_util__::convert_complex< T1 >::convert_to( in[ i ] ) );
180  data[ i ] = c.real( );
181  }
182 
183  ip[ 0 ] = 0;
184  data[ 0 ] *= 0.5;
185 
186  ooura_fft::ddct( static_cast<int>( in.size() ), 1, data, ip, w );
187 
188 
189  out.resize( in.size( ) );
190 
191  for( i = 0 ; i < out.size( ) ; i++ )
192  {
193  out[ i ] = __fft_util__::convert_complex< T2 >::convert_from( data[ i ] * 2.0 / out.size( ), 0.0 );
194  }
195 
196  __fft_util__::deallocate_memory( mem );
197 
198  return( true );
199 }
200 
201 
202 
214 template < class T1, class T2, class Allocator1, class Allocator2 >
215 bool dct( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out )
216 {
217  if( !__fft_util__::size_check( ( unsigned int ) in.width( ) ) || !__fft_util__::size_check( ( unsigned int ) in.height( ) ) )
218  {
219  return( false );
220  }
221 
222  typedef typename Allocator1::size_type size_type;
223  size_type i, j, size = in.width( ) > in.height( ) ? in.width( ) : in.height( );
224  __fft_util__::FFT_MEMORY2 mem;
225 
226  if( !__fft_util__::allocate_memory( mem,
227  in.width( ),
228  in.height( ),
229  4 * in.width( ) * FFT2D_MAX_THREADS,
230  static_cast< size_t >( std::sqrt( static_cast< double >( size / 2 ) ) + 3 ),
231  size * 3 / 2 ) )
232  {
233  __fft_util__::deallocate_memory( mem );
234  return( false );
235  }
236 
237  double **data = mem.data;
238  double *t = mem.t;
239  double *w = mem.w;
240  int *ip = mem.ip;
241 
242  for( i = 0 ; i < in.width( ) ; i++ )
243  {
244  for( j = 0 ; j < in.height( ) ; j++ )
245  {
246  std::complex< double > c( __fft_util__::convert_complex< T1 >::convert_to( in( i, j ) ) );
247  data[ i ][ j ] = c.real( );
248  }
249  }
250 
251  ip[ 0 ] = 0;
252 
253  ooura_fft::ddct2d( static_cast<int>( in.width( ) ), static_cast<int>( in.height( ) ), -1, data, t, ip, w );
254 
255 
256  out.resize( in.width( ), in.height( ) );
257 
258  for( i = 0 ; i < out.width( ) ; i++ )
259  {
260  for( j = 0 ; j < out.height( ) ; j++ )
261  {
262  out( i, j ) = __fft_util__::convert_complex< T2 >::convert_from( data[ i ][ j ], 0.0 );
263  }
264  }
265 
266  __fft_util__::deallocate_memory( mem );
267 
268  return( true );
269 }
270 
271 
272 
273 
284 template < class T1, class T2, class Allocator1, class Allocator2 >
285 bool idct( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out )
286 {
287  if( !__fft_util__::size_check( ( unsigned int ) in.width( ) ) || !__fft_util__::size_check( ( unsigned int ) in.height( ) ) )
288  {
289  return( false );
290  }
291 
292  typedef typename Allocator1::size_type size_type;
293  size_type i, j, size = in.width( ) > in.height( ) ? in.width( ) : in.height( );
294  __fft_util__::FFT_MEMORY2 mem;
295 
296  if( !__fft_util__::allocate_memory( mem,
297  in.width( ),
298  in.height( ),
299  4 * in.width( ) * FFT2D_MAX_THREADS,
300  static_cast< size_t >( std::sqrt( static_cast< double >( size / 2 ) ) + 3 ),
301  size * 3 / 2 ) )
302  {
303  __fft_util__::deallocate_memory( mem );
304  return( false );
305  }
306 
307  double **data = mem.data;
308  double *t = mem.t;
309  double *w = mem.w;
310  int *ip = mem.ip;
311 
312  for( i = 0 ; i < in.width( ) ; i++ )
313  {
314  for( j = 0 ; j < in.height( ) ; j++ )
315  {
316  std::complex< double > c( __fft_util__::convert_complex< T1 >::convert_to( in( i, j ) ) );
317  data[ i ][ j ] = c.real( );
318  }
319  }
320 
321  ip[ 0 ] = 0;
322 
323 
324  for( i = 0 ; i < in.width( ) ; i++ )
325  {
326  data[ i ][ 0 ] *= 0.5;
327  }
328 
329  for( j = 0 ; j < in.height( ) ; j++ )
330  {
331  data[ 0 ][ j ] *= 0.5;
332  }
333 
334 
335  ooura_fft::ddct2d( static_cast< int >( in.width( ) ), static_cast< int >( in.height( ) ), 1, data, t, ip, w );
336 
337 
338  out.resize( in.width( ), in.height( ) );
339 
340  for( i = 0 ; i < out.width( ) ; i++ )
341  {
342  for( j = 0 ; j < out.height( ) ; j++ )
343  {
344  out( i, j ) = __fft_util__::convert_complex< T2 >::convert_from( data[ i ][ j ] * 4.0 / in.size( ), 0.0 );
345  }
346  }
347 
348  __fft_util__::deallocate_memory( mem );
349 
350  return( true );
351 }
352 
353 
354 
355 
367 template < class T1, class T2, class Allocator1, class Allocator2 >
368 bool dct( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out )
369 {
370  if( !__fft_util__::size_check( ( unsigned int ) in.width( ) ) ||
371  !__fft_util__::size_check( ( unsigned int ) in.height( ) ) ||
372  !__fft_util__::size_check( ( unsigned int ) in.depth( ) ) )
373  {
374  return( false );
375  }
376 
377  typedef typename Allocator1::size_type size_type;
378  size_type i, j, k, size = in.width( ) > in.height( ) ? in.width( ) : in.height( );
379  __fft_util__::FFT_MEMORY3 mem;
380 
381  if( !__fft_util__::allocate_memory( mem,
382  in.width( ),
383  in.height( ),
384  in.depth( ),
385  4 * size * FFT3D_MAX_THREADS,
386  static_cast< size_t >( std::sqrt( static_cast< double >( size / 2 ) ) + 3 ),
387  size * 3 / 2 ) )
388  {
389  __fft_util__::deallocate_memory( mem );
390  return( false );
391  }
392 
393  double ***data = mem.data;
394  double *t = mem.t;
395  double *w = mem.w;
396  int *ip = mem.ip;
397 
398  for( i = 0 ; i < in.width( ) ; i++ )
399  {
400  for( j = 0 ; j < in.height( ) ; j++ )
401  {
402  for( k = 0 ; k < in.depth( ) ; k++ )
403  {
404  std::complex< double > c( __fft_util__::convert_complex< T1 >::convert_to( in( i, j, k ) ) );
405  data[ i ][ j ][ k ] = c.real( );
406  }
407  }
408  }
409 
410  ip[ 0 ] = 0;
411 
412  ooura_fft::ddct3d( static_cast<int>( in.width( ) ), static_cast<int>( in.height( ) ), static_cast<int>( in.depth( ) ), -1, data, t, ip, w );
413 
414 
415  out.resize( in.width( ), in.height( ), in.depth( ) );
416 
417  for( i = 0 ; i < out.width( ) ; i++ )
418  {
419  for( j = 0 ; j < out.height( ) ; j++ )
420  {
421  for( k = 0 ; k < out.depth( ) ; k++ )
422  {
423  out( i, j, k ) = __fft_util__::convert_complex< T2 >::convert_from( data[ i ][ j ][ k ], 0.0 );
424  }
425  }
426  }
427 
428  __fft_util__::deallocate_memory( mem );
429 
430  return( true );
431 }
432 
433 
434 
445 template < class T1, class T2, class Allocator1, class Allocator2 >
446 bool idct( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out )
447 {
448  if( !__fft_util__::size_check( ( unsigned int ) in.width( ) ) ||
449  !__fft_util__::size_check( ( unsigned int ) in.height( ) ) ||
450  !__fft_util__::size_check( ( unsigned int ) in.depth( ) ) )
451  {
452  return( false );
453  }
454 
455  typedef typename Allocator1::size_type size_type;
456  size_type i, j, k, size = in.width( ) > in.height( ) ? in.width( ) : in.height( );
457  __fft_util__::FFT_MEMORY3 mem;
458 
459  if( !__fft_util__::allocate_memory( mem,
460  in.width( ),
461  in.height( ),
462  in.depth( ),
463  4 * size * FFT3D_MAX_THREADS,
464  static_cast< size_t >( std::sqrt( static_cast< double >( size / 2 ) ) + 3 ),
465  size * 3 / 2 ) )
466  {
467  __fft_util__::deallocate_memory( mem );
468  return( false );
469  }
470 
471  double ***data = mem.data;
472  double *t = mem.t;
473  double *w = mem.w;
474  int *ip = mem.ip;
475 
476  for( i = 0 ; i < in.width( ) ; i++ )
477  {
478  for( j = 0 ; j < in.height( ) ; j++ )
479  {
480  for( k = 0 ; k < in.depth( ) ; k++ )
481  {
482  std::complex< double > c( __fft_util__::convert_complex< T1 >::convert_to( in( i, j, k ) ) );
483  data[ i ][ j ][ k ] = c.real( );
484  }
485  }
486  }
487 
488 
489  for( i = 0 ; i < in.width( ) ; i++ )
490  {
491  for( j = 0 ; j < in.height( ) ; j++ )
492  {
493  data[ i ][ j ][ 0 ] *= 0.5;
494  }
495 
496  for( k = 0 ; k < in.depth( ) ; k++ )
497  {
498  data[ i ][ 0 ][ k ] *= 0.5;
499  }
500  }
501 
502  for( j = 0 ; j < in.height( ) ; j++ )
503  {
504  for( k = 0 ; k < in.depth( ) ; k++ )
505  {
506  data[ 0 ][ j ][ k ] *= 0.5;
507  }
508  }
509 
510  ip[ 0 ] = 0;
511 
512  ooura_fft::ddct3d( static_cast<int>( in.width( ) ), static_cast<int>( in.height( ) ), static_cast<int>( in.depth( ) ), 1, data, t, ip, w );
513 
514  out.resize( in.width( ), in.height( ), in.depth( ) );
515 
516  for( i = 0 ; i < out.width( ) ; i++ )
517  {
518  for( j = 0 ; j < out.height( ) ; j++ )
519  {
520  for( k = 0 ; k < out.depth( ) ; k++ )
521  {
522  out( i, j, k ) = __fft_util__::convert_complex< T2 >::convert_from( data[ i ][ j ][ k ] * 8.0 / in.size( ), 0.0 );
523  }
524  }
525  }
526 
527  __fft_util__::deallocate_memory( mem );
528 
529  return( true );
530 }
531 
533 // DCT グループの終わり
534 
536 // Fourier グループの終わり
537 
538 
539 _MIST_END
540 
541 #endif
542 
543 
Generated on Wed Nov 12 2014 19:44:13 for MIST by doxygen 1.8.1.2