1 #include "mrilib.h"
2
3 #include "despike_inc.c"
4
5 /*-----------------------------------------------------------------*/
6 /* This function computes the "correlation" of vector x
7 with the pair of vectors a and b (presumed to be L2 norm 1).
8 *//*---------------------------------------------------------------*/
9
r2D(UAint n,float * a,float * b,float * x)10 static float r2D( UAint n , float *a , float *b , float *x )
11 {
12 float sax=0.0f , sbx=0.0f , sxx=0.0f ; UAint ii ;
13
14 for( ii=0 ; ii < n ; ii++ ){
15 sax += x[ii]*a[ii] ; /* a dot x */
16 sbx += x[ii]*b[ii] ; /* b dot x */
17 sxx += x[ii]*x[ii] ; /* x dot x */
18 }
19 if( sxx <= 0.0001f ) return 0.0f ;
20 sax = (sax*sax+sbx*sbx)/sxx ; /* multiple correlation squared */
21 #if 0
22 if( sax > 1.0f ) sax = 1.0f ; /* should be oompossibull */
23 #endif
24 return sax ;
25 }
26
27 /*----------------------------------------------------------------*/
28 /* Return the two principal vectors computed during the pvmapping */
29 /*----------------------------------------------------------------*/
30
31 static UAint nvec=0 ;
32 static float *uvec=NULL , *vvec=NULL ;
33 static float ulam=0.0f , vlam=0.0f ;
34
mri_pvmap_get_vecpair(void)35 MRI_IMAGE * mri_pvmap_get_vecpair(void)
36 {
37 MRI_IMAGE *uvim ;
38
39 if( nvec == 0 || uvec == NULL || vvec == NULL ) return NULL ;
40
41 uvim = mri_new( nvec , 2 , MRI_float ) ;
42 memcpy( MRI_FLOAT_PTR(uvim) , uvec , sizeof(float)*nvec ) ;
43 memcpy( MRI_FLOAT_PTR(uvim)+nvec , vvec , sizeof(float)*nvec ) ;
44 return uvim ;
45 }
46
47 /*------------------------------------*/
48 /* Return the pair of singular values */
49 /*------------------------------------*/
50
mri_pvmap_get_lampair(void)51 float_pair mri_pvmap_get_lampair(void)
52 {
53 float_pair uvlam ;
54 uvlam.a = ulam ; uvlam.b = vlam ; return uvlam ;
55 }
56
57 /*----------------------------------------------------------------*/
58 /* Inputs:
59 nx = number of points per time series (>= 9)
60 ny = number of time series (>= 9)
61 vlist[j] = j-th time series (j=0..ny-1)
62 Output:
63 newly malloc-ed array, length ny, of the pvmap intensity
64 (r2D) for each time series
65 *//*-------------------------------------------------------------*/
66
mri_veclist_to_pvmap(UAint nx,UAint ny,float ** vlist)67 static float * mri_veclist_to_pvmap( UAint nx , UAint ny , float **vlist )
68 {
69 float_pair svals ;
70 float *outar ;
71 unsigned short xran[3] ; /* seed for initializing iterations */
72 UAint ii ;
73 static int ncall=-109 ;
74
75 ENTRY("mri_veclist_to_pvmap") ;
76
77 if( vlist == NULL ) RETURN(NULL) ;
78 if( nx < 9 ) RETURN(NULL) ;
79 if( ny < 9 ) RETURN(NULL) ;
80
81 /* get space for the principal vector pair */
82
83 if( nx != nvec || uvec == NULL || vvec == NULL ){
84 uvec = (float *)realloc(uvec,sizeof(float)*nx) ;
85 vvec = (float *)realloc(vvec,sizeof(float)*nx) ;
86 nvec = nx ;
87 }
88
89 /* random number seed used to initialize
90 iterating vectors for the principal vector pair */
91
92 xran[0] = (unsigned short)(nx+ny+73) ;
93 xran[1] = (unsigned short)(nx-ny+473+ncall) ; ncall++ ;
94 xran[2] = (unsigned short)(nx*ny-777) ;
95
96 /* compute the first 2 principal vectors */
97
98 svals = principal_vector_pair( nx , ny , 1 , vlist , /* cf. cs_pv.c */
99 uvec , vvec , NULL , NULL , xran ) ;
100
101 ulam = svals.a ; vlam = svals.b ; /* singular values */
102
103 ININFO_message(" first 2 singular values = %g %g",ulam,vlam) ;
104
105 if( ulam < 0.0f || vlam < 0.0f ) RETURN(NULL) ; /* should not happen */
106
107 outar = (float *)malloc(sizeof(float)*ny) ;
108
109 THD_normalize(nx,uvec) ; /* sum of squares = 1 */
110 THD_normalize(nx,vvec) ;
111
112 /* compute pvmap value for each input vector */
113
114 for( ii=0 ; ii < ny ; ii++ ){
115 outar[ii] = r2D( nx , uvec , vvec , vlist[ii] ) ;
116 }
117
118 RETURN(outar) ;
119 }
120
121 /*----------------------------------------------------------------*/
122 /* Modified extensively to compute in time series length chunks
123 of maximum length NXMAX, then average (if more than one chunk).
124 Longer time series tend to take a very long time to compute.
125 [Rewritten Dec 2021, RWCox]
126 *//*--------------------------------------------------------------*/
127
128 #define NXMAX 1024
129
mri_vec_to_pvmap(MRI_IMAGE * inim)130 MRI_IMAGE * mri_vec_to_pvmap( MRI_IMAGE *inim )
131 {
132 Aint nx , ny , ii , numpart , lenpart , pp , xbot,xtop ;
133 float_pair svals ;
134 MRI_IMAGE *outim ;
135 float **vlist , *pvout=NULL , *pvpart=NULL , *iar ;
136
137 ENTRY("mri_vec_to_pvmap") ;
138
139 if( inim == NULL || inim->kind != MRI_float ) RETURN(NULL) ;
140
141 nx = (Aint)inim->nx ; if( nx < 9 ) RETURN(NULL) ;
142 ny = (Aint)inim->ny ; if( ny < 9 ) RETURN(NULL) ;
143
144 vlist = (float **)malloc(sizeof(float *)*ny) ;
145
146 numpart = (nx-1)/NXMAX + 1 ; /* number of parts */
147 lenpart = nx / numpart ; /* length of each part (except last) */
148
149 iar = MRI_FLOAT_PTR(inim) ;
150
151 for( pp=0 ; pp < numpart ; pp++ ){ /* loop over parts */
152 xbot = lenpart * pp ;
153 xtop = (pp < numpart-1) ? (xbot+lenpart) : nx ; /* last part goes to end */
154 /* pointers to start of each sub-vector for this part */
155 for( ii=0 ; ii < ny ; ii ++ ) vlist[ii] = iar + (xbot+ii*nx) ;
156 /* make PVmap for this part */
157 pvpart = mri_veclist_to_pvmap( (xtop-xbot) , ny , vlist ) ;
158 if( pp == 0 ){ /* first part = save as output */
159 pvout = pvpart ;
160 } else { /* later part = sum into output */
161 for( ii=0 ; ii < ny ; ii++) pvout[ii] += pvpart[ii] ;
162 free(pvpart) ;
163 }
164 }
165
166 free(vlist) ;
167
168 if( numpart > 1 ){ /* output = average of input part PVmaps */
169 float fac = 1.0f / numpart ;
170 for( ii=0 ; ii < ny ; ii++ ) pvout[ii] *= fac ;
171 }
172
173 /* manufacture output image with no data 'empty',
174 then put the output array into it as its data (sneaky, huh?) */
175
176 outim = mri_new_empty( ny , 1 , MRI_float ) ;
177 mri_set_data_pointer( outim , pvout ) ;
178
179 RETURN(outim) ;
180 }
181
182 /*-----------------------------------------------------------------*/
183
THD_dataset_to_pvmap(THD_3dim_dataset * dset,byte * mask)184 MRI_IMAGE * THD_dataset_to_pvmap( THD_3dim_dataset *dset , byte *mask )
185 {
186 int nvox, npt, nmask, ii,jj , polort ;
187 MRI_IMAGE *inim, *tim, *outim ;
188 float *inar, *tar, *outar, *dar ;
189
190 ENTRY("THD_dataset_to_pvmap") ;
191
192 if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
193
194 nvox = DSET_NVOX(dset) ;
195 npt = DSET_NVALS(dset) ;
196 if( nvox < 9 || npt < 9 ) RETURN(NULL) ;
197
198 if( mask != NULL ){
199 nmask = THD_countmask( nvox , mask ) ;
200 if( nmask < 9 ) RETURN(NULL) ;
201 } else {
202 nmask = nvox ;
203 }
204
205 inim = mri_new( npt , nmask , MRI_float ) ;
206 inar = MRI_FLOAT_PTR(inim) ;
207 dar = (float *)malloc(sizeof(float)*npt) ;
208 tar = (float *)malloc(sizeof(float)*npt*3) ;
209
210 DSET_load(dset) ;
211
212 polort = npt / 50 ; /* 24 Apr 2019 */
213 if( polort < 2 ) polort = 2 ; /* change detrending */
214 else if( polort > 20 ) polort = 20 ;
215
216 for( jj=ii=0 ; ii < nvox ; ii++ ){
217 if( mask == NULL || mask[ii] != 0 ){
218 THD_extract_array( ii , dset , 0 , dar ) ;
219 DES_despike25( npt , dar , tar ) ; /* despiking */
220 #if 0
221 THD_cubic_detrend( npt , dar ) ; /* detrending */
222 #else
223 THD_generic_detrend_LSQ( npt , dar , polort , 0,NULL,NULL ) ;
224 #endif
225 THD_normalize( npt , dar ) ; /* L2 normalize */
226 memcpy( inar+jj*npt , dar , sizeof(float)*npt ) ;
227 jj++ ;
228 }
229 }
230
231 free(tar) ; free(dar) ;
232
233 tim = mri_vec_to_pvmap( inim ) ;
234
235 mri_free(inim) ;
236
237 if( nmask == nvox ) RETURN(tim) ;
238
239 outim = mri_new( nvox , 1 , MRI_float ) ;
240 outar = MRI_FLOAT_PTR(outim) ;
241 tar = MRI_FLOAT_PTR(tim) ;
242
243 for( jj=ii=0 ; ii < nvox ; ii++ ){
244 if( mask == NULL || mask[ii] != 0 ) outar[ii] = tar[jj++] ;
245 else outar[ii] = 0.0f ;
246 }
247
248 mri_free(tim) ;
249
250 RETURN(outim) ;
251 }
252