1 /*************************************************************************
2 * recbuffer.c - tools for recursive filtering of 3D and 2D image buffers
3 *
4 * $Id$
5 *
6 * LICENSE:
7 * GPL v3.0 (see gpl-3.0.txt for details)
8 *
9 * DESCRIPTION:
10 *
11 * recursive filtering of a buffer (a [1,2,3]D array)
12 * according that the filtering is separable
13 *
14 *
15 *
16 * AUTHOR:
17 * Gregoire Malandain (gregoire.malandain@inria.fr)
18 *
19 * CREATION DATE:
20 * June, 9 1998
21 *
22 * ADDITIONS, CHANGES
23 *
24 * * Jul 6 1999 (Gregoire Malandain)
25 * a bug in RecursiveFilterOnBuffer (*&%^@$^ cut and paste)
26 *
27 */
28
29 #include <stdio.h>
30 #include <stdlib.h>
31 #include <math.h>
32
33
34
35 #include <convert.h>
36 #include <recbuffer.h>
37
38 static int _VERBOSE_ = 0;
39
40 #define EXIT_ON_FAILURE 0
41 #define EXIT_ON_SUCCESS 1
42
43
44
45
46
47
48
49
50
51
52 /*
53 *
54 * Gradient modulus
55 *
56 *
57 */
58
GradientModulus(void * bufferIn,bufferType typeIn,void * bufferOut,bufferType typeOut,int * bufferDims,int * borderLengths,float * filterCoefs,recursiveFilterType filterType)59 int GradientModulus( void *bufferIn,
60 bufferType typeIn,
61 void *bufferOut,
62 bufferType typeOut,
63 int *bufferDims,
64 int *borderLengths,
65 float *filterCoefs,
66 recursiveFilterType filterType )
67 {
68 char *proc = "GradientModulus";
69 float *auxBuf = NULL;
70 float *tmpBuf = NULL, *grdBuf = NULL;
71 int sizeAuxBuf = 0;
72 int derivatives[3];
73 int i;
74
75
76 sizeAuxBuf = bufferDims[0] * bufferDims[1] * bufferDims[2];
77 if ( typeOut != FLOAT || bufferIn == bufferOut )
78 sizeAuxBuf *= 2;
79
80
81 /* allocation des buffers de calcul
82 */
83 auxBuf = (float*)malloc( sizeAuxBuf * sizeof(float) );
84 if ( auxBuf == NULL ) {
85 if ( _VERBOSE_ > 0 )
86 fprintf( stderr, "%s: unable to allocate auxiliary buffer\n", proc );
87 return( EXIT_ON_FAILURE );
88 }
89 tmpBuf = auxBuf;
90 if ( typeOut != FLOAT || bufferIn == bufferOut ) {
91 grdBuf = tmpBuf;
92 grdBuf += bufferDims[0] * bufferDims[1] * bufferDims[2];
93 } else {
94 grdBuf = (float*)bufferOut;
95 }
96
97 /* cas 2D
98 */
99 if ( bufferDims[2] == 1 ) {
100
101 derivatives[0] = DERIVATIVE_1;
102 derivatives[1] = DERIVATIVE_0;
103 derivatives[2] = NODERIVATIVE;
104 if ( RecursiveFilterOnBuffer( bufferIn, typeIn, (void*)grdBuf, FLOAT,
105 bufferDims, borderLengths, derivatives,
106 filterCoefs, filterType ) != EXIT_ON_SUCCESS ) {
107 if ( _VERBOSE_ )
108 fprintf( stderr, "%s: unable to compute X derivative (2D)\n", proc );
109 free( auxBuf );
110 return( EXIT_ON_FAILURE );
111 }
112
113 derivatives[0] = DERIVATIVE_0;
114 derivatives[1] = DERIVATIVE_1;
115 derivatives[2] = NODERIVATIVE;
116 if ( RecursiveFilterOnBuffer( bufferIn, typeIn, (void*)tmpBuf, FLOAT,
117 bufferDims, borderLengths, derivatives,
118 filterCoefs, filterType ) != EXIT_ON_SUCCESS ) {
119 if ( _VERBOSE_ )
120 fprintf( stderr, "%s: unable to compute Y derivative (2D)\n", proc );
121 free( auxBuf );
122 return( EXIT_ON_FAILURE );
123 }
124
125 sizeAuxBuf = bufferDims[0] * bufferDims[1] * bufferDims[2];
126 for ( i = 0; i < sizeAuxBuf; i++ )
127 grdBuf[i] = (float)sqrt( grdBuf[i]*grdBuf[i] + tmpBuf[i]*tmpBuf[i] );
128
129 } else {
130
131 derivatives[0] = NODERIVATIVE;
132 derivatives[1] = NODERIVATIVE;
133 derivatives[2] = DERIVATIVE_0;
134 if ( RecursiveFilterOnBuffer( bufferIn, typeIn, (void*)tmpBuf, FLOAT,
135 bufferDims, borderLengths, derivatives,
136 filterCoefs, filterType ) != EXIT_ON_SUCCESS ) {
137 if ( _VERBOSE_ )
138 fprintf( stderr, "%s: unable to compute Z smoothing (3D)\n", proc );
139 free( auxBuf );
140 return( EXIT_ON_FAILURE );
141 }
142
143 derivatives[0] = DERIVATIVE_1;
144 derivatives[1] = DERIVATIVE_0;
145 derivatives[2] = NODERIVATIVE;
146 if ( RecursiveFilterOnBuffer( (void*)tmpBuf, FLOAT, (void*)grdBuf, FLOAT,
147 bufferDims, borderLengths, derivatives,
148 filterCoefs, filterType ) != EXIT_ON_SUCCESS ) {
149 if ( _VERBOSE_ )
150 fprintf( stderr, "%s: unable to compute X derivative (3D)\n", proc );
151 free( auxBuf );
152 return( EXIT_ON_FAILURE );
153 }
154
155 derivatives[0] = DERIVATIVE_0;
156 derivatives[1] = DERIVATIVE_1;
157 derivatives[2] = NODERIVATIVE;
158 if ( RecursiveFilterOnBuffer( (void*)tmpBuf, FLOAT, (void*)tmpBuf, FLOAT,
159 bufferDims, borderLengths, derivatives,
160 filterCoefs, filterType ) != EXIT_ON_SUCCESS ) {
161 if ( _VERBOSE_ )
162 fprintf( stderr, "%s: unable to compute Y derivative (3D)\n", proc );
163 free( auxBuf );
164 return( EXIT_ON_FAILURE );
165 }
166
167 sizeAuxBuf = bufferDims[0] * bufferDims[1] * bufferDims[2];
168 for ( i = 0; i < sizeAuxBuf; i++ )
169 grdBuf[i] = grdBuf[i]*grdBuf[i] + tmpBuf[i]*tmpBuf[i];
170
171 derivatives[0] = DERIVATIVE_0;
172 derivatives[1] = DERIVATIVE_0;
173 derivatives[2] = DERIVATIVE_1;
174 if ( RecursiveFilterOnBuffer( bufferIn, typeIn, (void*)tmpBuf, FLOAT,
175 bufferDims, borderLengths, derivatives,
176 filterCoefs, filterType ) != EXIT_ON_SUCCESS ) {
177 if ( _VERBOSE_ )
178 fprintf( stderr, "%s: unable to compute Z derivative (3D)\n", proc );
179 free( auxBuf );
180 return( EXIT_ON_FAILURE );
181 }
182
183 for ( i = 0; i < sizeAuxBuf; i++ )
184 grdBuf[i] = (float)sqrt( grdBuf[i] + tmpBuf[i]*tmpBuf[i] );
185
186 }
187
188 if ( grdBuf != bufferOut )
189 ConvertBuffer( grdBuf, FLOAT, bufferOut, typeOut,
190 bufferDims[0]*bufferDims[1]*bufferDims[2] );
191 free( auxBuf );
192 return( EXIT_ON_SUCCESS );
193 }
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218 /*
219 *
220 * Laplacian
221 *
222 *
223 */
Laplacian_2D(void * bufferIn,bufferType typeIn,void * bufferOut,bufferType typeOut,int * bufferDims,int * borderLengths,float * filterCoefs,recursiveFilterType filterType)224 int Laplacian_2D ( void *bufferIn,
225 bufferType typeIn,
226 void *bufferOut,
227 bufferType typeOut,
228 int *bufferDims,
229 int *borderLengths,
230 float *filterCoefs,
231 recursiveFilterType filterType )
232 {
233 char *proc = "Laplacian_2D";
234 float *theXX = NULL;
235 float *theYY = NULL;
236
237 derivativeOrder XXderiv[3] = { DERIVATIVE_2, SMOOTHING, NODERIVATIVE };
238 derivativeOrder YYderiv[3] = { SMOOTHING, DERIVATIVE_2, NODERIVATIVE };
239 int sliceDims[3];
240 int z, i, dimxXdimy;
241
242 void *sliceOut = NULL;
243
244
245
246 /*
247 * We check the buffers' dimensions.
248 */
249 if ( (bufferDims[0] <= 0) || (bufferDims[1] <= 0) || (bufferDims[2] <= 0) ) {
250 if ( _VERBOSE_ > 0 )
251 fprintf( stderr, " Fatal error in %s: improper buffer's dimension.\n", proc );
252 return( EXIT_ON_FAILURE );
253 }
254
255 /*
256 * test of the coefficients
257 */
258 if ( (filterCoefs[0] < 0.0) || (filterCoefs[1] < 0.0) ||
259 (filterCoefs[2] < 0.0) ) {
260 if ( _VERBOSE_ > 0 )
261 fprintf( stderr, " Error in %s: negative coefficient's value.\n", proc );
262 return( EXIT_ON_FAILURE );
263 }
264
265
266 /*
267 *
268 */
269 dimxXdimy = bufferDims[0] * bufferDims[1];
270 sliceDims[0] = bufferDims[0];
271 sliceDims[1] = bufferDims[1];
272 sliceDims[2] = 1;
273
274
275 if ( typeOut == FLOAT ) {
276 theXX = (float*)malloc( dimxXdimy * sizeof( float ) );
277 } else {
278 theXX = (float*)malloc( 2 * dimxXdimy * sizeof( float ) );
279 }
280
281 if ( theXX == NULL ) {
282 if ( _VERBOSE_ > 0 ) {
283 fprintf( stderr, " Fatal error in %s:", proc );
284 fprintf( stderr, " unable to allocate auxiliary buffer.\n" );
285 }
286 return( EXIT_ON_FAILURE );
287 }
288
289 if ( typeOut != FLOAT ) {
290 theYY = theXX;
291 theYY += dimxXdimy;
292 }
293
294
295
296 for ( z=0; z<bufferDims[2]; z++ ) {
297
298 if ( typeOut == FLOAT ) {
299 theYY = ((float*)bufferOut) + z * dimxXdimy;
300 }
301
302 if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theXX, FLOAT,
303 sliceDims, borderLengths,
304 XXderiv, filterCoefs, filterType ) == 0 ) {
305 if ( _VERBOSE_ > 0 ) {
306 fprintf( stderr, " Fatal error in %s:", proc );
307 fprintf( stderr, " unable to compute X^2 derivative.\n" );
308 }
309 free( theXX );
310 return( EXIT_ON_FAILURE );
311 }
312
313 if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theYY, FLOAT,
314 sliceDims, borderLengths,
315 YYderiv, filterCoefs, filterType ) == 0 ) {
316 if ( _VERBOSE_ > 0 ) {
317 fprintf( stderr, " Fatal error in %s:", proc );
318 fprintf( stderr, " unable to compute Y^2 derivative.\n" );
319 }
320 free( theXX );
321 return( EXIT_ON_FAILURE );
322 }
323
324
325 for ( i=0; i<dimxXdimy; i++ ) theYY[i] += theXX[i];
326 if ( typeOut != FLOAT ) {
327 switch ( typeOut ) {
328 case UCHAR :
329 sliceOut = (((u8*)bufferOut) + z * dimxXdimy);
330 break;
331 case SCHAR :
332 sliceOut = (((s8*)bufferOut) + z * dimxXdimy);
333 break;
334 case SSHORT :
335 sliceOut = (((s16*)bufferOut) + z * dimxXdimy);
336 break;
337 case DOUBLE :
338 sliceOut = (((r64*)bufferOut) + z * dimxXdimy);
339 break;
340 default :
341 if ( _VERBOSE_ > 0 )
342 fprintf( stderr, " Error in %s: such output type not handled.\n", proc );
343 free( theXX );
344 return( EXIT_ON_FAILURE );
345 }
346 ConvertBuffer( theYY, FLOAT, sliceOut, typeOut, dimxXdimy );
347 }
348 }
349
350 return( EXIT_ON_SUCCESS );
351 }
352
353
354
355
356
357
358
359
360
361
362
363
364
365
Laplacian(void * bufferIn,bufferType typeIn,void * bufferOut,bufferType typeOut,int * bufferDims,int * borderLengths,float * filterCoefs,recursiveFilterType filterType)366 int Laplacian ( void *bufferIn,
367 bufferType typeIn,
368 void *bufferOut,
369 bufferType typeOut,
370 int *bufferDims,
371 int *borderLengths,
372 float *filterCoefs,
373 recursiveFilterType filterType )
374 {
375 char *proc = "Laplacian";
376 float *theSL = NULL;
377 float *theZZ = NULL;
378 float *theZ0 = NULL;
379
380
381 derivativeOrder XXderiv[3] = { DERIVATIVE_2, SMOOTHING, NODERIVATIVE };
382 derivativeOrder YYderiv[3] = { SMOOTHING, DERIVATIVE_2, NODERIVATIVE };
383 derivativeOrder Zsmooth[3] = { NODERIVATIVE, NODERIVATIVE, SMOOTHING };
384 derivativeOrder ZZderiv[3] = { SMOOTHING, SMOOTHING, DERIVATIVE_2 };
385
386 int sliceDims[3];
387 int z, i, j, dimxXdimy;
388
389
390
391
392 /*
393 * We check the buffers' dimensions.
394 */
395 if ( bufferDims[2] == 1 ) {
396 return( Laplacian_2D ( bufferIn, typeIn, bufferOut, typeOut,
397 bufferDims, borderLengths, filterCoefs, filterType ) );
398 }
399
400 if ( (bufferDims[0] <= 0) || (bufferDims[1] <= 0) || (bufferDims[2] <= 0) ) {
401 if ( _VERBOSE_ > 0 )
402 fprintf( stderr, " Fatal error in %s: improper buffer's dimension.\n", proc );
403 return( EXIT_ON_FAILURE );
404 }
405
406 /*
407 * test of the coefficients
408 */
409 if ( (filterCoefs[0] < 0.0) || (filterCoefs[1] < 0.0) ||
410 (filterCoefs[2] < 0.0) ) {
411 if ( _VERBOSE_ > 0 )
412 fprintf( stderr, " Error in %s: negative coefficient's value.\n", proc );
413 return( EXIT_ON_FAILURE );
414 }
415
416
417 /*
418 *
419 */
420 dimxXdimy = bufferDims[0] * bufferDims[1];
421 sliceDims[0] = bufferDims[0];
422 sliceDims[1] = bufferDims[1];
423 sliceDims[2] = 1;
424
425
426
427 if ( typeOut == FLOAT ) {
428 theSL = (float*)malloc( (1+bufferDims[2]) * dimxXdimy * sizeof( float ) );
429 } else {
430 theSL = (float*)malloc( (1+2*bufferDims[2]) * dimxXdimy * sizeof( float ) );
431 }
432
433
434
435 if ( theSL == NULL ) {
436 if ( _VERBOSE_ > 0 ) {
437 fprintf( stderr, " Fatal error in %s:", proc );
438 fprintf( stderr, " unable to allocate auxiliary buffer.\n" );
439 }
440 return( EXIT_ON_FAILURE );
441 }
442
443
444
445 theZ0 = theSL;
446 theZ0 += dimxXdimy;
447
448
449
450 if ( typeOut == FLOAT ) {
451 theZZ = bufferOut;
452 } else {
453 theZZ = theZ0;
454 theZZ += dimxXdimy * bufferDims[2];
455 }
456
457
458
459 /*
460 *
461 * 3D filtering / filtering along Z
462 *
463 */
464
465 if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theZ0, FLOAT,
466 bufferDims, borderLengths,
467 Zsmooth, filterCoefs, filterType ) == 0 ) {
468 if ( _VERBOSE_ > 0 ) {
469 fprintf( stderr, " Fatal error in %s:", proc );
470 fprintf( stderr, " unable to compute Z^0 derivative.\n" );
471 }
472 free( theSL );
473 return( EXIT_ON_FAILURE );
474 }
475 if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theZZ, FLOAT,
476 bufferDims, borderLengths,
477 ZZderiv, filterCoefs, filterType ) == 0 ) {
478 if ( _VERBOSE_ > 0 ) {
479 fprintf( stderr, " Fatal error in %s:", proc );
480 fprintf( stderr, " unable to compute Z^2 derivative.\n" );
481 }
482 free( theSL );
483 return( EXIT_ON_FAILURE );
484 }
485
486
487
488
489
490
491 for ( z=0; z<bufferDims[2]; z++ ) {
492
493 /*
494 *
495 * 2D filtering / filtering along X and Y
496 *
497 */
498
499 if ( RecursiveFilterOnBuffer( theZ0+z*dimxXdimy, FLOAT, theSL, FLOAT,
500 sliceDims, borderLengths,
501 XXderiv, filterCoefs, filterType ) == 0 ) {
502 if ( _VERBOSE_ > 0 ) {
503 fprintf( stderr, " Fatal error in %s:", proc );
504 fprintf( stderr, " unable to compute X^2 derivative.\n" );
505 }
506 free( theSL );
507 return( EXIT_ON_FAILURE );
508 }
509
510 for ( j=z*dimxXdimy, i=0; i<dimxXdimy; j++, i++ ) {
511 theZZ[j] += theSL[i];
512 }
513
514 if ( RecursiveFilterOnBuffer( theZ0+z*dimxXdimy, FLOAT, theSL, FLOAT,
515 sliceDims, borderLengths,
516 YYderiv, filterCoefs, filterType ) == 0 ) {
517 if ( _VERBOSE_ > 0 ) {
518 fprintf( stderr, " Fatal error in %s:", proc );
519 fprintf( stderr, " unable to compute Y^2 derivative.\n" );
520 }
521 free( theSL );
522 return( EXIT_ON_FAILURE );
523 }
524
525 for ( j=z*dimxXdimy, i=0; i<dimxXdimy; j++, i++ ) {
526 theZZ[j] += theSL[i];
527 }
528
529 }
530
531 if ( typeOut != FLOAT ) {
532 ConvertBuffer( theZZ, FLOAT, bufferOut, typeOut, bufferDims[2]*dimxXdimy );
533 }
534
535 return( EXIT_ON_SUCCESS );
536 }
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562 /*
563 *
564 * Gradient . Hessian * Gradient
565 *
566 *
567 */
GradientHessianGradient_2D(void * bufferIn,bufferType typeIn,void * bufferOut,bufferType typeOut,int * bufferDims,int * borderLengths,float * filterCoefs,recursiveFilterType filterType)568 int GradientHessianGradient_2D ( void *bufferIn,
569 bufferType typeIn,
570 void *bufferOut,
571 bufferType typeOut,
572 int *bufferDims,
573 int *borderLengths,
574 float *filterCoefs,
575 recursiveFilterType filterType )
576 {
577 char *proc = "GradientHessianGradient_2D";
578 float *theXX = NULL;
579 float *theYY = NULL;
580 float *theXY = NULL;
581 float *theX = NULL;
582 float *theY = NULL;
583
584 derivativeOrder Xsmooth[3] = { SMOOTHING, NODERIVATIVE, NODERIVATIVE };
585 derivativeOrder Yderiv[3] = { NODERIVATIVE, DERIVATIVE_1_EDGES, NODERIVATIVE };
586 derivativeOrder YYderiv[3] = { NODERIVATIVE, DERIVATIVE_2, NODERIVATIVE };
587
588 derivativeOrder Ysmooth[3] = { NODERIVATIVE, SMOOTHING, NODERIVATIVE };
589 derivativeOrder Xderiv[3] = { DERIVATIVE_1_EDGES, NODERIVATIVE, NODERIVATIVE };
590 derivativeOrder XXderiv[3] = { DERIVATIVE_2, NODERIVATIVE, NODERIVATIVE };
591
592 derivativeOrder XYderiv[3] = { DERIVATIVE_1, DERIVATIVE_1, NODERIVATIVE };
593
594 int sliceDims[3];
595 int z, i, dimxXdimy;
596
597 void *sliceIn = NULL;
598 void *sliceOut = NULL;
599
600 double gx, gy, g;
601
602 /*
603 * We check the buffers' dimensions.
604 */
605 if ( (bufferDims[0] <= 0) || (bufferDims[1] <= 0) || (bufferDims[2] <= 0) ) {
606 if ( _VERBOSE_ > 0 )
607 fprintf( stderr, " Fatal error in %s: improper buffer's dimension.\n", proc );
608 return( EXIT_ON_FAILURE );
609 }
610
611 /*
612 * test of the coefficients
613 */
614 if ( (filterCoefs[0] < 0.0) || (filterCoefs[1] < 0.0) ||
615 (filterCoefs[2] < 0.0) ) {
616 if ( _VERBOSE_ > 0 )
617 fprintf( stderr, " Error in %s: negative coefficient's value.\n", proc );
618 return( EXIT_ON_FAILURE );
619 }
620
621
622 /*
623 *
624 */
625 dimxXdimy = bufferDims[0] * bufferDims[1];
626 sliceDims[0] = bufferDims[0];
627 sliceDims[1] = bufferDims[1];
628 sliceDims[2] = 1;
629
630
631 if ( typeOut == FLOAT ) {
632 theXX = (float*)malloc( 4 * dimxXdimy * sizeof( float ) );
633 } else {
634 theXX = (float*)malloc( 5 * dimxXdimy * sizeof( float ) );
635 }
636
637 if ( theXX == NULL ) {
638 if ( _VERBOSE_ > 0 ) {
639 fprintf( stderr, " Fatal error in %s:", proc );
640 fprintf( stderr, " unable to allocate auxiliary buffer.\n" );
641 }
642 return( EXIT_ON_FAILURE );
643 }
644
645
646
647 theX = theY = theYY = theXX;
648 theYY += dimxXdimy;
649 theX += 2*dimxXdimy;
650 theY += 3*dimxXdimy;
651
652
653
654 if ( typeOut != FLOAT ) {
655 theXY = theXX;
656 theXY += 4*dimxXdimy;
657 }
658
659
660
661 for ( z=0; z<bufferDims[2]; z++ ) {
662
663 switch( typeIn ) {
664 default :
665 break;
666 case UCHAR :
667 case SCHAR :
668 sliceIn = (void*)( ((u8*)bufferIn) + z*dimxXdimy ); break;
669 case USHORT :
670 case SSHORT :
671 sliceIn = (void*)( ((u16*)bufferIn) + z*dimxXdimy ); break;
672 case FLOAT :
673 sliceIn = (void*)( ((float*)bufferIn) + z*dimxXdimy ); break;
674 case DOUBLE :
675 sliceIn = (void*)( ((double*)bufferIn) + z*dimxXdimy ); break;
676 }
677 if ( typeOut == FLOAT ) {
678 theXY = ((float*)bufferOut) + z * dimxXdimy;
679 }
680
681 if ( RecursiveFilterOnBuffer( sliceIn, typeIn, theX, FLOAT,
682 sliceDims, borderLengths,
683 Ysmooth, filterCoefs, filterType ) == 0 ) {
684 if ( _VERBOSE_ > 0 ) {
685 fprintf( stderr, " Fatal error in %s:", proc );
686 fprintf( stderr, " unable to compute Y^0 derivative.\n" );
687 }
688 free( theXX );
689 return( EXIT_ON_FAILURE );
690 }
691
692 if ( RecursiveFilterOnBuffer( sliceIn, typeIn, theY, FLOAT,
693 sliceDims, borderLengths,
694 Xsmooth, filterCoefs, filterType ) == 0 ) {
695 if ( _VERBOSE_ > 0 ) {
696 fprintf( stderr, " Fatal error in %s:", proc );
697 fprintf( stderr, " unable to compute X^0 derivative.\n" );
698 }
699 free( theXX );
700 return( EXIT_ON_FAILURE );
701 }
702
703
704
705
706 if ( RecursiveFilterOnBuffer( sliceIn, typeIn, theXY, FLOAT,
707 sliceDims, borderLengths,
708 XYderiv, filterCoefs, filterType ) == 0 ) {
709 if ( _VERBOSE_ > 0 ) {
710 fprintf( stderr, " Fatal error in %s:", proc );
711 fprintf( stderr, " unable to compute X^1Y^1 derivative.\n" );
712 }
713 free( theXX );
714 return( EXIT_ON_FAILURE );
715 }
716
717
718
719
720 if ( RecursiveFilterOnBuffer( theX, FLOAT, theXX, FLOAT,
721 sliceDims, borderLengths,
722 XXderiv, filterCoefs, filterType ) == 0 ) {
723 if ( _VERBOSE_ > 0 ) {
724 fprintf( stderr, " Fatal error in %s:", proc );
725 fprintf( stderr, " unable to compute X^2 derivative.\n" );
726 }
727 free( theXX );
728 return( EXIT_ON_FAILURE );
729 }
730
731 if ( RecursiveFilterOnBuffer( theY, FLOAT, theYY, FLOAT,
732 sliceDims, borderLengths,
733 YYderiv, filterCoefs, filterType ) == 0 ) {
734 if ( _VERBOSE_ > 0 ) {
735 fprintf( stderr, " Fatal error in %s:", proc );
736 fprintf( stderr, " unable to compute Y^2 derivative.\n" );
737 }
738 free( theXX );
739 return( EXIT_ON_FAILURE );
740 }
741
742
743
744
745 if ( RecursiveFilterOnBuffer( theX, FLOAT, theX, FLOAT,
746 sliceDims, borderLengths,
747 Xderiv, filterCoefs, filterType ) == 0 ) {
748 if ( _VERBOSE_ > 0 ) {
749 fprintf( stderr, " Fatal error in %s:", proc );
750 fprintf( stderr, " unable to compute X^1 derivative.\n" );
751 }
752 free( theXX );
753 return( EXIT_ON_FAILURE );
754 }
755
756 if ( RecursiveFilterOnBuffer( theY, FLOAT, theY, FLOAT,
757 sliceDims, borderLengths,
758 Yderiv, filterCoefs, filterType ) == 0 ) {
759 if ( _VERBOSE_ > 0 ) {
760 fprintf( stderr, " Fatal error in %s:", proc );
761 fprintf( stderr, " unable to compute Y^1 derivative.\n" );
762 }
763 free( theXX );
764 return( EXIT_ON_FAILURE );
765 }
766
767
768
769
770 for ( i=0; i<dimxXdimy; i++ ) {
771 gx = theX[i];
772 gy = theY[i];
773 g = (gx*gx + gy*gy);
774 theXY[i] = (float)(gx * ( theXX[i] * gx + theXY[i] * gy ) +
775 gy * ( theXY[i] * gx + theYY[i] * gy ));
776 if ( g > 1e-10 ) theXY[i] = (float)(theXY[i] / g);
777 }
778
779 if ( typeOut != FLOAT ) {
780 switch ( typeOut ) {
781 case UCHAR :
782 sliceOut = (((u8*)bufferOut) + z * dimxXdimy);
783 break;
784 case SCHAR :
785 sliceOut = (((s8*)bufferOut) + z * dimxXdimy);
786 break;
787 case SSHORT :
788 sliceOut = (((s16*)bufferOut) + z * dimxXdimy);
789 break;
790 case DOUBLE :
791 sliceOut = (((r64*)bufferOut) + z * dimxXdimy);
792 break;
793 default :
794 if ( _VERBOSE_ > 0 )
795 fprintf( stderr, " Error in %s: such output type not handled.\n", proc );
796 free( theXX );
797 return( EXIT_ON_FAILURE );
798 }
799 ConvertBuffer( theXY, FLOAT, sliceOut, typeOut, dimxXdimy );
800 }
801 }
802
803 return( EXIT_ON_SUCCESS );
804 }
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
GradientHessianGradient(void * bufferIn,bufferType typeIn,void * bufferOut,bufferType typeOut,int * bufferDims,int * borderLengths,float * filterCoefs,recursiveFilterType filterType)827 int GradientHessianGradient ( void *bufferIn,
828 bufferType typeIn,
829 void *bufferOut,
830 bufferType typeOut,
831 int *bufferDims,
832 int *borderLengths,
833 float *filterCoefs,
834 recursiveFilterType filterType )
835 {
836 char *proc = "GradientHessianGradient";
837
838
839
840 float *theZZ = NULL;
841 float *theZ = NULL;
842 float *theZ1 = NULL;
843 float *theZ0 = NULL;
844
845 float *theXZ = NULL;
846 float *theYZ = NULL;
847
848 float *theXX = NULL;
849 float *theYY = NULL;
850 float *theXY = NULL;
851
852 float *theX = NULL;
853 float *theY = NULL;
854
855
856 derivativeOrder ZZderiv[3] = { SMOOTHING, SMOOTHING, DERIVATIVE_2 };
857 derivativeOrder Zderiv[3] = { SMOOTHING, SMOOTHING, DERIVATIVE_1 };
858 derivativeOrder Z1deriv[3] = { NODERIVATIVE, NODERIVATIVE, DERIVATIVE_1 };
859 derivativeOrder Z0deriv[3] = { NODERIVATIVE, NODERIVATIVE, SMOOTHING };
860
861 derivativeOrder XZderiv[3] = { DERIVATIVE_1, SMOOTHING, NODERIVATIVE };
862 derivativeOrder YZderiv[3] = { SMOOTHING, DERIVATIVE_1, NODERIVATIVE };
863
864 derivativeOrder XXderiv[3] = { DERIVATIVE_2, SMOOTHING, NODERIVATIVE };
865 derivativeOrder YYderiv[3] = { SMOOTHING, DERIVATIVE_2, NODERIVATIVE };
866 derivativeOrder XYderiv[3] = { DERIVATIVE_1, DERIVATIVE_1, NODERIVATIVE };
867
868 derivativeOrder Xderiv[3] = { DERIVATIVE_1, SMOOTHING, NODERIVATIVE };
869 derivativeOrder Yderiv[3] = { SMOOTHING, DERIVATIVE_1, NODERIVATIVE };
870
871 int sliceDims[3];
872 int z, i, j, dimxXdimy;
873
874 double gx, gy, gz, g;
875
876 /*
877 * We check the buffers' dimensions.
878 */
879 if ( bufferDims[2] == 1 ) {
880 return( GradientHessianGradient_2D ( bufferIn, typeIn, bufferOut, typeOut,
881 bufferDims, borderLengths, filterCoefs, filterType ) );
882 }
883
884 if ( (bufferDims[0] <= 0) || (bufferDims[1] <= 0) || (bufferDims[2] <= 0) ) {
885 if ( _VERBOSE_ > 0 )
886 fprintf( stderr, " Fatal error in %s: improper buffer's dimension.\n", proc );
887 return( EXIT_ON_FAILURE );
888 }
889
890 /*
891 * test of the coefficients
892 */
893 if ( (filterCoefs[0] < 0.0) || (filterCoefs[1] < 0.0) ||
894 (filterCoefs[2] < 0.0) ) {
895 if ( _VERBOSE_ > 0 )
896 fprintf( stderr, " Error in %s: negative coefficient's value.\n", proc );
897 return( EXIT_ON_FAILURE );
898 }
899
900
901 /*
902 *
903 */
904 dimxXdimy = bufferDims[0] * bufferDims[1];
905 sliceDims[0] = bufferDims[0];
906 sliceDims[1] = bufferDims[1];
907 sliceDims[2] = 1;
908
909
910 if ( typeOut == FLOAT ) {
911 theX = (float*)malloc( (7+3*bufferDims[2]) * dimxXdimy * sizeof( float ) );
912 } else {
913 theX = (float*)malloc( (7+4*bufferDims[2]) * dimxXdimy * sizeof( float ) );
914 }
915
916
917 if ( theX == NULL ) {
918 if ( _VERBOSE_ > 0 ) {
919 fprintf( stderr, " Fatal error in %s:", proc );
920 fprintf( stderr, " unable to allocate auxiliary buffer.\n" );
921 }
922 return( EXIT_ON_FAILURE );
923 }
924
925 /*
926 * BUFFERS
927 *
928 * slices : theX theY theXY theYY theXX theYZ theXZ
929 *
930 * volumes : theZ0 theZ1 theZ theZZ
931 *
932 */
933
934 theY = theXX = theXY = theYY = theYZ = theXZ = theX;
935 theZ0 = theZ1 = theZ = theX;
936
937 theY += dimxXdimy;
938 theXY += 2*dimxXdimy;
939 theYY += 3*dimxXdimy;
940 theXX += 4*dimxXdimy;
941 theYZ += 5*dimxXdimy;
942 theXZ += 6*dimxXdimy;
943
944 theZ0 += 7*dimxXdimy;
945 theZ1 += 7*dimxXdimy + bufferDims[2]*dimxXdimy;
946 theZ += 7*dimxXdimy + 2*bufferDims[2]*dimxXdimy;
947
948 if ( typeOut == FLOAT ) {
949 theZZ = (float*)bufferOut;
950 } else {
951 theZZ = theX;
952 theZZ += 7*dimxXdimy + 3*bufferDims[2]*dimxXdimy;
953 }
954
955
956
957 /*
958 *
959 * 3D filtering / filtering along Z
960 *
961 */
962
963 if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theZ0, FLOAT,
964 bufferDims, borderLengths,
965 Z0deriv, filterCoefs, filterType ) == 0 ) {
966 if ( _VERBOSE_ > 0 ) {
967 fprintf( stderr, " Fatal error in %s:", proc );
968 fprintf( stderr, " unable to compute Z^0 derivative.\n" );
969 }
970 free( theX );
971 return( EXIT_ON_FAILURE );
972 }
973
974 if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theZ1, FLOAT,
975 bufferDims, borderLengths,
976 Z1deriv, filterCoefs, filterType ) == 0 ) {
977 if ( _VERBOSE_ > 0 ) {
978 fprintf( stderr, " Fatal error in %s:", proc );
979 fprintf( stderr, " unable to compute Z^1 derivative.\n" );
980 }
981 free( theX );
982 return( EXIT_ON_FAILURE );
983 }
984
985 if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theZ, FLOAT,
986 bufferDims, borderLengths,
987 Zderiv, filterCoefs, filterType ) == 0 ) {
988 if ( _VERBOSE_ > 0 ) {
989 fprintf( stderr, " Fatal error in %s:", proc );
990 fprintf( stderr, " unable to compute Z^1 derivative (edge).\n" );
991 }
992 free( theX );
993 return( EXIT_ON_FAILURE );
994 }
995
996 if ( RecursiveFilterOnBuffer( bufferIn, typeIn, theZZ, FLOAT,
997 bufferDims, borderLengths,
998 ZZderiv, filterCoefs, filterType ) == 0 ) {
999 if ( _VERBOSE_ > 0 ) {
1000 fprintf( stderr, " Fatal error in %s:", proc );
1001 fprintf( stderr, " unable to compute Z^2 derivative.\n" );
1002 }
1003 free( theX );
1004 return( EXIT_ON_FAILURE );
1005 }
1006
1007
1008 /*
1009 * theZ0 : smoothed along Z
1010 * theZ1 : first derivative along Z
1011 * theZ : first derivative along Z, smoothed along X and Y
1012 * theZZ : second derivative along Z, smoothed along X and Y
1013 */
1014
1015
1016
1017 for ( z=0; z<bufferDims[2]; z++ ) {
1018 fprintf( stderr, "%s: processing slice %3d/%d\r",
1019 proc, z, bufferDims[2] );
1020 /*
1021 *
1022 * 2D filtering / filtering along X and Y
1023 *
1024 */
1025
1026 if ( RecursiveFilterOnBuffer( theZ1+z*dimxXdimy, FLOAT, theXZ, FLOAT,
1027 sliceDims, borderLengths,
1028 XZderiv, filterCoefs, filterType ) == 0 ) {
1029 if ( _VERBOSE_ > 0 ) {
1030 fprintf( stderr, " Fatal error in %s:", proc );
1031 fprintf( stderr, " unable to compute X^1Z^1 derivative.\n" );
1032 }
1033 free( theX );
1034 return( EXIT_ON_FAILURE );
1035 }
1036
1037 if ( RecursiveFilterOnBuffer( theZ1+z*dimxXdimy, FLOAT, theYZ, FLOAT,
1038 sliceDims, borderLengths,
1039 YZderiv, filterCoefs, filterType ) == 0 ) {
1040 if ( _VERBOSE_ > 0 ) {
1041 fprintf( stderr, " Fatal error in %s:", proc );
1042 fprintf( stderr, " unable to compute Y^1Z^1 derivative.\n" );
1043 }
1044 free( theX );
1045 return( EXIT_ON_FAILURE );
1046 }
1047
1048
1049
1050
1051 if ( RecursiveFilterOnBuffer( theZ0+z*dimxXdimy, FLOAT, theXX, FLOAT,
1052 sliceDims, borderLengths,
1053 XXderiv, filterCoefs, filterType ) == 0 ) {
1054 if ( _VERBOSE_ > 0 ) {
1055 fprintf( stderr, " Fatal error in %s:", proc );
1056 fprintf( stderr, " unable to compute X^2 derivative.\n" );
1057 }
1058 free( theX );
1059 return( EXIT_ON_FAILURE );
1060 }
1061
1062 if ( RecursiveFilterOnBuffer( theZ0+z*dimxXdimy, FLOAT, theYY, FLOAT,
1063 sliceDims, borderLengths,
1064 YYderiv, filterCoefs, filterType ) == 0 ) {
1065 if ( _VERBOSE_ > 0 ) {
1066 fprintf( stderr, " Fatal error in %s:", proc );
1067 fprintf( stderr, " unable to compute Y^2 derivative.\n" );
1068 }
1069 free( theX );
1070 return( EXIT_ON_FAILURE );
1071 }
1072
1073 if ( RecursiveFilterOnBuffer( theZ0+z*dimxXdimy, FLOAT, theXY, FLOAT,
1074 sliceDims, borderLengths,
1075 XYderiv, filterCoefs, filterType ) == 0 ) {
1076 if ( _VERBOSE_ > 0 ) {
1077 fprintf( stderr, " Fatal error in %s:", proc );
1078 fprintf( stderr, " unable to compute X^1Y^1 derivative.\n" );
1079 }
1080 free( theX );
1081 return( EXIT_ON_FAILURE );
1082 }
1083
1084
1085
1086 if ( RecursiveFilterOnBuffer( theZ0+z*dimxXdimy, FLOAT, theX, FLOAT,
1087 sliceDims, borderLengths,
1088 Xderiv, filterCoefs, filterType ) == 0 ) {
1089 if ( _VERBOSE_ > 0 ) {
1090 fprintf( stderr, " Fatal error in %s:", proc );
1091 fprintf( stderr, " unable to compute X^1 derivative (edge).\n" );
1092 }
1093 free( theX );
1094 return( EXIT_ON_FAILURE );
1095 }
1096
1097 if ( RecursiveFilterOnBuffer( theZ0+z*dimxXdimy, FLOAT, theY, FLOAT,
1098 sliceDims, borderLengths,
1099 Yderiv, filterCoefs, filterType ) == 0 ) {
1100 if ( _VERBOSE_ > 0 ) {
1101 fprintf( stderr, " Fatal error in %s:", proc );
1102 fprintf( stderr, " unable to compute Y^1 derivative (edge).\n" );
1103 }
1104 free( theX );
1105 return( EXIT_ON_FAILURE );
1106 }
1107
1108
1109
1110 for ( j=z*dimxXdimy, i=0; i<dimxXdimy; j++, i++ ) {
1111 gx = theX[i];
1112 gy = theY[i];
1113 gz = theZ[j];
1114 g = gx*gx + gy*gy + gz*gz;
1115 theZZ[j] = (float)(gx * ( theXX[i] * gx + theXY[i] * gy + theXZ[i] * gz ) +
1116 gy * ( theXY[i] * gx + theYY[i] * gy + theYZ[i] * gz ) +
1117 gz * ( theXZ[i] * gx + theYZ[i] * gy + theZZ[j] * gz ));
1118 if ( g > 1e-10 ) theZZ[j] = (float)(theZZ[j] / g);
1119
1120 }
1121
1122 }
1123
1124 if ( typeOut != FLOAT ) {
1125 ConvertBuffer( theZZ, FLOAT, bufferOut, typeOut, bufferDims[2]*dimxXdimy );
1126 }
1127
1128 free( theX );
1129
1130 return( EXIT_ON_SUCCESS );
1131 }
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169 /*
1170 *
1171 *
1172 *
1173 *
1174 */
1175
RecursiveFilterOnBuffer(void * bufferIn,bufferType typeIn,void * bufferOut,bufferType typeOut,int * bufferDims,int * borderLengths,derivativeOrder * derivatives,float * filterCoefs,recursiveFilterType filterType)1176 int RecursiveFilterOnBuffer( void *bufferIn,
1177 bufferType typeIn,
1178 void *bufferOut,
1179 bufferType typeOut,
1180 int *bufferDims,
1181 int *borderLengths,
1182 derivativeOrder *derivatives,
1183 float *filterCoefs,
1184 recursiveFilterType filterType )
1185 {
1186 char *proc = "RecursiveFilterOnBuffer";
1187 register int dimx, dimxXdimy;
1188 int dimy, dimz;
1189 register int x, y, z;
1190 /*
1191 *obviously, we need to perform the computation
1192 * with float or double values. For this reason,
1193 * we allocate an auxiliary buffer if the output buffer
1194 * is not of type float or double.
1195 */
1196 void *bufferToBeProcessed = (void*)NULL;
1197 bufferType typeToBeProcessed = TYPE_UNKNOWN;
1198 void *bufferResult = (void*)NULL;
1199 bufferType typeResult = TYPE_UNKNOWN;
1200 /*
1201 * lines' lengths
1202 */
1203 int lengthX = 0;
1204 int lengthY = 0;
1205 int lengthZ = 0;
1206 int maxLengthline = 0;
1207 int borderXlength = 0;
1208 int borderYlength = 0;
1209 int borderZlength = 0;
1210 /*
1211 * 1D arrays for computations.
1212 */
1213 double *theLine = (double*)NULL;
1214 double *resLine = (double*)NULL;
1215 double *tmpLine = (double*)NULL;
1216 /*
1217 * pointers for computations;
1218 */
1219 register r32 *r32firstPoint = (r32*)NULL;
1220 register r64 *r64firstPoint = (r64*)NULL;
1221 register r32 *r32_pt = (r32*)NULL;
1222 register r64 *r64_pt = (r64*)NULL;
1223 register double *dbl_pt1 = (double*)NULL;
1224 register double *dbl_pt2 = (double*)NULL;
1225 register double dbl_first = 0.0;
1226 register double dbl_last = 0.0;
1227 int offsetLastPoint = 0;
1228 int offsetNextFirstPoint = 0;
1229 register r32 *r32firstPointResult = (r32*)NULL;
1230 register r64 *r64firstPointResult = (r64*)NULL;
1231 double *theLinePlusBorder = (double*)NULL;
1232 double *resLinePlusBorder = (double*)NULL;
1233
1234 RFcoefficientType *RFC = NULL;
1235
1236 /*
1237 * We check the buffers' dimensions.
1238 */
1239 dimx = bufferDims[0]; dimy = bufferDims[1]; dimz = bufferDims[2];
1240 dimxXdimy = dimx * dimy;
1241 if ( (dimx <= 0) || (dimy <= 0) || (dimz <= 0) ) {
1242 if ( _VERBOSE_ > 0 )
1243 fprintf( stderr, " Fatal error in %s: improper buffer's dimension.\n", proc );
1244 return( EXIT_ON_FAILURE );
1245 }
1246 /*
1247 * We check the pointers.
1248 */
1249 if ( (bufferIn == (void*)NULL) || (bufferOut == (void*)NULL) ) {
1250 if ( _VERBOSE_ > 0 )
1251 fprintf( stderr, " Fatal error in %s: NULL pointer on buffer.\n", proc );
1252 return( EXIT_ON_FAILURE );
1253 }
1254
1255 /*
1256 * May we use the buffer bufferOut as the bufferResult?
1257 * If its type is FLOAT or DOUBLE, then yes.
1258 * If not, we have to allocate an auxiliary buffer.
1259 */
1260 if ( (typeOut == FLOAT) || (typeOut == DOUBLE) ) {
1261 bufferResult = bufferOut;
1262 typeResult = typeOut;
1263 } else {
1264 bufferResult = (void*)malloc( (dimx*dimy*dimz) * sizeof(r32) );
1265 if ( bufferResult == (void*)NULL ) {
1266 if ( _VERBOSE_ > 0 )
1267 fprintf( stderr, " Fatal error in %s: unable to allocate auxiliary buffer.\n", proc );
1268 return( EXIT_ON_FAILURE );
1269 }
1270 typeResult = FLOAT;
1271 }
1272
1273 /*
1274 * May we consider the buffer bufferIn as the bufferToBeProcessed?
1275 * If its type is FLOAT or DOUBLE, then yes.
1276 * If not, we convert it into the buffer bufferResult, and this
1277 * last buffer is the bufferToBeProcessed.
1278 */
1279 if ( (typeIn == FLOAT) || (typeIn == DOUBLE) ) {
1280 bufferToBeProcessed = bufferIn;
1281 typeToBeProcessed = typeIn;
1282 } else {
1283 ConvertBuffer( bufferIn, typeIn, bufferResult, typeResult, (dimx*dimy*dimz) );
1284 bufferToBeProcessed = bufferResult;
1285 typeToBeProcessed = typeResult;
1286 }
1287
1288 /*
1289 * Estimation of the lines' length along each direction.
1290 */
1291 if ( borderLengths != NULL ) {
1292 borderXlength = borderLengths[0];
1293 borderYlength = borderLengths[1];
1294 borderZlength = borderLengths[2];
1295 if ( borderXlength < 0 ) borderXlength = 0;
1296 if ( borderYlength < 0 ) borderYlength = 0;
1297 if ( borderZlength < 0 ) borderZlength = 0;
1298 }
1299
1300 /*
1301 * Tue Jul 6 19:15:15 MET DST 1999 (gregoire.malandainoire Malandain)
1302 * changes 3 x dimx -> dimx, dimy, dimz
1303 */
1304 lengthX = dimx + 2 * borderXlength;
1305 lengthY = dimy + 2 * borderYlength;
1306 lengthZ = dimz + 2 * borderZlength;
1307 maxLengthline = lengthX;
1308 if ( maxLengthline < lengthY ) maxLengthline = lengthY;
1309 if ( maxLengthline < lengthZ ) maxLengthline = lengthZ;
1310 if ( maxLengthline <= 0 ) {
1311 if ( _VERBOSE_ > 0 )
1312 fprintf( stderr, " Error in %s: unable to deal with dimensions = 0.\n", proc );
1313 if ( (typeOut != FLOAT) && (typeOut != DOUBLE) )
1314 free( bufferResult );
1315 return( EXIT_ON_FAILURE );
1316 }
1317 /*
1318 * Allocations of work arrays.
1319 * We will use them to process each line.
1320 */
1321 theLine = (double*)malloc( 3 * maxLengthline * sizeof(double) );
1322 if ( theLine == (double*)NULL ) {
1323 if ( _VERBOSE_ > 0 )
1324 fprintf( stderr, " Fatal error in %s: unable to allocate auxiliary work arrays.\n", proc );
1325 if ( (typeOut != FLOAT) && (typeOut != DOUBLE) )
1326 free( bufferResult );
1327 return( EXIT_ON_FAILURE );
1328 }
1329 resLine = theLine + maxLengthline;
1330 tmpLine = resLine + maxLengthline;
1331
1332 /*
1333 * From now,
1334 * typeToBeProcessed is either FLOAT or DOUBLE
1335 * so is typeResult.
1336 */
1337
1338
1339 /*
1340 * Processing along X.
1341 */
1342 if ( dimx > 4 )
1343 if (derivatives[0] != NODERIVATIVE)
1344 if (filterCoefs[0] > 0.0) {
1345
1346 if ( _VERBOSE_ != 0 )
1347 fprintf( stderr, " %s: processing along X.\n", proc );
1348
1349 RFC = InitRecursiveCoefficients( (double)filterCoefs[0], filterType, derivatives[0] );
1350
1351 if ( RFC == NULL ) {
1352 if ( _VERBOSE_ != 0 )
1353 fprintf( stderr, " %s: unable to allocate coefficients\n", proc );
1354 if ( (typeOut != FLOAT) && (typeOut != DOUBLE) )
1355 free( bufferResult );
1356 return( EXIT_ON_FAILURE );
1357 }
1358
1359 r64firstPoint = (r64*)bufferToBeProcessed;
1360 r32firstPoint = (r32*)bufferToBeProcessed;
1361
1362 r64firstPointResult = (r64*)bufferResult;
1363 r32firstPointResult = (r32*)bufferResult;
1364
1365 offsetLastPoint = borderXlength + dimx - 1;
1366
1367 theLinePlusBorder = theLine + borderXlength;
1368 resLinePlusBorder = resLine + borderXlength;
1369
1370 /*
1371 * There are dimz*dimy X lines to be processed.
1372 */
1373 for ( z=0; z<dimz; z++ )
1374 for ( y=0; y<dimy; y++ ) {
1375 /*
1376 * Acquiring a X line.
1377 */
1378 dbl_pt1 = theLinePlusBorder;
1379 switch ( typeToBeProcessed ) {
1380 case DOUBLE :
1381 (void)memcpy( (void*)dbl_pt1, (void*)r64firstPoint, dimx * sizeof(r64) );
1382 r64firstPoint += dimx;
1383 break;
1384 case FLOAT :
1385 default :
1386 for ( x=0; x<dimx; x++, dbl_pt1++, r32firstPoint++ ) *dbl_pt1 = *r32firstPoint;
1387 }
1388 /*
1389 * Adding points at both ends of the line.
1390 */
1391 if ( borderXlength > 0 ) {
1392 dbl_pt1 = theLine + borderXlength; dbl_first = *dbl_pt1;
1393 dbl_pt2 = theLine + offsetLastPoint; dbl_last = *dbl_pt2;
1394 for ( x=0; x<borderXlength; x++ ) {
1395 *--dbl_pt1 = dbl_first;
1396 *++dbl_pt2 = dbl_last;
1397 }
1398 }
1399 /*
1400 * Processing the line.
1401 */
1402 if ( RecursiveFilter1D( RFC, theLine, resLine, tmpLine, resLine, lengthX ) == 0 ) {
1403 if ( _VERBOSE_ != 0 )
1404 fprintf(stderr," Error in %s: unable to process X line (y=%d,z=%d).\n", proc, y, z);
1405 if ( (typeOut != FLOAT) && (typeOut != DOUBLE) )
1406 free( bufferResult );
1407 free( (void*)theLine );
1408 return( EXIT_ON_FAILURE );
1409 }
1410 /*
1411 * Copy the result into the buffer bufferResult.
1412 */
1413 dbl_pt1 = resLinePlusBorder;
1414 switch ( typeResult ) {
1415 case DOUBLE :
1416 (void)memcpy( (void*)r64firstPointResult, (void*)dbl_pt1, dimx * sizeof(r64) );
1417 r64firstPointResult += dimx;
1418 break;
1419 case FLOAT :
1420 default :
1421 for ( x=0; x<dimx; x++, dbl_pt1++, r32firstPointResult++ )
1422 *r32firstPointResult = (r32)(*dbl_pt1);
1423 }
1424 }
1425
1426 /*
1427 * The next buffer to be processed is the buffer
1428 * bufferResult.
1429 */
1430 bufferToBeProcessed = bufferResult;
1431 typeToBeProcessed = typeResult;
1432
1433 free( RFC );
1434 RFC = NULL;
1435
1436 } /* end of Processing along X. */
1437
1438 /*
1439 * Processing along Y.
1440 */
1441 if ( dimy > 4 )
1442 if (derivatives[1] != NODERIVATIVE)
1443 if (filterCoefs[1] > 0.0) {
1444
1445 if ( _VERBOSE_ != 0 )
1446 fprintf( stderr, " %s: processing along Y.\n", proc );
1447
1448 RFC = InitRecursiveCoefficients( (double)filterCoefs[1], filterType, derivatives[1] );
1449
1450 if ( RFC == NULL ) {
1451 if ( _VERBOSE_ != 0 )
1452 fprintf( stderr, " %s: unable to allocate coefficients\n", proc );
1453 if ( (typeOut != FLOAT) && (typeOut != DOUBLE) )
1454 free( bufferResult );
1455 return( EXIT_ON_FAILURE );
1456 }
1457
1458 r64firstPoint = (r64*)bufferToBeProcessed;
1459 r32firstPoint = (r32*)bufferToBeProcessed;
1460
1461 r64firstPointResult = (r64*)bufferResult;
1462 r32firstPointResult = (r32*)bufferResult;
1463
1464 offsetLastPoint = borderYlength + dimy - 1;
1465 offsetNextFirstPoint = dimx * dimy - dimx;
1466
1467 theLinePlusBorder = theLine + borderYlength;
1468 resLinePlusBorder = resLine + borderYlength;
1469
1470 /*
1471 * There are dimz*dimx Y lines to be processed.
1472 */
1473 for ( z=0; z<dimz; z++ ) {
1474 for ( x=0; x<dimx; x++ ) {
1475 /*
1476 * Acquiring a Y line.
1477 */
1478 dbl_pt1 = theLinePlusBorder;
1479 switch ( typeToBeProcessed ) {
1480 case DOUBLE :
1481 r64_pt = r64firstPoint;
1482 for ( y=0; y<dimy; y++, dbl_pt1++, r64_pt += dimx ) *dbl_pt1 = *r64_pt;
1483 /*
1484 * Going to the first point of the next Y line
1485 */
1486 r64firstPoint ++;
1487 break;
1488 case FLOAT :
1489 default :
1490 r32_pt = r32firstPoint;
1491 for ( y=0; y<dimy; y++, dbl_pt1++, r32_pt += dimx ) *dbl_pt1 = *r32_pt;
1492 r32firstPoint ++;
1493 }
1494 /*
1495 * Adding points at both ends of the line.
1496 */
1497 if ( borderYlength > 0 ) {
1498 dbl_pt1 = theLine + borderYlength; dbl_first = *dbl_pt1;
1499 dbl_pt2 = theLine + offsetLastPoint; dbl_last = *dbl_pt2;
1500 for ( y=0; y<borderYlength; y++ ) {
1501 *--dbl_pt1 = dbl_first;
1502 *++dbl_pt2 = dbl_last;
1503 }
1504 }
1505 /*
1506 * Processing the line.
1507 */
1508 if ( RecursiveFilter1D( RFC, theLine, resLine, tmpLine, resLine, lengthY ) == 0 ) {
1509 if ( _VERBOSE_ != 0 )
1510 fprintf(stderr," Error in %s: unable to process Y line (x=%d,z=%d).\n", proc, x, z);
1511 if ( (typeOut != FLOAT) && (typeOut != DOUBLE) )
1512 free( bufferResult );
1513 free( (void*)theLine );
1514 return( EXIT_ON_FAILURE );
1515 }
1516 /*
1517 * Copy the result into the buffer bufferResult.
1518 */
1519 dbl_pt1 = resLinePlusBorder;
1520 switch ( typeResult ) {
1521 case DOUBLE :
1522 r64_pt = r64firstPointResult;
1523 for ( y=0; y<dimy; y++, dbl_pt1++, r64_pt += dimx ) *r64_pt = *dbl_pt1;
1524 r64firstPointResult ++;
1525 break;
1526 case FLOAT :
1527 default :
1528 r32_pt = r32firstPointResult;
1529 for ( y=0; y<dimy; y++, dbl_pt1++, r32_pt += dimx )
1530 *r32_pt = (float)*dbl_pt1;
1531 r32firstPointResult ++;
1532 }
1533 }
1534 /*
1535 * Going to the first point of the next Y line
1536 * which is the first Y line of the next slice.
1537 *
1538 * The pointer r[32,64]firstPoint[Result] has
1539 * already been increased by dimx. To reach
1540 * the first point of the next slice, we
1541 * have to increase it by (dimx*dimy)-dimx.
1542 */
1543 switch ( typeToBeProcessed ) {
1544 case DOUBLE :
1545 r64firstPoint += offsetNextFirstPoint;
1546 break;
1547 case FLOAT :
1548 default :
1549 r32firstPoint += offsetNextFirstPoint;
1550 }
1551 switch ( typeResult ) {
1552 case DOUBLE :
1553 r64firstPointResult += offsetNextFirstPoint;
1554 break;
1555 case FLOAT :
1556 default :
1557 r32firstPointResult += offsetNextFirstPoint;
1558 }
1559 }
1560
1561 /*
1562 * The next buffer to be processed is the buffer
1563 * bufferResult.
1564 */
1565 bufferToBeProcessed = bufferResult;
1566 typeToBeProcessed = typeResult;
1567
1568 free( RFC );
1569 RFC = NULL;
1570
1571 } /* end of Processing along Y. */
1572
1573
1574 /*
1575 * Processing along Z.
1576 */
1577 if ( dimz > 4 )
1578 if (derivatives[2] != NODERIVATIVE)
1579 if (filterCoefs[2] > 0.0) {
1580
1581 if ( _VERBOSE_ != 0 )
1582 fprintf( stderr, " %s: processing along Z.\n", proc );
1583
1584 RFC = InitRecursiveCoefficients( (double)filterCoefs[2], filterType, derivatives[2] );
1585
1586 if ( RFC == NULL ) {
1587 if ( _VERBOSE_ != 0 )
1588 fprintf( stderr, " %s: unable to allocate coefficients\n", proc );
1589 if ( (typeOut != FLOAT) && (typeOut != DOUBLE) )
1590 free( bufferResult );
1591 return( EXIT_ON_FAILURE );
1592 }
1593
1594 r64firstPoint = (r64*)bufferToBeProcessed;
1595 r32firstPoint = (r32*)bufferToBeProcessed;
1596
1597 offsetLastPoint = borderZlength + dimz - 1;
1598
1599 r64firstPointResult = (r64*)bufferResult;
1600 r32firstPointResult = (r32*)bufferResult;
1601
1602 offsetLastPoint = borderZlength + dimz - 1;
1603
1604 theLinePlusBorder = theLine + borderYlength;
1605 resLinePlusBorder = resLine + borderYlength;
1606
1607 /*
1608 * There are dimy*dimx Z lines to be processed.
1609 */
1610 for ( y=0; y<dimy; y++ )
1611 for ( x=0; x<dimx; x++ ) {
1612 /*
1613 * Acquiring a Z line.
1614 */
1615 dbl_pt1 = theLinePlusBorder;
1616 switch ( typeToBeProcessed ) {
1617 case DOUBLE :
1618 r64_pt = r64firstPoint;
1619 for ( z=0; z<dimz; z++, dbl_pt1++, r64_pt += dimxXdimy ) *dbl_pt1 = *r64_pt;
1620 /*
1621 * Going to the first point of the next Z line
1622 */
1623 r64firstPoint ++;
1624 break;
1625 case FLOAT :
1626 default :
1627 r32_pt = r32firstPoint;
1628 for ( z=0; z<dimz; z++, dbl_pt1++, r32_pt += dimxXdimy ) *dbl_pt1 = *r32_pt;
1629 r32firstPoint ++;
1630 }
1631 /*
1632 * Adding points at both ends of the line.
1633 */
1634 if ( borderZlength > 0 ) {
1635 dbl_pt1 = theLine + borderZlength; dbl_first = *dbl_pt1;
1636 dbl_pt2 = theLine + offsetLastPoint; dbl_last = *dbl_pt2;
1637 for ( z=0; z<borderZlength; z++ ) {
1638 *--dbl_pt1 = dbl_first;
1639 *++dbl_pt2 = dbl_last;
1640 }
1641 }
1642 /*
1643 * Processing the line.
1644 */
1645 if ( RecursiveFilter1D( RFC, theLine, resLine, tmpLine, resLine, lengthZ ) == 0 ) {
1646 if ( _VERBOSE_ != 0 )
1647 fprintf(stderr," Error in %s: unable to process Z line (x=%d,y=%d).\n", proc, x, y);
1648 if ( (typeOut != FLOAT) && (typeOut != DOUBLE) )
1649 free( bufferResult );
1650 free( (void*)theLine );
1651 return( EXIT_ON_FAILURE );
1652 }
1653
1654 /*
1655 * Copy the result into the buffer bufferResult.
1656 */
1657 dbl_pt1 = resLinePlusBorder;
1658 switch ( typeResult ) {
1659 case DOUBLE :
1660 r64_pt = r64firstPointResult;
1661 for ( z=0; z<dimz; z++, dbl_pt1++, r64_pt += dimxXdimy )
1662 *r64_pt = *dbl_pt1;
1663 r64firstPointResult ++;
1664 break;
1665 case FLOAT :
1666 default :
1667 r32_pt = r32firstPointResult;
1668 for ( z=0; z<dimz; z++, dbl_pt1++, r32_pt += dimxXdimy )
1669 *r32_pt = (float)*dbl_pt1;
1670 r32firstPointResult ++;
1671 }
1672 }
1673
1674 free( RFC );
1675 RFC = NULL;
1676
1677 } /* end of Processing along Z. */
1678
1679
1680
1681
1682 /*
1683 * From bufferResult to bufferOut
1684 */
1685 ConvertBuffer( bufferResult, typeResult, bufferOut, typeOut, (dimx*dimy*dimz) );
1686
1687 /*
1688 * Releasing the buffers.
1689 */
1690 if ( (typeOut != FLOAT) && (typeOut != DOUBLE) )
1691 free( bufferResult );
1692 free( (void*)theLine );
1693
1694 return( EXIT_ON_SUCCESS );
1695 }
1696
1697
1698
1699
1700
1701
Recbuffer_verbose()1702 void Recbuffer_verbose ( )
1703 {
1704 _VERBOSE_ = 1;
1705 Recline_verbose ( );
1706 }
1707
1708
1709
1710
1711
Recbuffer_noverbose()1712 void Recbuffer_noverbose ( )
1713 {
1714 _VERBOSE_ = 0;
1715 Recline_noverbose ( );
1716 }
1717