1 #include "mrilib.h"
2
3 /*------------------------------------------------------------------------*/
4 /*! Compute first principal component of a bunch of vectors
5 * each vector has its own mean removed
6 * the overall mean of the vectors is not removed
7 * only time series points from ibot..itop (inclusive) are processed,
8 and the output timeseries will be itop-ibot+1 points long
9 * NULL is returned if something stupid happens
10 --------------------------------------------------------------------------*/
11
mri_pcvector(MRI_IMARR * imar,int ibot,int itop)12 MRI_IMAGE * mri_pcvector( MRI_IMARR *imar , int ibot, int itop )
13 {
14 int nx , nvec , ii,jj , npos,nneg ;
15 float *amat , *umat , sum ;
16 float *far ; MRI_IMAGE *tim ;
17
18 if( imar == NULL ) return NULL ;
19 nvec = IMARR_COUNT(imar) ; if( nvec < 1 ) return NULL ;
20 nx = IMARR_SUBIM(imar,0)->nx ; if( nx < 1 ) return NULL ;
21 if( ibot < 0 ) ibot = 0 ;
22 if( itop <= ibot || itop >= nx ) itop = nx-1 ;
23 nx = itop-ibot+1 ; if( nx < 2 ) return NULL ;
24
25 #define A(i,j) amat[(i)+(j)*nx] /* nx X nvec matrix */
26
27 amat = (float *)malloc( sizeof(float)*nx*nvec ) ;
28 umat = (float *)malloc( sizeof(float)*nx ) ;
29
30 for( jj=0 ; jj < nvec ; jj++ ){
31 tim = IMARR_SUBIM(imar,jj) ;
32 far = MRI_FLOAT_PTR(tim) ;
33 for( sum=ii=0 ; ii < nx ; ii++ ){
34 A(ii,jj) = far[ii+ibot] ; sum += A(ii,jj) ;
35 }
36 sum /= nx ; /* remove mean from each vector */
37 for( ii=0 ; ii < nx ; ii++ ) A(ii,jj) -= sum ;
38 }
39
40 jj = first_principal_vectors( nx , nvec , amat , 1 , NULL , umat ) ;
41
42 if( jj <= 0 ){ free(umat); free(amat); return NULL; } /* bad */
43
44 tim = mri_new( nx , 1 , MRI_float ) ; /* all zero */
45 far = MRI_FLOAT_PTR(tim) ;
46 for( ii=0 ; ii < nx ; ii++ ) far[ii] = umat[ii] ;
47
48 /* compute dot products with original vectors,
49 and flip sign if more of them are negative than positive */
50
51 for( npos=nneg=jj=0 ; jj < nvec ; jj++ ){
52 for( sum=ii=0 ; ii < nx ; ii++ ) sum += A(ii,jj)*far[ii] ;
53 if( sum > 0.0f ) npos++ ; else if( sum < 0.0f ) nneg++ ;
54 }
55 if( nneg > npos ){
56 for( ii=0 ; ii < nx ; ii++ ) far[ii] = -far[ii] ;
57 }
58
59 free(umat); free(amat); return tim;
60 }
61
62 /*------------------------------------------------------------------------*/
63
mri_meanvector(MRI_IMARR * imar,int ibot,int itop)64 MRI_IMAGE * mri_meanvector( MRI_IMARR *imar , int ibot, int itop )
65 {
66 float *qar , *far ;
67 int nx,nv,jj,kk , nxout ;
68 MRI_IMAGE *im ;
69
70 if( imar == NULL ) return NULL ;
71
72 nx = IMARR_SUBIM(imar,0)->nx ; nv = IMARR_COUNT(imar) ;
73 if( ibot < 0 ) ibot = 0 ;
74 if( itop <= ibot || itop >= nx ) itop = nx-1 ;
75 nxout = itop-ibot+1 ;
76 im = mri_new( nxout , 1 , MRI_float ) ; far = MRI_FLOAT_PTR(im) ;
77 for( jj=0 ; jj < nv ; jj++ ){
78 qar = MRI_FLOAT_PTR(IMARR_SUBIM(imar,jj)) ;
79 for( kk=0 ; kk < nxout ; kk++ ) far[kk] += qar[kk+ibot] ;
80 }
81 for( kk=0 ; kk < nxout ; kk++ ) far[kk] /= nv ;
82 return im ;
83 }
84
85 /*------------------------------------------------------------------------*/
86 /* code = 0 ==> median
87 1 ==> MAD * 1.4826
88 2 ==> biweight midvariance */
89
mri_MMBvector(MRI_IMARR * imar,int ibot,int itop,int code)90 MRI_IMAGE * mri_MMBvector( MRI_IMARR *imar , int ibot, int itop , int code )
91 {
92 float *qar , *far , *var , med,mad,bmv ;
93 int nx,nv,jj,kk , nxout ;
94 MRI_IMAGE *im ;
95
96 if( imar == NULL ) return NULL ;
97
98 nx = IMARR_SUBIM(imar,0)->nx ; nv = IMARR_COUNT(imar) ;
99 if( nv < 2 ) return NULL ;
100
101 if( ibot < 0 ) ibot = 0 ;
102 if( itop <= ibot || itop >= nx ) itop = nx-1 ;
103 nxout = itop-ibot+1 ;
104 im = mri_new( nxout , 1 , MRI_float ) ; far = MRI_FLOAT_PTR(im) ;
105 var = (float *)malloc(sizeof(float)*nv) ;
106 for( kk=0 ; kk < nxout ; kk++ ){
107 for( jj=0 ; jj < nv ; jj++ ){
108 qar = MRI_FLOAT_PTR(IMARR_SUBIM(imar,jj)) ; var[jj] = qar[kk+ibot] ;
109 }
110 qmedmadbmv_float( nv , var , &med , &mad , &bmv ) ;
111 if( code <= 0 ) far[kk] = med ;
112 else if( code == 1 ) far[kk] = mad * 1.4826f ;
113 else if( code >= 2 ) far[kk] = bmv ;
114 }
115 free(var) ; return im ;
116 }
117