1 /*M///////////////////////////////////////////////////////////////////////////////////////
2 //
3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
4 //
5 // By downloading, copying, installing or using the software you agree to this license.
6 // If you do not agree to this license, do not download, install,
7 // copy or use the software.
8 //
9 //
10 // Intel License Agreement
11 // For Open Source Computer Vision Library
12 //
13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
14 // Third party copyrights are property of their respective owners.
15 //
16 // Redistribution and use in source and binary forms, with or without modification,
17 // are permitted provided that the following conditions are met:
18 //
19 // * Redistribution's of source code must retain the above copyright notice,
20 // this list of conditions and the following disclaimer.
21 //
22 // * Redistribution's in binary form must reproduce the above copyright notice,
23 // this list of conditions and the following disclaimer in the documentation
24 // and/or other materials provided with the distribution.
25 //
26 // * The name of Intel Corporation may not be used to endorse or promote products
27 // derived from this software without specific prior written permission.
28 //
29 // This software is provided by the copyright holders and contributors "as is" and
30 // any express or implied warranties, including, but not limited to, the implied
31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
32 // In no event shall the Intel Corporation or contributors be liable for any direct,
33 // indirect, incidental, special, exemplary, or consequential damages
34 // (including, but not limited to, procurement of substitute goods or services;
35 // loss of use, data, or profits; or business interruption) however caused
36 // and on any theory of liability, whether in contract, strict liability,
37 // or tort (including negligence or otherwise) arising in any way out of
38 // the use of this software, even if advised of the possibility of such damage.
39 //
40 //M*/
41 #include "_cvaux.h"
42
43 //*F///////////////////////////////////////////////////////////////////////////////////////
44 // Name: icvImgToObs_DCT_8u32f_C1R
45 // Purpose: The function takes as input an image and returns the sequnce of observations
46 // to be used with an embedded HMM; Each observation is top-left block of DCT
47 // coefficient matrix.
48 // Context:
49 // Parameters: img - pointer to the original image ROI
50 // imgStep - full row width of the image in bytes
51 // roi - width and height of ROI in pixels
52 // obs - pointer to resultant observation vectors
53 // dctSize - size of the block for which DCT is calculated
54 // obsSize - size of top-left block of DCT coeffs matrix, which is treated
55 // as observation. Each observation vector consists of
56 // obsSize.width * obsSize.height floats.
57 // The following conditions should be satisfied:
58 // 0 < objSize.width <= dctSize.width,
59 // 0 < objSize.height <= dctSize.height.
60 // delta - dctBlocks are overlapped and this parameter specifies horizontal
61 // and vertical shift.
62 // Returns:
63 // CV_NO_ERR or error code
64 // Notes:
65 // The algorithm is following:
66 // 1. First, number of observation vectors per row and per column are calculated:
67 //
68 // Nx = floor((roi.width - dctSize.width + delta.width)/delta.width);
69 // Ny = floor((roi.height - dctSize.height + delta.height)/delta.height);
70 //
71 // So, total number of observation vectors is Nx*Ny, and total size of
72 // array obs must be >= Nx*Ny*obsSize.width*obsSize.height*sizeof(float).
73 // 2. Observation vectors are calculated in the following loop
74 // ( actual implementation may be different ), where
75 // I[x1:x2,y1:y2] means block of pixels from source image with
76 // x1 <= x < x2, y1 <= y < y2,
77 // D[x1:x2,y1:y2] means sub matrix of DCT matrix D.
78 // O[x,y] means observation vector that corresponds to position
79 // (x*delta.width,y*delta.height) in the source image
80 // ( all indices are counted from 0 ).
81 //
82 // for( y = 0; y < Ny; y++ )
83 // {
84 // for( x = 0; x < Nx; x++ )
85 // {
86 // D = DCT(I[x*delta.width : x*delta.width + dctSize.width,
87 // y*delta.height : y*delta.height + dctSize.height]);
88 // O[x,y] = D[0:obsSize.width, 0:obsSize.height];
89 // }
90 // }
91 //F*/
92
93 /*comment out the following line to make DCT be calculated in floating-point arithmetics*/
94 //#define _CV_INT_DCT
95
96 /* for integer DCT only */
97 #define DCT_SCALE 15
98
99 #ifdef _CV_INT_DCT
100 typedef int work_t;
101
102 #define DESCALE CV_DESCALE
103 #define SCALE(x) CV_FLT_TO_FIX((x),DCT_SCALE)
104 #else
105 typedef float work_t;
106
107 #define DESCALE(x,n) (float)(x)
108 #define SCALE(x) (float)(x)
109 #endif
110
111 /* calculate dct transform matrix */
112 static void icvCalcDCTMatrix( work_t * cfs, int n );
113
114 #define MAX_DCT_SIZE 32
115
116 static CvStatus CV_STDCALL
icvImgToObs_DCT_8u32f_C1R(uchar * img,int imgStep,CvSize roi,float * obs,CvSize dctSize,CvSize obsSize,CvSize delta)117 icvImgToObs_DCT_8u32f_C1R( uchar * img, int imgStep, CvSize roi,
118 float *obs, CvSize dctSize,
119 CvSize obsSize, CvSize delta )
120 {
121 /* dct transform matrices: horizontal and vertical */
122 work_t tab_x[MAX_DCT_SIZE * MAX_DCT_SIZE / 2 + 2];
123 work_t tab_y[MAX_DCT_SIZE * MAX_DCT_SIZE / 2 + 2];
124
125 /* temporary buffers for dct */
126 work_t temp0[MAX_DCT_SIZE * 4];
127 work_t temp1[MAX_DCT_SIZE * 4];
128 work_t *buffer = 0;
129 work_t *buf_limit;
130
131 double s;
132
133 int y;
134 int Nx, Ny;
135
136 int n1 = dctSize.height, m1 = n1 / 2;
137 int n2 = dctSize.width, m2 = n2 / 2;
138
139 if( !img || !obs )
140 return CV_NULLPTR_ERR;
141
142 if( roi.width <= 0 || roi.height <= 0 )
143 return CV_BADSIZE_ERR;
144
145 if( delta.width <= 0 || delta.height <= 0 )
146 return CV_BADRANGE_ERR;
147
148 if( obsSize.width <= 0 || dctSize.width < obsSize.width ||
149 obsSize.height <= 0 || dctSize.height < obsSize.height )
150 return CV_BADRANGE_ERR;
151
152 if( dctSize.width > MAX_DCT_SIZE || dctSize.height > MAX_DCT_SIZE )
153 return CV_BADRANGE_ERR;
154
155 Nx = (roi.width - dctSize.width + delta.width) / delta.width;
156 Ny = (roi.height - dctSize.height + delta.height) / delta.height;
157
158 if( Nx <= 0 || Ny <= 0 )
159 return CV_BADRANGE_ERR;
160
161 buffer = (work_t *)cvAlloc( roi.width * obsSize.height * sizeof( buffer[0] ));
162 if( !buffer )
163 return CV_OUTOFMEM_ERR;
164
165 icvCalcDCTMatrix( tab_x, dctSize.width );
166 icvCalcDCTMatrix( tab_y, dctSize.height );
167
168 buf_limit = buffer + obsSize.height * roi.width;
169
170 for( y = 0; y < Ny; y++, img += delta.height * imgStep )
171 {
172 int x, i, j, k;
173 work_t k0 = 0;
174
175 /* do transfroms for each column. Calc only first obsSize.height DCT coefficients */
176 for( x = 0; x < roi.width; x++ )
177 {
178 float is = 0;
179 work_t *buf = buffer + x;
180 work_t *tab = tab_y + 2;
181
182 if( n1 & 1 )
183 {
184 is = img[x + m1 * imgStep];
185 k0 = ((work_t) is) * tab[-1];
186 }
187
188 /* first coefficient */
189 for( j = 0; j < m1; j++ )
190 {
191 float t0 = img[x + j * imgStep];
192 float t1 = img[x + (n1 - 1 - j) * imgStep];
193 float t2 = t0 + t1;
194
195 t0 -= t1;
196 temp0[j] = (work_t) t2;
197 is += t2;
198 temp1[j] = (work_t) t0;
199 }
200
201 buf[0] = DESCALE( is * tab[-2], PASS1_SHIFT );
202 if( (buf += roi.width) >= buf_limit )
203 continue;
204
205 /* other coefficients */
206 for( ;; )
207 {
208 s = 0;
209
210 for( k = 0; k < m1; k++ )
211 s += temp1[k] * tab[k];
212
213 buf[0] = DESCALE( s, PASS1_SHIFT );
214 if( (buf += roi.width) >= buf_limit )
215 break;
216
217 tab += m1;
218 s = 0;
219
220 if( n1 & 1 )
221 {
222 k0 = -k0;
223 s = k0;
224 }
225 for( k = 0; k < m1; k++ )
226 s += temp0[k] * tab[k];
227
228 buf[0] = DESCALE( s, PASS1_SHIFT );
229 tab += m1;
230
231 if( (buf += roi.width) >= buf_limit )
232 break;
233 }
234 }
235
236 k0 = 0;
237
238 /* do transforms for rows. */
239 for( x = 0; x + dctSize.width <= roi.width; x += delta.width )
240 {
241 for( i = 0; i < obsSize.height; i++ )
242 {
243 work_t *buf = buffer + x + roi.width * i;
244 work_t *tab = tab_x + 2;
245 float *obs_limit = obs + obsSize.width;
246
247 s = 0;
248
249 if( n2 & 1 )
250 {
251 s = buf[m2];
252 k0 = (work_t) (s * tab[-1]);
253 }
254
255 /* first coefficient */
256 for( j = 0; j < m2; j++ )
257 {
258 work_t t0 = buf[j];
259 work_t t1 = buf[n2 - 1 - j];
260 work_t t2 = t0 + t1;
261
262 t0 -= t1;
263 temp0[j] = (work_t) t2;
264 s += t2;
265 temp1[j] = (work_t) t0;
266 }
267
268 *obs++ = (float) DESCALE( s * tab[-2], PASS2_SHIFT );
269
270 if( obs == obs_limit )
271 continue;
272
273 /* other coefficients */
274 for( ;; )
275 {
276 s = 0;
277
278 for( k = 0; k < m2; k++ )
279 s += temp1[k] * tab[k];
280
281 obs[0] = (float) DESCALE( s, PASS2_SHIFT );
282 if( ++obs == obs_limit )
283 break;
284
285 tab += m2;
286
287 s = 0;
288
289 if( n2 & 1 )
290 {
291 k0 = -k0;
292 s = k0;
293 }
294 for( k = 0; k < m2; k++ )
295 s += temp0[k] * tab[k];
296 obs[0] = (float) DESCALE( s, PASS2_SHIFT );
297
298 tab += m2;
299 if( ++obs == obs_limit )
300 break;
301 }
302 }
303 }
304 }
305
306 cvFree( &buffer );
307 return CV_NO_ERR;
308 }
309
310
311 static CvStatus CV_STDCALL
icvImgToObs_DCT_32f_C1R(float * img,int imgStep,CvSize roi,float * obs,CvSize dctSize,CvSize obsSize,CvSize delta)312 icvImgToObs_DCT_32f_C1R( float * img, int imgStep, CvSize roi,
313 float *obs, CvSize dctSize,
314 CvSize obsSize, CvSize delta )
315 {
316 /* dct transform matrices: horizontal and vertical */
317 work_t tab_x[MAX_DCT_SIZE * MAX_DCT_SIZE / 2 + 2];
318 work_t tab_y[MAX_DCT_SIZE * MAX_DCT_SIZE / 2 + 2];
319
320 /* temporary buffers for dct */
321 work_t temp0[MAX_DCT_SIZE * 4];
322 work_t temp1[MAX_DCT_SIZE * 4];
323 work_t *buffer = 0;
324 work_t *buf_limit;
325
326 double s;
327
328 int y;
329 int Nx, Ny;
330
331 int n1 = dctSize.height, m1 = n1 / 2;
332 int n2 = dctSize.width, m2 = n2 / 2;
333
334 if( !img || !obs )
335 return CV_NULLPTR_ERR;
336
337 if( roi.width <= 0 || roi.height <= 0 )
338 return CV_BADSIZE_ERR;
339
340 if( delta.width <= 0 || delta.height <= 0 )
341 return CV_BADRANGE_ERR;
342
343 if( obsSize.width <= 0 || dctSize.width < obsSize.width ||
344 obsSize.height <= 0 || dctSize.height < obsSize.height )
345 return CV_BADRANGE_ERR;
346
347 if( dctSize.width > MAX_DCT_SIZE || dctSize.height > MAX_DCT_SIZE )
348 return CV_BADRANGE_ERR;
349
350 Nx = (roi.width - dctSize.width + delta.width) / delta.width;
351 Ny = (roi.height - dctSize.height + delta.height) / delta.height;
352
353 if( Nx <= 0 || Ny <= 0 )
354 return CV_BADRANGE_ERR;
355
356 buffer = (work_t *)cvAlloc( roi.width * obsSize.height * sizeof( buffer[0] ));
357 if( !buffer )
358 return CV_OUTOFMEM_ERR;
359
360 icvCalcDCTMatrix( tab_x, dctSize.width );
361 icvCalcDCTMatrix( tab_y, dctSize.height );
362
363 buf_limit = buffer + obsSize.height * roi.width;
364
365 imgStep /= sizeof(img[0]);
366
367 for( y = 0; y < Ny; y++, img += delta.height * imgStep )
368 {
369 int x, i, j, k;
370 work_t k0 = 0;
371
372 /* do transfroms for each column. Calc only first obsSize.height DCT coefficients */
373 for( x = 0; x < roi.width; x++ )
374 {
375 float is = 0;
376 work_t *buf = buffer + x;
377 work_t *tab = tab_y + 2;
378
379 if( n1 & 1 )
380 {
381 is = img[x + m1 * imgStep];
382 k0 = ((work_t) is) * tab[-1];
383 }
384
385 /* first coefficient */
386 for( j = 0; j < m1; j++ )
387 {
388 float t0 = img[x + j * imgStep];
389 float t1 = img[x + (n1 - 1 - j) * imgStep];
390 float t2 = t0 + t1;
391
392 t0 -= t1;
393 temp0[j] = (work_t) t2;
394 is += t2;
395 temp1[j] = (work_t) t0;
396 }
397
398 buf[0] = DESCALE( is * tab[-2], PASS1_SHIFT );
399 if( (buf += roi.width) >= buf_limit )
400 continue;
401
402 /* other coefficients */
403 for( ;; )
404 {
405 s = 0;
406
407 for( k = 0; k < m1; k++ )
408 s += temp1[k] * tab[k];
409
410 buf[0] = DESCALE( s, PASS1_SHIFT );
411 if( (buf += roi.width) >= buf_limit )
412 break;
413
414 tab += m1;
415 s = 0;
416
417 if( n1 & 1 )
418 {
419 k0 = -k0;
420 s = k0;
421 }
422 for( k = 0; k < m1; k++ )
423 s += temp0[k] * tab[k];
424
425 buf[0] = DESCALE( s, PASS1_SHIFT );
426 tab += m1;
427
428 if( (buf += roi.width) >= buf_limit )
429 break;
430 }
431 }
432
433 k0 = 0;
434
435 /* do transforms for rows. */
436 for( x = 0; x + dctSize.width <= roi.width; x += delta.width )
437 {
438 for( i = 0; i < obsSize.height; i++ )
439 {
440 work_t *buf = buffer + x + roi.width * i;
441 work_t *tab = tab_x + 2;
442 float *obs_limit = obs + obsSize.width;
443
444 s = 0;
445
446 if( n2 & 1 )
447 {
448 s = buf[m2];
449 k0 = (work_t) (s * tab[-1]);
450 }
451
452 /* first coefficient */
453 for( j = 0; j < m2; j++ )
454 {
455 work_t t0 = buf[j];
456 work_t t1 = buf[n2 - 1 - j];
457 work_t t2 = t0 + t1;
458
459 t0 -= t1;
460 temp0[j] = (work_t) t2;
461 s += t2;
462 temp1[j] = (work_t) t0;
463 }
464
465 *obs++ = (float) DESCALE( s * tab[-2], PASS2_SHIFT );
466
467 if( obs == obs_limit )
468 continue;
469
470 /* other coefficients */
471 for( ;; )
472 {
473 s = 0;
474
475 for( k = 0; k < m2; k++ )
476 s += temp1[k] * tab[k];
477
478 obs[0] = (float) DESCALE( s, PASS2_SHIFT );
479 if( ++obs == obs_limit )
480 break;
481
482 tab += m2;
483
484 s = 0;
485
486 if( n2 & 1 )
487 {
488 k0 = -k0;
489 s = k0;
490 }
491 for( k = 0; k < m2; k++ )
492 s += temp0[k] * tab[k];
493 obs[0] = (float) DESCALE( s, PASS2_SHIFT );
494
495 tab += m2;
496 if( ++obs == obs_limit )
497 break;
498 }
499 }
500 }
501 }
502
503 cvFree( &buffer );
504 return CV_NO_ERR;
505 }
506
507
508 static void
icvCalcDCTMatrix(work_t * cfs,int n)509 icvCalcDCTMatrix( work_t * cfs, int n )
510 {
511 static const double sqrt2 = 1.4142135623730950488016887242097;
512 static const double pi = 3.1415926535897932384626433832795;
513
514 static const double sincos[16 * 2] = {
515 1.00000000000000000, 0.00000000000000006,
516 0.70710678118654746, 0.70710678118654757,
517 0.49999999999999994, 0.86602540378443871,
518 0.38268343236508978, 0.92387953251128674,
519 0.30901699437494740, 0.95105651629515353,
520 0.25881904510252074, 0.96592582628906831,
521 0.22252093395631439, 0.97492791218182362,
522 0.19509032201612825, 0.98078528040323043,
523 0.17364817766693033, 0.98480775301220802,
524 0.15643446504023087, 0.98768834059513777,
525 0.14231483827328514, 0.98982144188093268,
526 0.13052619222005157, 0.99144486137381038,
527 0.12053668025532305, 0.99270887409805397,
528 0.11196447610330786, 0.99371220989324260,
529 0.10452846326765346, 0.99452189536827329,
530 0.09801714032956060, 0.99518472667219693,
531 };
532
533 #define ROTATE( c, s, dc, ds ) \
534 { \
535 t = c*dc - s*ds; \
536 s = c*ds + s*dc; \
537 c = t; \
538 }
539
540 #define WRITE2( j, a, b ) \
541 { \
542 cfs[j] = SCALE(a); \
543 cfs2[j] = SCALE(b); \
544 }
545
546 double t, scale = 1. / sqrt( (double)n );
547 int i, j, m = n / 2;
548
549 cfs[0] = SCALE( scale );
550 scale *= sqrt2;
551 cfs[1] = SCALE( scale );
552 cfs += 2 - m;
553
554 if( n > 1 )
555 {
556 double a0, b0;
557 double da0, db0;
558 work_t *cfs2 = cfs + m * n;
559
560 if( n <= 16 )
561 {
562 da0 = a0 = sincos[2 * n - 1];
563 db0 = b0 = sincos[2 * n - 2];
564 }
565 else
566 {
567 t = pi / (2 * n);
568 da0 = a0 = cos( t );
569 db0 = b0 = sin( t );
570 }
571
572 /* other rows */
573 for( i = 1; i <= m; i++ )
574 {
575 double a = a0 * scale;
576 double b = b0 * scale;
577 double da = a0 * a0 - b0 * b0;
578 double db = a0 * b0 + a0 * b0;
579
580 cfs += m;
581 cfs2 -= m;
582
583 for( j = 0; j < m; j += 2 )
584 {
585 WRITE2( j, a, b );
586 ROTATE( a, b, da, db );
587 if( j + 1 < m )
588 {
589 WRITE2( j + 1, a, -b );
590 ROTATE( a, b, da, db );
591 }
592 }
593
594 ROTATE( a0, b0, da0, db0 );
595 }
596 }
597 #undef ROTATE
598 #undef WRITE2
599 }
600
601
602 CV_IMPL void
cvImgToObs_DCT(const void * arr,float * obs,CvSize dctSize,CvSize obsSize,CvSize delta)603 cvImgToObs_DCT( const void* arr, float *obs, CvSize dctSize,
604 CvSize obsSize, CvSize delta )
605 {
606 CV_FUNCNAME( "cvImgToObs_DCT" );
607
608 __BEGIN__;
609
610 CvMat stub, *mat = (CvMat*)arr;
611
612 CV_CALL( mat = cvGetMat( arr, &stub ));
613
614 switch( CV_MAT_TYPE( mat->type ))
615 {
616 case CV_8UC1:
617 IPPI_CALL( icvImgToObs_DCT_8u32f_C1R( mat->data.ptr, mat->step,
618 cvGetMatSize(mat), obs,
619 dctSize, obsSize, delta ));
620 break;
621 case CV_32FC1:
622 IPPI_CALL( icvImgToObs_DCT_32f_C1R( mat->data.fl, mat->step,
623 cvGetMatSize(mat), obs,
624 dctSize, obsSize, delta ));
625 break;
626 default:
627 CV_ERROR( CV_StsUnsupportedFormat, "" );
628 }
629
630 __END__;
631 }
632
633
634 /* End of file. */
635