1 #include "mrilib.h"
2
3 /*---------------------------------------------------------------------------*/
4
THD_svdblur(THD_3dim_dataset * inset,byte * mask,float rad,int pdim,int nort,float ** ort)5 THD_3dim_dataset * THD_svdblur( THD_3dim_dataset *inset, byte *mask,
6 float rad, int pdim, int nort, float **ort )
7 {
8 THD_3dim_dataset *outset=NULL ;
9 MCW_cluster *nbhd=NULL ;
10 char *prefix="./SVDblur" ;
11 int kk,nx,ny,nz,nxy,nxyz,nt , nmask ;
12 float na,nb,nc , dx,dy,dz ;
13
14 ENTRY("THD_svdblur") ;
15
16 /*---- error checking and basic setup ----*/
17
18 if( !ISVALID_DSET(inset) ) RETURN(NULL) ;
19 if( mask == NULL ) RETURN(NULL) ;
20 if( rad == 0.0f ) RETURN(NULL) ;
21 if( nort > 0 && ort == NULL ) RETURN(NULL) ;
22
23 nt = DSET_NVALS(inset) ; if( nt < 9+nort ) RETURN(NULL) ;
24
25 nx = DSET_NX(inset) ;
26 ny = DSET_NY(inset) ; nxy = nx*ny ;
27 nz = DSET_NZ(inset) ; nxyz = nxy*nz ;
28
29 nmask = THD_countmask( nxyz , mask ) ; if( nmask < 9 ) RETURN(NULL) ;
30
31 DSET_load(inset) ; if( !DSET_LOADED(inset) ) RETURN(NULL) ;
32
33 /*---- create spherical neighborhood (as a cluster) -----*/
34
35 if( rad < 0.0f ){ dx = dy = dz = 1.0f; rad = -rad; rad = MAX(rad,1.01f); }
36 else { dx = fabsf(DSET_DX(inset));
37 dy = fabsf(DSET_DY(inset));
38 dz = fabsf(DSET_DZ(inset)); }
39 nbhd = MCW_spheremask( dx,dy,dz , rad ) ;
40 MCW_radsort_cluster( nbhd , dx,dy,dz ) ; /* ensures first value is centroid */
41
42 /*---- create output dataset ----*/
43
44 outset = EDIT_empty_copy(inset) ;
45 EDIT_dset_items( outset, ADN_prefix,prefix, ADN_brick_fac,NULL, ADN_none );
46 for( kk=0 ; kk < nt ; kk++ ) /* create zero-filled bricks */
47 EDIT_substitute_brick( outset , kk , MRI_float , NULL ) ;
48
49 /**** the real work now begins ****/
50
51 { int kk , xx,yy,zz , vv,ii ;
52 MRI_IMARR *imar ; MRI_IMAGE *pim ;
53 float *tsar ;
54
55 for( kk=0 ; kk < nxyz ; kk++ ){ /* voxel loop */
56 if( !mask[kk] ) continue ;
57 IJK_TO_THREE( kk , xx,yy,zz , nx,nxy ) ;
58 imar = THD_get_dset_nbhd_array( inset , mask , xx,yy,zz , nbhd ) ;
59 if( imar == NULL ) continue ;
60 if( nort > 0 ){
61 for( vv=0 ; vv < IMARR_COUNT(imar) ; vv++ ){
62 tsar = MRI_FLOAT_PTR(IMARR_SUBIM(imar,vv)) ;
63 THD_generic_detrend_LSQ( nt , tsar , -1 , nort , ort , NULL ) ;
64 }
65 }
66 pim = mri_svdproj( imar , pdim ) ;
67 DESTROY_IMARR(imar) ;
68 if( pim != NULL ){
69 THD_insert_series( kk, outset, nt, MRI_float, MRI_FLOAT_PTR(pim), 0 ) ;
70 mri_free(pim) ;
71 }
72 } /* end of voxel loop */
73 }
74
75 /*** cleanup and exit ***/
76
77 RETURN(outset) ;
78 }
79
80 /*------------------------------------------------------------------------*/
81
THD_get_dset_nbhd_array(THD_3dim_dataset * dset,byte * mask,int xx,int yy,int zz,MCW_cluster * nbhd)82 MRI_IMARR * THD_get_dset_nbhd_array( THD_3dim_dataset *dset, byte *mask,
83 int xx, int yy, int zz, MCW_cluster *nbhd )
84 {
85 MRI_IMARR *imar ;
86 int nvox, *ivox , nx,ny,nz , nxy,nxyz , npt, aa,bb,cc,kk,ii ;
87
88 nx = DSET_NX(dset) ;
89 ny = DSET_NY(dset) ; nxy = nx*ny ;
90 nz = DSET_NZ(dset) ; nxyz = nxy*nz ; npt = nbhd->num_pt ;
91
92 kk = xx + yy*nx + zz*nxy ; if( kk < 0 || kk >= nxyz ) return(NULL) ;
93
94 ivox = (int *)malloc(sizeof(int)*npt) ; nvox = 0 ;
95 for( ii=0 ; ii < npt ; ii++ ){
96 aa = xx + nbhd->i[ii] ; if( aa < 0 || aa >= nx ) continue ;
97 bb = yy + nbhd->j[ii] ; if( bb < 0 || bb >= ny ) continue ;
98 cc = zz + nbhd->k[ii] ; if( cc < 0 || cc >= nz ) continue ;
99 kk = aa + bb*nx + cc*nxy ;
100 if( mask == NULL || mask[kk] ) ivox[nvox++] = kk ;
101 }
102 if( nvox == 0 ){ free(ivox) ; return(NULL) ; } /* no voxels to extract */
103
104 imar = THD_extract_many_series( nvox, ivox, dset ) ;
105 free(ivox) ; return(imar) ;
106 }
107
108 /*------------------------------------------------------------------------*/
109 /* Project the vectors in imar along the principal nev-dimensional
110 subspace; if nev <= 0, just return the first eigenvector instead.
111 --------------------------------------------------------------------------*/
112
mri_svdproj(MRI_IMARR * imar,int nev)113 MRI_IMAGE * mri_svdproj( MRI_IMARR *imar , int nev )
114 {
115 int nx , nvec , ii,jj , itop , doproj=(nev > 0) ;
116 float *far , *xar ; MRI_IMAGE *tim ;
117 float *vnorm=NULL ;
118 register float sum ;
119 float *uvec ;
120
121 if( imar == NULL ) return(NULL) ;
122 nvec = IMARR_COUNT(imar) ; if( nvec < 1 ) return(NULL) ;
123 nx = IMARR_SUBIM(imar,0)->nx ; if( nx < 1 ) return(NULL) ;
124
125 if( nvec == 1 ){ /* trivial case */
126 tim = mri_to_float( IMARR_SUBIM(imar,0) ) ;
127 if( nev <= 0 ) THD_normalize( tim->nx , MRI_FLOAT_PTR(tim) ) ;
128 return(tim) ;
129 }
130
131 #undef U
132 #define U(i,j) uvec[(i)+(j)*nx] /* ditto */
133
134 if( !doproj ) nev = 1 ;
135 uvec = (float *)malloc( sizeof(float)*nx*nev ) ;
136
137 nev = mri_principal_vectors( imar , nev , NULL , uvec ) ;
138 if( nev <= 0 ){ free(uvec); return(NULL); }
139
140 /** create output **/
141
142 tim = mri_new( nx , 1 , MRI_float ) ;
143 far = MRI_FLOAT_PTR(tim) ; /* zero filled */
144 xar = MRI_FLOAT_PTR(IMARR_SUBIM(imar,0)) ; /* central vector */
145
146 if( doproj ){ /* project input time series (1st vector in imar) onto subspace */
147
148 for( jj=0 ; jj < nev ; jj++ ){
149 sum = 0.0f ;
150 for( ii=0 ; ii < nx ; ii++ ) sum += U(ii,jj) * xar[ii] ;
151 for( ii=0 ; ii < nx ; ii++ ) far[ii] += sum * U(ii,jj) ;
152 }
153
154 } else { /* just save first eigenvector */
155
156 sum = 0.0f ;
157 for( ii=0 ; ii < nx ; ii++ ) sum += U(ii,0) * xar[ii] ;
158 if( sum < 0.0f )
159 for( ii=0 ; ii < nx ; ii++ ) far[ii] = -U(ii,0) ;
160 else
161 for( ii=0 ; ii < nx ; ii++ ) far[ii] = U(ii,0) ;
162
163 }
164
165 /** finito **/
166
167 free(uvec); return(tim);
168 }
169 #undef U
170
171 /*------------------------------------------------------------------------*/
172
mri_first_principal_vector(MRI_IMARR * imar)173 MRI_IMAGE * mri_first_principal_vector( MRI_IMARR *imar )
174 {
175 return mri_svdproj( imar , 0 ) ;
176 }
177
178 /*----------------------------------------------------------------------------*/
179 /*! Compute the nvec principal singular vectors of a set of m columns, each
180 of length n, stored in the image array imar.
181
182 The singular values (largest to smallest) are stored in sval, and
183 the left singular vectors [first nvec columns of U in X = U S V'] are
184 stored into uvec[i+j*n] for i=0..n-1, j=0..nvec-1. These columns
185 are L2-normalized and orthogonal.
186
187 The return value is the number of vectors computed. If the return
188 value is not positive, something bad happened. Normally, the return
189 value would be the same as nvec, but it cannot be larger than MIN(n,m).
190
191 If sval==NULL, then the output into sval is skipped.
192 If uval==NULL, then the output into uval is skipped.
193 If both are NULL, exactly why did you want to call this function?
194 *//*--------------------------------------------------------------------------*/
195
mri_principal_vectors(MRI_IMARR * imar,int nvec,float * sval,float * uvec)196 int mri_principal_vectors( MRI_IMARR *imar, int nvec, float *sval, float *uvec )
197 {
198 int nn , mm , nsym , ii,jj,kk,qq ;
199 double *asym , *deval ;
200 float **xpt ;
201 register double sum , qsum ; register float *xj , *xk ;
202
203 if( imar == NULL || IMARR_COUNT(imar) < 1 ) return(-555) ;
204
205 nn = IMARR_SUBIM(imar,0)->nx ;
206 mm = IMARR_COUNT(imar) ;
207
208 nsym = MIN(nn,mm) ; /* size of the symmetric matrix to create */
209
210 if( nsym < 1 || (uvec == NULL && sval == NULL) ) return(-666) ;
211
212 /* pointers to each vector */
213
214 xpt = (float **)malloc(sizeof(float *)*mm) ;
215 for( kk=0 ; kk < mm ; kk++ ) xpt[kk] = MRI_FLOAT_PTR(IMARR_SUBIM(imar,kk)) ;
216
217 /* number of eigenpairs to compute */
218
219 if( nvec > nsym ) nvec = nsym ; /* can't compute more vectors than nsym */
220 else if( nvec <= 0 ) nvec = 1 ;
221
222 #pragma omp critical (MALLOC)
223 { asym = (double *)malloc(sizeof(double)*nsym*nsym) ; /* symmetric matrix */
224 deval = (double *)malloc(sizeof(double)*nsym) ; /* its eigenvalues */
225 }
226
227 #undef A
228 #define A(i,j) asym[(i)+(j)*nsym]
229
230 /** setup matrix to eigensolve: choose smaller of [X]'[X] and [X][X]' **/
231 /** since [X] is n x m, [X]'[X] is m x m and [X][X]' is n x n **/
232
233 if( nn > mm ){ /* more rows than columns: */
234 /* so [A] = [X]'[X] = m x m */
235 int n1 = nn-1 ;
236 for( jj=0 ; jj < mm ; jj++ ){
237 xj = xpt[jj] ;
238 for( kk=0 ; kk <= jj ; kk++ ){
239 sum = 0.0 ; xk = xpt[kk] ;
240 for( ii=0 ; ii < n1 ; ii+=2 ) sum += xj[ii]*xk[ii] + xj[ii+1]*xk[ii+1];
241 if( ii == n1 ) sum += xj[ii]*xk[ii] ;
242 A(jj,kk) = sum ; if( kk < jj ) A(kk,jj) = sum ;
243 }
244 }
245
246 } else { /* more columns than rows: */
247 /* so [A] = [X][X]' = n x n */
248 float *xt ; int m1=mm-1 ;
249 #pragma omp critical (MALLOC)
250 xt = (float *)malloc(sizeof(float)*nn*mm) ;
251 for( jj=0 ; jj < mm ; jj++ ){ /* form [X]' into array xt */
252 for( ii=0 ; ii < nn ; ii++ ) xt[jj+ii*mm] = xpt[jj][ii] ;
253 }
254
255 for( jj=0 ; jj < nn ; jj++ ){
256 xj = xt + jj*mm ;
257 for( kk=0 ; kk <= jj ; kk++ ){
258 sum = 0.0 ; xk = xt + kk*mm ;
259 for( ii=0 ; ii < m1 ; ii+=2 ) sum += xj[ii]*xk[ii] + xj[ii+1]*xk[ii+1];
260 if( ii == m1 ) sum += xj[ii]*xk[ii] ;
261 A(jj,kk) = sum ; if( kk < jj ) A(kk,jj) = sum ;
262 }
263 }
264
265 #pragma omp critical (MALLOC)
266 free(xt) ; /* don't need this no more */
267 }
268
269 /** compute the nvec eigenvectors corresponding to largest eigenvalues **/
270 /** these eigenvectors are stored on top of first nvec columns of asym **/
271
272 ii = symeig_irange( nsym, asym, deval, nsym-nvec, nsym-1, (uvec==NULL) ) ;
273
274 if( ii != 0 ){
275 #pragma omp critical (MALLOC)
276 { free(deval) ; free(asym) ; free(xpt) ; }
277 return(-33333) ; /* eigensolver failed!? */
278 }
279
280 /** Store singular values (sqrt of eigenvalues), if desired:
281 Note that symeig_irange returns things smallest to largest,
282 but we want largest to smallest, so have to reverse the order **/
283
284 if( sval != NULL ){
285 for( jj=0 ; jj < nvec ; jj++ ){
286 sum = deval[nvec-1-jj] ;
287 sval[jj] = (sum <= 0.0) ? 0.0 : sqrt(sum) ;
288 }
289 }
290
291 /** if no output vectors desired, we are done done done!!! **/
292
293 if( uvec == NULL ){
294 #pragma omp critical (MALLOC)
295 { free(deval) ; free(asym) ; free(xpt) ; }
296 return(nvec) ;
297 }
298
299 /** SVD is [X] = [U] [S] [V]', where [U] = desired output vectors
300
301 case n <= m: [A] = [X][X]' = [U] [S][S]' [U]'
302 so [A][U] = [U] [S][S]'
303 so eigenvectors of [A] are just [U]
304
305 case n > m: [A] = [X]'[X] = [V] [S]'[S] [V]'
306 so [A][V] = [V] [S'][S]
307 so eigenvectors of [A] are [V], but we want [U]
308 note that [X][V] = [U] [S]
309 so pre-multiplying each column vector in [V] by matrix [X]
310 will give the corresponding column in [U], but scaled;
311 below, just L2-normalize the column to get output vector **/
312
313 if( nn <= mm ){ /* copy eigenvectors into output directly */
314 /* (e.g., more vectors than time points) */
315 for( jj=0 ; jj < nvec ; jj++ ){
316 qq = nvec-1-jj ; /* eigenvalues are in reversed order */
317 for( ii=0 ; ii < nn ; ii++ )
318 uvec[ii+jj*nn] = (float)asym[ii+qq*nn] ;
319 }
320
321 } else { /* n > m: transform eigenvectors to get left singular vectors */
322 /* (e.g., more time points than vectors) */
323
324 for( jj=0 ; jj < nvec ; jj++ ){
325 qq = nvec-1-jj ; qsum = 0.0 ; /* eigenvalues are in reversed order */
326 for( ii=0 ; ii < nn ; ii++ ){
327 sum = 0.0 ;
328 for( kk=0 ; kk < mm ; kk++ ) sum += xpt[kk][ii] * asym[kk+qq*mm] ;
329 uvec[ii+jj*nn] = sum ; qsum += sum*sum ;
330 }
331 if( qsum > 0.0 ){ /* L2 normalize */
332 register float fac ;
333 fac = (float)(1.0/sqrt(qsum)) ;
334 for( ii=0 ; ii < nn ; ii++ ) uvec[ii+jj*nn] *= fac ;
335 }
336 }
337 }
338
339 /** free at last!!! **/
340
341 #pragma omp critical (MALLOC)
342 { free(deval) ; free(asym) ; free(xpt) ; }
343
344 return(nvec) ;
345 }
346 #undef A
347
348 /*----------------------------------------------------------------------------*/
349
mri_average_vector(MRI_IMARR * imar)350 MRI_IMAGE * mri_average_vector( MRI_IMARR *imar )
351 {
352 int nx , nvec , jj ;
353 float *far , *xar ; MRI_IMAGE *tim ;
354 register int ii ;
355
356 if( imar == NULL ) return NULL ;
357 nvec = IMARR_COUNT(imar) ; if( nvec < 1 ) return NULL ;
358 nx = IMARR_SUBIM(imar,0)->nx ; if( nx < 1 ) return NULL ;
359
360 /** create output **/
361
362 tim = mri_new( nx , 1 , MRI_float ) ;
363 far = MRI_FLOAT_PTR(tim) ; /* zero filled */
364
365 /** sum them up, then average **/
366
367 for( jj=0 ; jj < nvec ; jj++ ){
368 xar = MRI_FLOAT_PTR( IMARR_SUBIM(imar,jj) ) ;
369 for( ii=0 ; ii < nx ; ii++ ) far[ii] += xar[ii] ;
370 }
371 if( nvec > 1 ){
372 register float fac ;
373 fac = 1.0f / nvec ;
374 for( ii=0 ; ii < nx ; ii++ ) far[ii] *= fac ;
375 }
376
377 return tim;
378 }
379