1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8
9 /*** NOT 7D SAFE ***/
10
11 /* Collect stats on sequence of images.
12 Usage: call repeatedly, once for each image; will return void.
13 after last image, call with NULL as argument, will return
14 pointer to array of two pointers to (float) images, the
15 first the mean, the second the stdev.
16 */
17
mri_stat_seq(MRI_IMAGE * imin)18 MRI_IMAGE ** mri_stat_seq( MRI_IMAGE * imin )
19 {
20
21 /* static variables (exist between calls) */
22
23 static MRI_IMAGE * imsum , * imsumq ;
24 static int nim = 0 , npix , nx , ny ;
25 static double * dbsum , * dbsumq ;
26
27 /* local variables */
28
29 register int ii ;
30 MRI_IMAGE * imfl , * imsd , ** retval ;
31 float * flar , * sdar ;
32 double scl , vscl , vvv ;
33
34 /*** case: set up new problem ***/
35
36 if( nim == 0 ){
37
38 if( imin == NULL ){
39 fprintf(stderr,"mri_stat_seq: NULL argument for initial call!\n") ;
40 EXIT(1) ;
41 }
42
43 if( ! MRI_IS_2D(imin) ){
44 fprintf(stderr,"\n*** mri_stat_seq: only works on 2D images!\n") ;
45 EXIT(1) ;
46 }
47
48 nx = imin->nx ;
49 ny = imin->ny ;
50 npix = nx * ny ;
51
52 imsum = mri_new( nx , ny , MRI_double ) ;
53 imsumq = mri_new( nx , ny , MRI_double ) ;
54 dbsum = mri_data_pointer( imsum ) ;
55 dbsumq = mri_data_pointer( imsumq ) ;
56
57 MRI_COPY_AUX(imsum,imin) ;
58
59 for( ii=0 ; ii < npix ; ii++ ) dbsum[ii] = dbsumq[ii] = 0.0 ;
60
61 }
62
63 /*** case: add new image into problem ***/
64
65 if( imin != NULL ){
66
67 if( imin->nx != nx || imin->ny != ny ){
68 fprintf(stderr,"mri_stat_seq: input image size mismatch!\n") ;
69 EXIT(1) ;
70 }
71
72 if( ! MRI_IS_2D(imin) ){
73 fprintf(stderr,"\n*** mri_stat_seq: only works on 2D images!\n") ;
74 EXIT(1) ;
75 }
76
77 if( imin->kind == MRI_float ){
78 imfl = imin ;
79 } else {
80 imfl = mri_to_float( imin ) ;
81 }
82 flar = mri_data_pointer( imfl ) ;
83
84 for( ii=0 ; ii < npix ; ii++ ){
85 dbsum[ii] += flar[ii] ;
86 dbsumq[ii] += flar[ii] * flar[ii] ;
87 }
88
89 if( imin != imfl ) mri_free( imfl ) ;
90 nim++ ;
91 return NULL ;
92 }
93
94 /*** case: make report on results, and reset static data ***/
95
96 if( nim < 1 ){
97 fprintf(stderr,"mri_stat_seq: # images input=%d; too small!\n",nim) ;
98 EXIT(1) ;
99 }
100
101 imfl = mri_new( nx , ny , MRI_float ) ;
102 imsd = mri_new( nx , ny , MRI_float ) ;
103 flar = mri_data_pointer( imfl ) ;
104 sdar = mri_data_pointer( imsd ) ;
105
106 MRI_COPY_AUX(imfl,imsum) ;
107
108 scl = 1.0 / nim ;
109 vscl = (nim==1) ? 1.0 : sqrt( ((double)(nim))/(nim-1.0) ) ;
110
111 for( ii=0 ; ii < npix ; ii++ ){
112 flar[ii] = scl * dbsum[ii] ;
113 vvv = scl * dbsumq[ii] - flar[ii] * flar[ii] ;
114 sdar[ii] = (vvv > 0.0) ? (sqrt(vvv)*vscl) : 0.0 ;
115 }
116
117 mri_free( imsum ) ; mri_free( imsumq ) ;
118 nim = 0 ;
119
120 retval = (MRI_IMAGE **) malloc( 2 * sizeof(MRI_IMAGE *) ) ;
121 if( retval == NULL ){
122 fprintf(stderr,"mri_stat_seq: malloc error for retval!\n") ;
123 EXIT(1) ;
124 }
125
126 retval[0] = imfl ;
127 retval[1] = imsd ;
128
129 return retval ;
130 }
131