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 /*********************************************************************
8 These are the template routines to form functional images
9 recursively from various data types.
10
11 To create actual routines for this purpose, you must compile
12 the file with the preprocessor symbol DTYPE set to one of
13 the following types:
14
15 byte short int float
16
17 cc -c -DDTYPE=short afni_pcor.c
18 mv -f afni_pcor.o afni_pcor_short.o
19
20 In the example above, the resulting routine will be named
21 PCOR_update_short and will take as input a "short *"
22 (plus the other stuff, which isn't DTYPE dependent).
23 **********************************************************************/
24
25 #ifndef DTYPE
26 #error "Cannot compile, since DTYPE is undefined."
27 #endif
28
29 #undef MAIN
30 #include "afni_pcor.h"
31
32 /** macros for function names defined in this file **/
33
34 #define PCOR_UPDATE TWO_TWO(PCOR_update_,DTYPE)
35
36 /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
37
38 /*** update all voxels with a new array of data
39 inputs: vdata = array[nvox] of new data for each voxel
40 ref = pointer to references structure to use
41 vc = pointer to correlation data structure
42 output: updated vc
43 ***/
44
PCOR_UPDATE(DTYPE * vdata,PCOR_references * ref,PCOR_voxel_corr * vc)45 void PCOR_UPDATE( DTYPE * vdata , PCOR_references * ref , PCOR_voxel_corr * vc )
46 {
47 int vox , jj , /* loop indices */
48 nv = vc->nvox , /* number of voxels */
49 nr = vc->nref ; /* number of references */
50
51 float *aaa = ref->alp ,
52 *fff = ref->ff ,
53 *ggg = ref->gg ;
54
55 float zz , bq = ref->betasq ;
56
57 /*** check inputs for OK-ness ***/
58
59 if( vc->nref != ref->nref ){
60 fprintf( stderr , "PCOR_UPDATE: reference size mismatch!\n" ) ;
61 EXIT(1) ;
62 }
63
64 /**----------------------------------------------------------------------
65 innermost loop expansion is for speedup if nref is small, if enabled
66 UPZZ updates zz for each row element > This pair is performed for
67 UPCH updates the Cholesky row element > each non-diagonal element
68 UPLL updates the Cholesky diagonal element
69 -------------------------------------------------------------------------**/
70
71 #ifdef EXPAND_UPDATE
72 # define UPZZ(j) zz -= aaa[j] * VCH(vc,vox,j)
73 # define UPCH(j) VCH(vc,vox,j) = fff[j] * VCH(vc,vox,j) + ggg[j] * zz
74 # define UPLL(j) VCH(vc,vox,j) += bq * zz * zz
75
76 switch( nr ){
77 default: /*** generic case: nested loops ***/
78 #endif
79
80 /*** for each voxel ***/
81
82 for( vox=0 ; vox < nv ; vox++ ){
83
84 /*** update last row of each Cholesky factor ***/
85
86 zz = (float) vdata[vox] ;
87 for( jj=0 ; jj < nr ; jj++ ){
88 zz -= aaa[jj] * VCH(vc,vox,jj) ;
89 VCH(vc,vox,jj) = fff[jj] * VCH(vc,vox,jj) + ggg[jj] * zz ;
90 }
91 VCH(vc,vox,nr) += bq * zz * zz ; /* square of true Cholesky diagonal */
92 }
93
94 #ifdef EXPAND_UPDATE
95 break ;
96
97 /***------------------------------------------------------------
98 Below here are the special cases for 1-9 reference waveforms:
99 The above loop over jj is just unrolled manually
100 (using the UP?? macros) in order to gain speed.
101 ----------------------------------------------------------------***/
102
103 case 1:
104 for( vox=0 ; vox < nv ; vox++ ){
105 zz = (float) vdata[vox] ;
106 UPZZ(0) ; UPCH(0) ; UPLL(1) ;
107 }
108 break ;
109
110 case 2:
111 for( vox=0 ; vox < nv ; vox++ ){
112 zz = (float) vdata[vox] ;
113 UPZZ(0) ; UPCH(0) ;
114 UPZZ(1) ; UPCH(1) ; UPLL(2) ;
115 }
116 break ;
117
118 case 3:
119 for( vox=0 ; vox < nv ; vox++ ){
120 zz = (float) vdata[vox] ;
121 UPZZ(0) ; UPCH(0) ;
122 UPZZ(1) ; UPCH(1) ;
123 UPZZ(2) ; UPCH(2) ; UPLL(3) ;
124 }
125 break ;
126
127 case 4:
128 for( vox=0 ; vox < nv ; vox++ ){
129 zz = (float) vdata[vox] ;
130 UPZZ(0) ; UPCH(0) ;
131 UPZZ(1) ; UPCH(1) ;
132 UPZZ(2) ; UPCH(2) ;
133 UPZZ(3) ; UPCH(3) ; UPLL(4) ;
134 }
135 break ;
136
137 case 5:
138 for( vox=0 ; vox < nv ; vox++ ){
139 zz = (float) vdata[vox] ;
140 UPZZ(0) ; UPCH(0) ;
141 UPZZ(1) ; UPCH(1) ;
142 UPZZ(2) ; UPCH(2) ;
143 UPZZ(3) ; UPCH(3) ;
144 UPZZ(4) ; UPCH(4) ; UPLL(5) ;
145 }
146 break ;
147
148 case 6:
149 for( vox=0 ; vox < nv ; vox++ ){
150 zz = (float) vdata[vox] ;
151 UPZZ(0) ; UPCH(0) ;
152 UPZZ(1) ; UPCH(1) ;
153 UPZZ(2) ; UPCH(2) ;
154 UPZZ(3) ; UPCH(3) ;
155 UPZZ(4) ; UPCH(4) ;
156 UPZZ(5) ; UPCH(5) ; UPLL(6) ;
157 }
158 break ;
159
160 case 7:
161 for( vox=0 ; vox < nv ; vox++ ){
162 zz = (float) vdata[vox] ;
163 UPZZ(0) ; UPCH(0) ;
164 UPZZ(1) ; UPCH(1) ;
165 UPZZ(2) ; UPCH(2) ;
166 UPZZ(3) ; UPCH(3) ;
167 UPZZ(4) ; UPCH(4) ;
168 UPZZ(5) ; UPCH(5) ;
169 UPZZ(6) ; UPCH(6) ; UPLL(7) ;
170 }
171 break ;
172
173 case 8:
174 for( vox=0 ; vox < nv ; vox++ ){
175 zz = (float) vdata[vox] ;
176 UPZZ(0) ; UPCH(0) ;
177 UPZZ(1) ; UPCH(1) ;
178 UPZZ(2) ; UPCH(2) ;
179 UPZZ(3) ; UPCH(3) ;
180 UPZZ(4) ; UPCH(4) ;
181 UPZZ(5) ; UPCH(5) ;
182 UPZZ(6) ; UPCH(6) ;
183 UPZZ(7) ; UPCH(7) ; UPLL(8) ;
184 }
185 break ;
186
187 case 9:
188 for( vox=0 ; vox < nv ; vox++ ){
189 zz = (float) vdata[vox] ;
190 UPZZ(0) ; UPCH(0) ;
191 UPZZ(1) ; UPCH(1) ;
192 UPZZ(2) ; UPCH(2) ;
193 UPZZ(3) ; UPCH(3) ;
194 UPZZ(4) ; UPCH(4) ;
195 UPZZ(5) ; UPCH(5) ;
196 UPZZ(6) ; UPCH(6) ;
197 UPZZ(7) ; UPCH(7) ;
198 UPZZ(8) ; UPCH(8) ; UPLL(9) ;
199 }
200 break ;
201
202 }
203 #endif /* EXPAND_UPDATE */
204
205 (vc->nupdate)++ ; /* June 1995: another update completed! */
206 return ;
207 }
208