1 #include "mrilib.h"
2
3 /****** functions used mostly in 1dmatcalc.c, but elsewhere as well ******/
4 /***** these functions operate on matrices stored in 2D float images *****/
5
6 /*----------------------------------------------------------------------------*/
7
8 static int verb_rpn = 0 ;
mri_matrix_evalrpn_verb(int i)9 void mri_matrix_evalrpn_verb(int i){ verb_rpn = i ; }
10
11 /*-----------------------------------------------------------------------*/
12
mri_matrix_print(FILE * fp,MRI_IMAGE * ima,char * label)13 void mri_matrix_print( FILE *fp , MRI_IMAGE *ima , char *label )
14 {
15 int ii,jj , nr,nc , ipr ; float *amat , val ;
16
17 if( ima == NULL ) return ;
18
19 nr = ima->nx ; nc = ima->ny ; amat = MRI_FLOAT_PTR(ima) ;
20
21 #undef A
22 #define A(i,j) amat[(i)+(j)*nr] /* nr X nc */
23
24 for( ii=0 ; ii < ima->nvox ; ii++ ){
25 val = (int)amat[ii] ;
26 if( val != amat[ii] || fabsf(val) > 99.0f ) break ;
27 }
28 ipr = (ii == ima->nvox ) ;
29
30 if( fp == NULL ) fp = stdout ;
31 if( label != NULL ) fprintf(fp,"Matrix [%dX%d] %s\n",nr,nc,label) ;
32 for( ii=0 ; ii < nr ; ii++ ){
33 for( jj=0 ; jj < nc ; jj++ ){
34 if( ipr ) fprintf(fp," %3d" , (int)A(ii,jj) ) ;
35 else fprintf(fp," %11.5g" , A(ii,jj) ) ;
36 }
37 fprintf(fp,"\n") ;
38 }
39 fprintf(fp,"\n") ; fflush(fp) ; return ;
40 }
41
42 /*-----------------------------------------------------------------------*/
43 /*! Compute the product of two matrices, stored in 2D float images.
44 -------------------------------------------------------------------------*/
45
mri_matrix_mult(MRI_IMAGE * ima,MRI_IMAGE * imb)46 MRI_IMAGE * mri_matrix_mult( MRI_IMAGE *ima , MRI_IMAGE *imb )
47 {
48 int nr , nc , mm ;
49 MRI_IMAGE *imc ;
50 float *amat , *bmat , *cmat , sum ;
51 int ii,jj,kk ;
52
53 ENTRY("mri_matrix_mult") ;
54
55 if( ima == NULL || imb == NULL ) RETURN( NULL );
56 if( ima->kind != MRI_float || imb->kind != MRI_float ) RETURN( NULL );
57
58 if( ima->nx == 1 && ima->ny == 1 ){ /** scalar multiply? **/
59 amat = MRI_FLOAT_PTR(ima) ;
60 imc = mri_matrix_scale( amat[0] , imb ) ;
61 RETURN(imc) ;
62 } else if( imb->nx == 1 && imb->ny == 1 ){
63 bmat = MRI_FLOAT_PTR(imb) ;
64 imc = mri_matrix_scale( bmat[0] , ima ) ;
65 RETURN(imc) ;
66 }
67
68 nr = ima->nx ; mm = ima->ny ; nc = imb->ny ;
69 if( imb->nx != mm ){
70 ERROR_message("mri_matrix_mult( %d X %d , %d X %d )?",
71 ima->nx , ima->ny , imb->nx , imb->ny ) ;
72 RETURN( NULL );
73 }
74
75 #undef A
76 #undef B
77 #undef C
78 #define A(i,j) amat[(i)+(j)*nr] /* nr X mm */
79 #define B(i,j) bmat[(i)+(j)*mm] /* mm X nc */
80 #define C(i,j) cmat[(i)+(j)*nr] /* nr X nc */
81
82 imc = mri_new( nr , nc , MRI_float ) ;
83 amat = MRI_FLOAT_PTR(ima); bmat = MRI_FLOAT_PTR(imb);
84 cmat = MRI_FLOAT_PTR(imc);
85
86 for( jj=0 ; jj < nc ; jj++ ){
87 for( ii=0 ; ii < nr ; ii++ ){
88 sum = 0.0f ;
89 for( kk=0 ; kk < mm ; kk++ ) sum += A(ii,kk)*B(kk,jj) ;
90 C(ii,jj) = sum ;
91 }}
92
93 RETURN( imc );
94 }
95
96 /*-----------------------------------------------------------------------*/
97 /*! Compute the product of two matrices, stored in 2D float images.
98 The first matrix is transposed.
99 -------------------------------------------------------------------------*/
100
mri_matrix_multranA(MRI_IMAGE * ima,MRI_IMAGE * imb)101 MRI_IMAGE * mri_matrix_multranA( MRI_IMAGE *ima , MRI_IMAGE *imb )
102 {
103 int nr , nc , mm ;
104 MRI_IMAGE *imc ;
105 float *amat , *bmat , *cmat , sum ;
106 int ii,jj,kk ;
107
108 ENTRY("mri_matrix_multranA") ;
109
110 if( ima == NULL || imb == NULL ) RETURN( NULL );
111 if( ima->kind != MRI_float || imb->kind != MRI_float ) RETURN( NULL );
112
113 nr = ima->ny ; mm = ima->nx ; nc = imb->ny ;
114 if( imb->nx != mm ){
115 ERROR_message("mri_matrix_multranA( %d X %d , %d X %d )?",
116 ima->nx , ima->ny , imb->nx , imb->ny ) ;
117 RETURN( NULL );
118 }
119
120 #undef A
121 #undef B
122 #undef C
123 #define A(i,j) amat[(i)+(j)*mm] /* mm X nr */
124 #define B(i,j) bmat[(i)+(j)*mm] /* mm X nc */
125 #define C(i,j) cmat[(i)+(j)*nr] /* nr X nc */
126
127 imc = mri_new( nr , nc , MRI_float ) ;
128 amat = MRI_FLOAT_PTR(ima); bmat = MRI_FLOAT_PTR(imb);
129 cmat = MRI_FLOAT_PTR(imc);
130
131 for( jj=0 ; jj < nc ; jj++ ){
132 for( ii=0 ; ii < nr ; ii++ ){
133 sum = 0.0f ;
134 for( kk=0 ; kk < mm ; kk++ ) sum += A(kk,ii)*B(kk,jj) ;
135 C(ii,jj) = sum ;
136 }}
137
138 RETURN( imc );
139 }
140
141 /*-----------------------------------------------------------------------*/
142 /*! Compute the product of two matrices, stored in 2D float images.
143 The second matrix is transposed.
144 -------------------------------------------------------------------------*/
145
mri_matrix_multranB(MRI_IMAGE * ima,MRI_IMAGE * imb)146 MRI_IMAGE * mri_matrix_multranB( MRI_IMAGE *ima , MRI_IMAGE *imb )
147 {
148 int nr , nc , mm ;
149 MRI_IMAGE *imc ;
150 float *amat , *bmat , *cmat , sum ;
151 int ii,jj,kk ;
152
153 ENTRY("mri_matrix_multranB") ;
154
155 if( ima == NULL || imb == NULL ) RETURN( NULL );
156 if( ima->kind != MRI_float || imb->kind != MRI_float ) RETURN( NULL );
157
158 nr = ima->nx ; mm = ima->ny ; nc = imb->nx ;
159 if( imb->ny != mm ){
160 ERROR_message("mri_matrix_multranB( %d X %d , %d X %d )?",
161 ima->nx , ima->ny , imb->nx , imb->ny ) ;
162 RETURN( NULL );
163 }
164
165 #undef A
166 #undef B
167 #undef C
168 #define A(i,j) amat[(i)+(j)*nr] /* nr X mm */
169 #define B(i,j) bmat[(i)+(j)*nc] /* nc X mm */
170 #define C(i,j) cmat[(i)+(j)*nr] /* nr X nc */
171
172 imc = mri_new( nr , nc , MRI_float ) ;
173 amat = MRI_FLOAT_PTR(ima); bmat = MRI_FLOAT_PTR(imb);
174 cmat = MRI_FLOAT_PTR(imc);
175
176 #if 0
177 INFO_message("mri_matrix_multranB( %d X %d , %d X %d )",
178 ima->nx , ima->ny , imb->nx , imb->ny ) ;
179 #endif
180
181 for( jj=0 ; jj < nc ; jj++ ){
182 for( ii=0 ; ii < nr ; ii++ ){
183 sum = 0.0f ;
184 for( kk=0 ; kk < mm ; kk++ ) sum += A(ii,kk)*B(jj,kk) ;
185 C(ii,jj) = sum ;
186 }}
187
188 RETURN( imc );
189 }
190
191 /*-----------------------------------------------------------------------*/
192 /*! Multiply matrix in ima by scalar fa, matrix in imb by scalar fb,
193 and add the results. Used for general purpose matrix add/subtract.
194 -------------------------------------------------------------------------*/
195
mri_matrix_sadd(float fa,MRI_IMAGE * ima,float fb,MRI_IMAGE * imb)196 MRI_IMAGE * mri_matrix_sadd( float fa, MRI_IMAGE *ima, float fb, MRI_IMAGE *imb )
197 {
198 int ii , nrc ;
199 MRI_IMAGE *imc ;
200 float *amat , *bmat , *cmat ;
201
202 ENTRY("mri_matrix_sadd") ;
203
204 if( ima == NULL || imb == NULL ) RETURN( NULL );
205 if( ima->kind != MRI_float || imb->kind != MRI_float ) RETURN( NULL );
206
207 nrc = ima->nvox ;
208 if( imb->nvox != nrc ){
209 ERROR_message("mri_matrix_sadd( %d X %d , %d X %d ) ?",
210 ima->nx,ima->ny , imb->nx,imb->ny ) ;
211 RETURN(NULL) ;
212 }
213
214 imc = mri_new_conforming(ima,MRI_float) ;
215 amat = MRI_FLOAT_PTR(ima); bmat = MRI_FLOAT_PTR(imb);
216 cmat = MRI_FLOAT_PTR(imc);
217 for( ii=0 ; ii < nrc ; ii++ ) cmat[ii] = fa*amat[ii] + fb*bmat[ii] ;
218
219 RETURN(imc) ;
220 }
221
222 /*-----------------------------------------------------------------------*/
223 /*! Multiply matrix in ima by scalar fa. */
224
mri_matrix_scale(float fa,MRI_IMAGE * ima)225 MRI_IMAGE * mri_matrix_scale( float fa, MRI_IMAGE *ima )
226 {
227 int ii , nrc ;
228 MRI_IMAGE *imc ;
229 float *amat , *cmat ;
230
231 ENTRY("mri_matrix_scale") ;
232
233 if( ima == NULL ) RETURN( NULL );
234 if( ima->kind != MRI_float ) RETURN( NULL );
235
236 nrc = ima->nvox ;
237 imc = mri_new_conforming(ima,MRI_float ) ;
238 amat = MRI_FLOAT_PTR(ima);
239 cmat = MRI_FLOAT_PTR(imc);
240 for( ii=0 ; ii < nrc ; ii++ ) cmat[ii] = fa*amat[ii] ;
241
242 RETURN(imc) ;
243 }
244
245 /*----------------------------------------------------------------------------*/
246
mri_matrix_svals(MRI_IMAGE * imc)247 MRI_IMAGE * mri_matrix_svals( MRI_IMAGE *imc ) /* 05 May 2020 */
248 {
249 float *rmat ;
250 int m , n , ii,jj,kk ;
251 double *amat , *sval , smax ;
252 MRI_IMAGE *imp=NULL ; float *pmat ;
253
254 ENTRY("mri_matrix_svals") ;
255
256 if( imc == NULL || imc->kind != MRI_float ) RETURN( NULL );
257 m = imc->nx ; /* number of rows in input */
258 n = imc->ny ; /* number of columns */
259
260 #undef R
261 #undef A
262 #define R(i,j) rmat[(i)+(j)*m] /* i=0..m-1 , j=0..n-1 */
263 #define A(i,j) amat[(i)+(j)*m] /* i=0..m-1 , j=0..n-1 */
264
265 rmat = MRI_FLOAT_PTR(imc) ;
266
267 imp = mri_new( n , 1 , MRI_float ) ;
268 pmat = MRI_FLOAT_PTR(imp) ;
269
270 if( n == 1 && m == 1 ){
271 pmat[0] = fabsf(rmat[0]) ; RETURN(imp) ; /* trivial case */
272 }
273
274 amat = (double *)calloc( sizeof(double),m*n ) ; /* input matrix */
275
276 for( ii=0 ; ii < m ; ii++ )
277 for( jj=0 ; jj < n ; jj++ ) A(ii,jj) = R(ii,jj) ;
278
279 sval = (double *)calloc( sizeof(double),n );
280 svd_double( m , n , amat , sval , NULL,NULL ) ;
281
282 free(amat) ; /* done with this */
283
284 for( ii=0 ; ii < n ; ii++ ) pmat[ii] = (float)sval[ii] ;
285
286 free(sval) ;
287 RETURN(imp) ;
288 }
289
290 /*----------------------------------------------------------------------------*/
291
mri_matrix_svprod(MRI_IMAGE * imc,int do_log)292 MRI_IMAGE * mri_matrix_svprod( MRI_IMAGE *imc , int do_log ) /* 05 May 2020 */
293 {
294 float *rmat ;
295 int m , n , ii,jj,kk ;
296 double *amat , *sval , smax ;
297 MRI_IMAGE *imp=NULL ; float *pmat ;
298
299 ENTRY("mri_matrix_svprod") ;
300
301 if( imc == NULL || imc->kind != MRI_float ) RETURN( NULL );
302 m = imc->nx ; /* number of rows in input */
303 n = imc->ny ; /* number of columns */
304
305 #undef R
306 #undef A
307 #define R(i,j) rmat[(i)+(j)*m] /* i=0..m-1 , j=0..n-1 */
308 #define A(i,j) amat[(i)+(j)*m] /* i=0..m-1 , j=0..n-1 */
309 #define P(i,j) pmat[(i)+(j)*n] /* i=0..n-1 , j=0..m-1 */
310
311 rmat = MRI_FLOAT_PTR(imc) ;
312
313 imp = mri_new( 1 , 1 , MRI_float ) ;
314 pmat = MRI_FLOAT_PTR(imp) ;
315
316 if( n == 1 && m == 1 ){
317 pmat[0] = fabsf(rmat[0]) ; RETURN(imp) ; /* trivial case */
318 }
319
320 amat = (double *)calloc( sizeof(double),m*n ) ; /* input matrix */
321
322 for( ii=0 ; ii < m ; ii++ )
323 for( jj=0 ; jj < n ; jj++ ) A(ii,jj) = R(ii,jj) ;
324
325 sval = (double *)calloc( sizeof(double),n );
326 svd_double( m , n , amat , sval , NULL,NULL ) ;
327
328 free(amat) ; /* done with this */
329
330 /* find largest singular value */
331
332 smax = sval[0] ;
333 for( ii=1 ; ii < n ; ii++ ) if( sval[ii] > smax ) smax = sval[ii] ;
334
335 if( smax <= 0.0 ){ /* this is bad */
336 static int first = 1 ;
337 #pragma omp critical (STDERR)
338 { if( first ) ERROR_message("SVD fails in mri_matrix_svprod()!\n"); }
339 free(sval); first = 0;
340 RETURN(imp);
341 }
342
343 if( do_log ){
344 smax = 0.0 ;
345 for( ii=0 ; ii < n ; ii++ )
346 if( sval[ii] > 0.0 ) smax += log(sval[ii]) ;
347 } else {
348 smax = 1.0 ;
349 for( ii=0 ; ii < n ; ii++ )
350 if( sval[ii] > 0.0 ) smax *= sval[ii] ;
351 }
352
353 free(sval) ;
354 pmat[0] = (float)smax ; RETURN(imp) ;
355 }
356
357 /*----------------------------------------------------------------------------*/
358 static int force_svd = 0 ;
mri_matrix_psinv_svd(int i)359 void mri_matrix_psinv_svd( int i ){ force_svd = i; }
360 /*----------------------------------------------------------------------------*/
361 /*! Compute the pseudo-inverse of a matrix stored in a 2D float image.
362 If the input is mXn, the output is nXm. wt[] is an optional array
363 of positive weights, m of them. The result can be used to solve
364 the weighted least squares problem
365 [imc] [b] = [v]
366 where [b] is an n-vector and [v] is an m-vector, where m >= n.
367 If alpha > 0, then the actual matrix calculated is
368 -1
369 [imc' imc + alpha I] imc' (where ' = transpose)
370
371 This result can be used to solve the penalized least squares problem
372
373 ( )
374 min ( LQ{ [imc] [b] - [v] } + alpha LQ{ [b] } )
375 [b] ( )
376
377 where LQ{ [x] } is the sum of squares of the elements of vector [x].
378 The 'penalty' consists of trying to keep the elements of [b] small.
379
380 Note that matrices are stored in column-major order in the 2D image arrays!
381 If m < n, then the SVD solution is tried immediately, since the Choleski
382 method makes no sense in this case. You can force the SVD method to
383 be used by calling mri_matrix_psinv_svd().
384 *//*--------------------------------------------------------------------------*/
385
mri_matrix_psinv(MRI_IMAGE * imc,float * wt,float alpha)386 MRI_IMAGE * mri_matrix_psinv( MRI_IMAGE *imc , float *wt , float alpha )
387 {
388 float *rmat ;
389 int m , n , ii,jj,kk ;
390 double *amat , *umat , *vmat , *sval , *xfac , smax,del,ww , alp ;
391 MRI_IMAGE *imp=NULL ; float *pmat ;
392 register double sum ;
393 int do_svd = (force_svd || AFNI_yesenv("AFNI_PSINV_SVD")) ;
394 int *kbot=NULL,*ktop=NULL , ibot,itop,jbot,jtop , mn ; /* 12 Feb 2009 */
395
396 ENTRY("mri_matrix_psinv") ;
397
398 if( imc == NULL || imc->kind != MRI_float ) RETURN( NULL );
399 m = imc->nx ; /* number of rows in input */
400 n = imc->ny ; /* number of columns */
401 if( PRINT_TRACING ){ char str[222]; sprintf(str,"m=%d n=%d",m,n); STATUS(str); }
402
403 /* deal with a single vector (of length m) [30 Apr 2009] */
404
405 rmat = MRI_FLOAT_PTR(imc) ;
406
407 if( n == 1 ){
408 for( sum=0.0,ii=0 ; ii < m ; ii++ ) sum += rmat[ii]*rmat[ii] ;
409 imp = mri_new( 1 , m , MRI_float ) ;
410 if( sum > 0.0 ){
411 sum = 1.0 / sum ; pmat = MRI_FLOAT_PTR(imp) ;
412 for( ii=0 ; ii < m ; ii++ ) pmat[ii] = sum * rmat[ii] ;
413 }
414 RETURN(imp) ;
415 }
416
417 /* OK, have a real matrix to handle here */
418
419 amat = (double *)calloc( sizeof(double),m*n ) ; /* input matrix */
420 xfac = (double *)calloc( sizeof(double),n ) ; /* column norms of [a] */
421
422 if( amat == NULL || xfac == NULL ){ /* 29 Dec 2008 */
423 ERROR_message("mri_matrix_psinv: can't malloc matrix workspace!") ;
424 if( amat != NULL ) free(amat) ;
425 if( xfac != NULL ) free(xfac) ;
426 RETURN(NULL) ;
427 }
428
429 #undef PSINV_EPS
430 #define PSINV_EPS 1.e-12
431
432 alp = (alpha <= 0.0f) ? PSINV_EPS : (double)alpha ;
433
434 #undef R
435 #undef A
436 #undef P
437 #undef U
438 #undef V
439 #define R(i,j) rmat[(i)+(j)*m] /* i=0..m-1 , j=0..n-1 */
440 #define A(i,j) amat[(i)+(j)*m] /* i=0..m-1 , j=0..n-1 */
441 #define P(i,j) pmat[(i)+(j)*n] /* i=0..n-1 , j=0..m-1 */
442 #define U(i,j) umat[(i)+(j)*m]
443 #define V(i,j) vmat[(i)+(j)*n]
444
445 /* copy input matrix into amat */
446
447 STATUS("copy matrix") ;
448 for( ii=0 ; ii < m ; ii++ )
449 for( jj=0 ; jj < n ; jj++ ) A(ii,jj) = R(ii,jj) ;
450
451 /* weight rows? */
452
453 if( wt != NULL ){
454 for( ii=0 ; ii < m ; ii++ ){
455 ww = wt[ii] ;
456 if( ww > 0.0 ) for( jj=0 ; jj < n ; jj++ ) A(ii,jj) *= ww ;
457 }
458 }
459
460 /* scale each column to have norm 1 */
461
462 STATUS("scale matrix") ;
463 for( jj=0 ; jj < n ; jj++ ){
464 for( sum=0.0,ii=0 ; ii < m ; ii++ ) sum += A(ii,jj)*A(ii,jj) ;
465 if( sum > 0.0 ) sum = 1.0/sqrt(sum) ; else do_svd = 1 ;
466 xfac[jj] = sum ;
467 if( sum > 0.0 ){
468 for( ii=0 ; ii < m ; ii++ ) A(ii,jj) *= sum ;
469 }
470 }
471
472 /*** computations follow, via SVD or Choleski ***/
473
474 vmat = (double *)calloc( sizeof(double),n*n );
475
476 if( vmat == NULL ){ /* 29 Dec 2008 */
477 ERROR_message("mri_matrix_psinv: can't malloc vmat workspace!") ;
478 free(amat) ; free(xfac) ; RETURN(NULL) ;
479 }
480
481 if( do_svd || m < n ) goto SVD_PLACE ;
482
483 /*** Try the Choleski method first ***/
484
485 mn = MAX(m,n) ;
486 kbot = (int *)malloc(sizeof(int)*mn) ; /* start index of each column */
487 ktop = (int *)malloc(sizeof(int)*mn) ; /* end index */
488 for( ii=0 ; ii < n ; ii++ ){
489 for( kk=0 ; kk < m && A(kk,ii) == 0.0 ; kk++ ) ; /*nada*/
490 kbot[ii] = kk ;
491 for( kk=m-1 ; kk >= 0 && A(kk,ii) == 0.0 ; kk-- ) ; /*nada*/
492 ktop[ii] = kk ;
493 }
494
495 STATUS("form normal eqns") ;
496 for( ii=0 ; ii < n ; ii++ ){ /* form normal equations */
497 if( ii%1000==999 ) STATUS("999") ;
498 ibot = kbot[ii] ; itop = ktop[ii] ;
499 for( jj=0 ; jj <= ii ; jj++ ){
500 jbot = kbot[jj] ; if( ibot > jbot ) jbot = ibot ;
501 jtop = ktop[jj] ; if( itop < jtop ) jtop = itop ;
502 sum = 0.0 ;
503 for( kk=jbot ; kk <= jtop ; kk++ ) sum += A(kk,ii) * A(kk,jj) ;
504 V(ii,jj) = sum ;
505 }
506 V(ii,ii) += alp ; /* note V(ii,ii)==1 before this */
507 }
508
509 free(ktop) ;
510
511 #if 0
512 fprintf(stderr,"mri_psinv: V matrix (%dx%d)\n",n,n) ;
513 for( ii=0 ; ii < n ; ii++ ){
514 fprintf(stderr,"%2d:",ii) ;
515 for( jj=0 ; jj <= ii ; jj++ ) fprintf(stderr," %11.4g",V(ii,jj)) ;
516 fprintf(stderr,"\n") ;
517 }
518 #endif
519
520 /* Choleski factor V in place */
521
522 for( ii=0 ; ii < n ; ii++ ){
523 for( kk=0 ; kk < ii && V(ii,kk) == 0.0 ; kk++ ) ; /*nada*/
524 kbot[ii] = kk ;
525 }
526
527 STATUS("Choleski") ;
528 for( ii=0 ; ii < n ; ii++ ){
529 if( ii%1000==999 ) STATUS("999") ;
530 ibot = kbot[ii] ;
531 for( jj=0 ; jj < ii ; jj++ ){
532 jbot = kbot[jj] ; if( ibot > jbot ) jbot = ibot ;
533 sum = V(ii,jj) ;
534 for( kk=jbot ; kk < jj ; kk++ ) sum -= V(ii,kk) * V(jj,kk) ;
535 V(ii,jj) = sum / V(jj,jj) ;
536 }
537 sum = V(ii,ii) ;
538 for( kk=ibot ; kk < ii ; kk++ ) sum -= V(ii,kk) * V(ii,kk) ;
539 if( sum < PSINV_EPS ){
540 static int first=1 ;
541 #pragma omp critical (STDERR)
542 { if( first ) WARNING_message("Choleski fails in mri_matrix_psinv()!\n");
543 first = 0 ; do_svd = 1 ; }
544 goto SVD_PLACE ;
545 }
546 V(ii,ii) = sqrt(sum) ;
547 }
548
549 free(kbot) ;
550
551 /* create pseudo-inverse from what's now in V */
552
553 imp = mri_new( n , m , MRI_float ) ; /* recall that m > n */
554 pmat = MRI_FLOAT_PTR(imp) ;
555
556 sval = (double *)calloc( sizeof(double),n ) ; /* row #jj of A */
557
558 STATUS("psinv from Choleski") ;
559 for( jj=0 ; jj < m ; jj++ ){
560 if( jj%1000==999 ) STATUS("999") ;
561 for( ii=0 ; ii < n ; ii++ ) sval[ii] = A(jj,ii) ; /* extract row */
562
563 for( ii=0 ; ii < n ; ii++ ){ /* forward solve */
564 sum = sval[ii] ;
565 for( kk=0 ; kk < ii ; kk++ ) sum -= V(ii,kk) * sval[kk] ;
566 sval[ii] = sum / V(ii,ii) ;
567 }
568 for( ii=n-1 ; ii >= 0 ; ii-- ){ /* backward solve */
569 sum = sval[ii] ;
570 for( kk=ii+1 ; kk < n ; kk++ ) sum -= V(kk,ii) * sval[kk] ;
571 sval[ii] = sum / V(ii,ii) ;
572 }
573
574 for( ii=0 ; ii < n ; ii++ ) P(ii,jj) = (float)sval[ii] ;
575 }
576 free(amat); free(vmat); free(sval);
577 goto RESCALE_PLACE ;
578
579 SVD_PLACE:
580
581 #if 0
582 vmat = (double *)calloc( sizeof(double),n*n ); /* right singular vectors */
583 #endif
584 umat = (double *)calloc( sizeof(double),m*n ); /* left singular vectors */
585
586 if( umat == NULL ){ /* 29 Dec 2008 */
587 ERROR_message("mri_matrix_psinv: can't malloc umat workspace!") ;
588 free(vmat) ; free(amat) ; free(xfac) ; RETURN(NULL) ;
589 }
590
591 sval = (double *)calloc( sizeof(double),n ); /* singular values */
592
593 /* compute SVD of scaled matrix */
594
595 STATUS("SVD") ;
596 svd_double( m , n , amat , sval , umat , vmat ) ;
597
598 free(amat) ; /* done with this */
599
600 /* find largest singular value */
601
602 smax = sval[0] ;
603 for( ii=1 ; ii < n ; ii++ ) if( sval[ii] > smax ) smax = sval[ii] ;
604
605 if( smax <= 0.0 ){ /* this is bad */
606 static int first = 1 ;
607 #pragma omp critical (STDERR)
608 { if( first ) ERROR_message("SVD fails in mri_matrix_psinv()!\n"); }
609 free(xfac); free(sval); first = 0;
610 free(vmat); free(umat); RETURN( NULL);
611 }
612
613 for( ii=0 ; ii < n ; ii++ )
614 if( sval[ii] < 0.0 ) sval[ii] = 0.0 ; /* should not happen */
615
616 if( verb_rpn ){
617 fprintf(stderr," singular values:") ;
618 for( ii=0 ; ii < n ; ii++ )
619 fprintf(stderr," %g",sval[ii]) ;
620 fprintf(stderr,"\n") ;
621 }
622
623 /* "reciprocals" of singular values: 1/s is actually s/(s^2+del) */
624
625 del = PSINV_EPS * smax*smax ; if( del < alp ) del = alp ;
626 for( ii=0 ; ii < n ; ii++ )
627 sval[ii] = sval[ii] / ( sval[ii]*sval[ii] + del ) ;
628
629 /* create pseudo-inverse */
630
631 imp = mri_new( n , m , MRI_float ) ; /* recall that m > n */
632 pmat = MRI_FLOAT_PTR(imp) ;
633
634 STATUS("psinv from SVD") ;
635 for( ii=0 ; ii < n ; ii++ ){
636 for( jj=0 ; jj < m ; jj++ ){
637 sum = 0.0 ;
638 for( kk=0 ; kk < n ; kk++ ) sum += sval[kk] * V(ii,kk) * U(jj,kk) ;
639 P(ii,jj) = (float)sum ;
640 }
641 }
642 free(sval); free(vmat); free(umat);
643
644 RESCALE_PLACE:
645 /** from either method, must now rescale rows from norming */
646
647 STATUS("rescale") ;
648 for( ii=0 ; ii < n ; ii++ ){
649 for( jj=0 ; jj < m ; jj++ ) P(ii,jj) *= xfac[ii] ;
650 }
651 free(xfac);
652
653 /* rescale cols for weight? */
654
655 if( wt != NULL ){
656 for( ii=0 ; ii < m ; ii++ ){
657 ww = wt[ii] ;
658 if( ww > 0.0 ) for( jj=0 ; jj < n ; jj++ ) P(jj,ii) *= ww ;
659 }
660 }
661
662 RETURN( imp );
663 }
664
665 /*----------------------------------------------------------------------------*/
666 /*! The output matrix is the orthogonal projection onto the linear space
667 spanned by the columns of the input imc. If pout != 0, then instead
668 it is the orthogonal projection onto the complement of this space.
669 If the input is NxM, the output is NxN. Note that the matrix output
670 by this function will be symmetric. [10 Apr 2006]
671 ------------------------------------------------------------------------------*/
672
mri_matrix_ortproj(MRI_IMAGE * imc,int pout)673 MRI_IMAGE * mri_matrix_ortproj( MRI_IMAGE *imc , int pout )
674 {
675 MRI_IMAGE *imp , *imt ;
676
677 ENTRY("mri_matrix_ortproj") ;
678
679 if( imc == NULL || imc->kind != MRI_float ) RETURN( NULL );
680
681 imp = mri_matrix_psinv( imc , NULL , 0.0f ) ; /* inv[C'C] C' */
682 if( imp == NULL ) RETURN(NULL) ;
683 imt = mri_matrix_mult( imc , imp ) ; /* C inv[C'C] C' */
684 mri_free(imp) ;
685
686 if( pout ){ /* I - C inv[C'C] C' */
687 int nn , nq , ii ; float *tar ;
688 nn = imt->nx ; nq = nn*nn ; tar = MRI_FLOAT_PTR(imt) ;
689 for( ii=0 ; ii < nq ; ii+=(nn+1) ) tar[ii] -= 1.0f ;
690 for( ii=0 ; ii < nq ; ii++ ) tar[ii] = -tar[ii] ;
691 }
692
693 RETURN(imt) ;
694 }
695
696 /*----------------------------------------------------------------------------*/
697 /*! Return both the pseudo-inverse and the orthogonal projection.
698 If [C] is N x M, then the psinv is M x N and the ortproj is N x N. */
699
mri_matrix_psinv_ortproj(MRI_IMAGE * imc,int pout)700 MRI_IMARR * mri_matrix_psinv_ortproj( MRI_IMAGE *imc , int pout )
701 {
702 MRI_IMARR *imar ; MRI_IMAGE *imp , *imt ;
703
704 ENTRY("mri_matrix_psinv_ortproj") ;
705
706 if( imc == NULL || imc->kind != MRI_float ) RETURN( NULL );
707
708 imp = mri_matrix_psinv( imc , NULL , 0.0f ) ; /* inv[C'C] C' */
709 if( imp == NULL ) RETURN(NULL) ;
710 imt = mri_matrix_mult( imc , imp ) ; /* C inv[C'C] C' */
711
712 if( pout ){ /* I - C inv[C'C] C' */
713 int nn , nq , ii ; float *tar ;
714 nn = imt->nx ; nq = nn*nn ; tar = MRI_FLOAT_PTR(imt) ;
715 for( ii=0 ; ii < nq ; ii+=(nn+1) ) tar[ii] -= 1.0f ; /* diagonal */
716 for( ii=0 ; ii < nq ; ii++ ) tar[ii] = -tar[ii] ;
717 }
718
719 INIT_IMARR(imar) ; ADDTO_IMARR(imar,imp) ; ADDTO_IMARR(imar,imt) ; RETURN(imar) ;
720 }
721
722 /*----------------------------------------------------------------------------*/
723 /*! Return both inv[C'C]C' and inv[C'C] -- RWCox -- 19 May 2010 */
724
mri_matrix_psinv_pair(MRI_IMAGE * imc,float alpha)725 MRI_IMARR * mri_matrix_psinv_pair( MRI_IMAGE *imc , float alpha )
726 {
727 float *rmat ;
728 int m , n , ii,jj,kk ;
729 double *amat , *umat , *vmat , *sval , *xfac , smax,del,ww , alp ;
730 MRI_IMAGE *imp=NULL , *imq=NULL ; float *pmat,*qmat ; MRI_IMARR *imar ;
731 register double sum ;
732
733 ENTRY("mri_matrix_psinv_pair") ;
734
735 if( imc == NULL || imc->kind != MRI_float ) RETURN(NULL);
736 m = imc->nx ; /* number of rows in input */
737 n = imc->ny ; /* number of columns */
738 if( m < 1 || n < 1 ) RETURN(NULL) ; /* bad users should be shot! */
739
740 rmat = MRI_FLOAT_PTR(imc) ;
741
742 /* deal with a single vector (of length m) */
743
744 if( n == 1 ){
745 for( sum=0.0,ii=0 ; ii < m ; ii++ ) sum += rmat[ii]*rmat[ii] ;
746 imp = mri_new( 1 , m , MRI_float ) ;
747 imq = mri_new( 1 , 1 , MRI_float ) ; /* a very boring image */
748 if( sum > 0.0 ){
749 sum = 1.0 / sum ; pmat = MRI_FLOAT_PTR(imp) ;
750 for( ii=0 ; ii < m ; ii++ ) pmat[ii] = sum * rmat[ii] ;
751 qmat = MRI_FLOAT_PTR(imq) ; qmat[0] = sum ;
752 }
753 INIT_IMARR(imar); ADDTO_IMARR(imar,imp); ADDTO_IMARR(imar,imq); RETURN(imar);
754 }
755
756 /* OK, have a real matrix to handle here (at least 2 columns) */
757
758 amat = (double *)calloc( sizeof(double),m*n ) ; /* input matrix */
759 xfac = (double *)calloc( sizeof(double),n ) ; /* column norms of [a] */
760
761 if( amat == NULL || xfac == NULL ){
762 ERROR_message("mri_matrix_psinv_pair: can't malloc matrix workspace!") ;
763 if( amat != NULL ) free(amat) ;
764 if( xfac != NULL ) free(xfac) ;
765 RETURN(NULL) ;
766 }
767
768 #undef PSINV_EPS
769 #define PSINV_EPS 1.e-12
770
771 alp = (alpha <= 0.0f) ? PSINV_EPS : (double)alpha ;
772
773 #undef R
774 #undef A
775 #undef P
776 #undef U
777 #undef V
778 #undef Q
779 #define R(i,j) rmat[(i)+(j)*m] /* i=0..m-1 , j=0..n-1 */
780 #define A(i,j) amat[(i)+(j)*m] /* i=0..m-1 , j=0..n-1 */
781 #define P(i,j) pmat[(i)+(j)*n] /* i=0..n-1 , j=0..m-1 */
782 #define U(i,j) umat[(i)+(j)*m]
783 #define V(i,j) vmat[(i)+(j)*n]
784 #define Q(i,j) qmat[(i)+(j)*n] /* i=0..n-1 , j=0..n-1 */
785
786 /* copy input matrix into amat */
787
788 STATUS("copy matrix") ;
789 for( ii=0 ; ii < m ; ii++ )
790 for( jj=0 ; jj < n ; jj++ ) A(ii,jj) = R(ii,jj) ;
791
792 /* scale each column to have norm 1 */
793
794 STATUS("scale matrix") ;
795 for( jj=0 ; jj < n ; jj++ ){
796 for( sum=0.0,ii=0 ; ii < m ; ii++ ) sum += A(ii,jj)*A(ii,jj) ;
797 if( sum > 0.0 ) sum = 1.0/sqrt(sum) ;
798 xfac[jj] = sum ;
799 if( sum > 0.0 ){ for( ii=0 ; ii < m ; ii++ ) A(ii,jj) *= sum ; }
800 }
801
802 /*** real computations follow, via SVD ***/
803
804 vmat = (double *)calloc( sizeof(double),n*n );
805
806 if( vmat == NULL ){ /* 29 Dec 2008 */
807 ERROR_message("mri_matrix_psinv_pair: can't malloc vmat workspace!") ;
808 free(amat) ; free(xfac) ; RETURN(NULL) ;
809 }
810
811 umat = (double *)calloc( sizeof(double),m*n ); /* left singular vectors */
812
813 if( umat == NULL ){ /* 29 Dec 2008 */
814 ERROR_message("mri_matrix_psinv_pair: can't malloc umat workspace!") ;
815 free(vmat) ; free(amat) ; free(xfac) ; RETURN(NULL) ;
816 }
817
818 sval = (double *)calloc( sizeof(double),n ); /* singular values */
819
820 /* compute SVD of scaled matrix */
821
822 STATUS("SVD") ;
823 svd_double( m , n , amat , sval , umat , vmat ) ;
824
825 free(amat) ; /* done with this */
826
827 /* find largest singular value */
828
829 smax = sval[0] ;
830 for( ii=1 ; ii < n ; ii++ ) if( sval[ii] > smax ) smax = sval[ii] ;
831
832 if( smax <= 0.0 ){ /* this is bad */
833 static int first = 1 ;
834 #pragma omp critical (STDERR)
835 { if( first ) ERROR_message("SVD fails in mri_matrix_psinv_pair()!\n"); }
836 free(xfac); free(sval); first = 0;
837 free(vmat); free(umat); RETURN(NULL);
838 }
839
840 for( ii=0 ; ii < n ; ii++ )
841 if( sval[ii] < 0.0 ) sval[ii] = 0.0 ; /* should not happen */
842
843 /* "reciprocals" of singular values: 1/s is actually s/(s^2+del) */
844
845 del = PSINV_EPS * smax*smax ; if( del < alp ) del = alp ;
846 for( ii=0 ; ii < n ; ii++ )
847 sval[ii] = sval[ii] / ( sval[ii]*sval[ii] + del ) ;
848
849 /* create pseudo-inverse */
850
851 STATUS("psinv from SVD") ;
852 imp = mri_new( n , m , MRI_float ) ; /* recall that m > n */
853 pmat = MRI_FLOAT_PTR(imp) ;
854
855 for( ii=0 ; ii < n ; ii++ ){
856 for( jj=0 ; jj < m ; jj++ ){
857 sum = 0.0 ;
858 for( kk=0 ; kk < n ; kk++ ) sum += sval[kk] * V(ii,kk) * U(jj,kk) ;
859 P(ii,jj) = (float)(sum*xfac[ii]) ;
860 }
861 }
862
863 STATUS("inv[XtX] from SVD") ;
864 imq = mri_new( n , n , MRI_float ) ;
865 qmat = MRI_FLOAT_PTR(imq) ;
866
867 for( kk=0 ; kk < n ; kk++ ) sval[kk] *= sval[kk] ; /* 1/s^2 */
868
869 for( ii=0 ; ii < n ; ii++ ){
870 for( jj=0 ; jj <= ii ; jj++ ){
871 sum = 0.0 ;
872 for( kk=0 ; kk < n ; kk++ )
873 sum += sval[kk] * V(ii,kk) * V(jj,kk) ;
874 Q(ii,jj) = (float)(sum * xfac[ii] * xfac[jj]) ;
875 if( jj < ii ) Q(jj,ii) = Q(ii,jj) ;
876 }
877 }
878
879 free(sval); free(vmat); free(umat); free(xfac) ;
880
881 INIT_IMARR(imar); ADDTO_IMARR(imar,imp); ADDTO_IMARR(imar,imq); RETURN(imar);
882 }
883
884 /*---------------------------------------------------------------------------*/
885 /*! Return the sqrt(eigenvalues) of NxN matrix [C'C], scaled to diagonal 1.
886 The output points to a 1D image as vector of floats, of length N=C->ny.
887 This image should be mri_free()-ed when you are done with it.
888 -----------------------------------------------------------------------------*/
889
mri_matrix_singvals(MRI_IMAGE * imc)890 MRI_IMAGE * mri_matrix_singvals( MRI_IMAGE *imc ) /* 24 May 2010 */
891 {
892 int i,j,k , M , N ;
893 double *a , *e , sum ;
894 float *rmat ; MRI_IMAGE *ime ;
895
896 ENTRY("mri_matrix_singvals") ;
897
898 if( imc == NULL || imc->kind != MRI_float ) RETURN(NULL);
899
900 M = imc->nx ; N = imc->ny ;
901 a = (double *)malloc( sizeof(double)*N*N ) ;
902 e = (double *)malloc( sizeof(double)*N ) ;
903
904 #undef R
905 #define R(i,j) rmat[(i)+(j)*M]
906 rmat = MRI_FLOAT_PTR(imc) ;
907
908 for( i=0 ; i < N ; i++ ){
909 for( j=0 ; j <= i ; j++ ){
910 sum = 0.0 ;
911 for( k=0 ; k < M ; k++ ) sum += R(k,i)*R(k,j) ;
912 a[j+N*i] = sum ; if( j < i ) a[i+N*j] = sum ;
913 }
914 }
915
916 for( i=0 ; i < N ; i++ ){
917 if( a[i+N*i] > 0.0 ) e[i] = 1.0 / sqrt(a[i+N*i]) ;
918 else e[i] = 1.0 ;
919 }
920 for( i=0 ; i < N ; i++ ){
921 for( j=0 ; j < N ; j++ ) a[j+N*i] *= e[i]*e[j] ;
922 }
923
924 symeigval_double( N , a , e ) ; free(a) ;
925 ime = mri_new( N , 1 , MRI_float ) ; rmat = MRI_FLOAT_PTR(ime) ;
926 for( j=0 ; j < N ; j++ )
927 rmat[j] = (float) ( (e[j] <= 0.0) ? 0.0 : sqrt(e[j]) ) ;
928 free(e) ; RETURN(ime) ;
929 }
930
931 /*----------------------------------------------------------------------------*/
932 /*! imt = time series to detrend in place (can be more than 1 vector in here)
933 imq = array of orts to remove (imq->nx must equal imt->nx)
934 imp = mri_matrix_psinv(imq) (imp is imq->ny X imq->nx)
935 *//*--------------------------------------------------------------------------*/
936
mri_matrix_detrend(MRI_IMAGE * imt,MRI_IMAGE * imq,MRI_IMAGE * imp)937 void mri_matrix_detrend( MRI_IMAGE *imt , MRI_IMAGE *imq , MRI_IMAGE *imp )
938 {
939 int nlen , nort , nvec , kk,jj,iv ;
940 float *par , *qar , *tar , *xar , *rar , *pt , *qt , xt ;
941
942 ENTRY("mri_matrix_detrend") ;
943
944 if( imt == NULL || imq == NULL || imp == NULL ) EXRETURN ;
945
946 nlen = imt->nx ; if( nlen != imq->nx ) EXRETURN ;
947 nort = imq->ny ; if( imp->nx != nort || imp->ny != nlen ) EXRETURN ;
948 nvec = imt->ny ;
949
950 tar = MRI_FLOAT_PTR(imt) ; /* nlen X nvec */
951 qar = MRI_FLOAT_PTR(imq) ; /* nlen X nort */
952 par = MRI_FLOAT_PTR(imp) ; /* nort X nlen */
953 rar = (float *)malloc(sizeof(float)*nort) ; /* workspace */
954
955 for( iv=0 ; iv < nvec ; iv++ ){ /* [xar] <- [xar] - [Q][P][xar] */
956 xar = tar + iv*nlen ;
957 for( jj=0 ; jj < nort ; jj++ ) rar[jj] = 0.0f ;
958 for( kk=0 ; kk < nlen ; kk++ ){ /* [rar] <- [P][xar] */
959 pt = par + kk*nort ; xt = xar[kk] ;
960 for( jj=0 ; jj < nort ;jj++ ) rar[jj] += pt[jj]*xt ;
961 }
962 for( jj=0 ; jj < nort ; jj++ ){ /* [xar] <- [Q][rar] */
963 qt = qar + jj*nlen ; xt = rar[jj] ;
964 for( kk=0 ; kk < nlen ; kk++ ) xar[kk] -= qt[kk]*xt ;
965 }
966 }
967
968 free(rar) ; EXRETURN ;
969 }
970
971 /*----------------------------------------------------------------------------*/
972 /*! Return the average absolute value of the entries of the matrix. */
973
mri_matrix_size(MRI_IMAGE * imc)974 float mri_matrix_size( MRI_IMAGE *imc ) /* 30 Jul 2007 */
975 {
976 int nxy , ii ;
977 float sum , *car ;
978
979 if( imc == NULL || imc->kind != MRI_float ) return(-1.0f) ;
980 nxy = imc->nx * imc->ny ; car = MRI_FLOAT_PTR(imc) ;
981 sum = 0.0f ;
982 for( ii=0 ; ii < nxy ; ii++ ) sum += fabs(car[ii]) ;
983 sum /= nxy ; return sum ;
984 }
985
986 /*----------------------------------------------------------------------------*/
987 /*! Compute the square root of a matrix, iteratively. */
988
mri_matrix_sqrt(MRI_IMAGE * imc)989 MRI_IMAGE * mri_matrix_sqrt( MRI_IMAGE *imc ) /* 30 Jul 2007 */
990 {
991 float gam , fa,fb , csiz ;
992 int nn , ite , ii ;
993 MRI_IMAGE *imy , *imz , *imyinv,*imzinv , *tim ;
994 float *yar , *zar , *car ;
995
996 if( imc == NULL || imc->kind != MRI_float ) RETURN(NULL) ;
997 nn = imc->nx ; if( nn != imc->ny ) RETURN(NULL) ;
998 if( nn == 1 ){
999 car = MRI_FLOAT_PTR(imc); if( car[0] < 0.0f ) RETURN(NULL) ;
1000 imz = mri_new(1,1,MRI_float); zar = MRI_FLOAT_PTR(imz);
1001 zar[0] = sqrt(car[0]) ; RETURN(imz) ;
1002 }
1003
1004 imz = mri_new(nn,nn,MRI_float) ; zar = MRI_FLOAT_PTR(imz) ;
1005 csiz = mri_matrix_size(imc) ; if( csiz <= 0.0f ) RETURN(imz) ;
1006 imy = mri_copy(imc) ;
1007 for( ii=0 ; ii < nn ; ii++ ) zar[ii+ii*nn] = 1.0f ; /* identity */
1008
1009 /** start iterations: usually 3-5 is enough **/
1010
1011 for( ite=0 ; ite < 50 ; ite++ ){
1012 imyinv = mri_matrix_psinv( imy , NULL , 0.0f ) ;
1013 if( imyinv == NULL ){
1014 ERROR_message("mri_matrix_sqrt() fails at ite=%d",ite);
1015 mri_free(imz); mri_free(imy); RETURN(NULL);
1016 }
1017 if( ite == 0 ) imzinv = mri_copy(imz) ;
1018 else imzinv = mri_matrix_psinv( imz , NULL , 0.0f ) ;
1019 if( imzinv == NULL ){
1020 ERROR_message("mri_matrix_sqrt() fails at ite=%d",ite);
1021 mri_free(imz); mri_free(imy); mri_free(imyinv); RETURN(NULL);
1022 }
1023 gam = 1.0f ; /* someday, could put a scaling calculation in for gam */
1024 fa = 0.5f*gam ; fb = 0.5f/gam ;
1025 tim = mri_matrix_sadd( fa,imy , fb,imzinv ) ; mri_free(imy) ; imy = tim ;
1026 tim = mri_matrix_sadd( fa,imz , fb,imyinv ) ; mri_free(imz) ; imz = tim ;
1027 mri_free(imyinv) ; mri_free(imzinv) ;
1028 if( ite > 2 ){ /* check for convergence: see if imy*imy-imc is small */
1029 imyinv = mri_matrix_mult( imy , imy ) ;
1030 tim = mri_matrix_sadd( 1.0f,imyinv , -1.0f,imc ) ; mri_free(imyinv) ;
1031 fa = mri_matrix_size(tim) ; mri_free(tim) ;
1032 if( fa <= 5.e-5*csiz ){ mri_free(imz); RETURN(imy); }
1033 }
1034 }
1035 #if 0
1036 ERROR_message("mri_matrix_sqrt() fails to converge: err=%g",fa/csiz);
1037 #endif
1038 mri_free(imz) ; mri_free(imy) ; RETURN(NULL) ;
1039 }
1040
1041 /*----------------------------------------------------------------------------*/
1042 /*! The first principal component (singular) vector of a matrix [16 Jul 2009] */
1043
mri_matrix_pcvector(MRI_IMAGE * imx)1044 MRI_IMAGE * mri_matrix_pcvector( MRI_IMAGE *imx )
1045 {
1046 int jj ; MRI_IMAGE *imu ;
1047
1048 ENTRY("mri_matrix_pcvector") ;
1049
1050 if( imx == NULL || imx->kind != MRI_float ) RETURN(NULL) ;
1051
1052 imu = mri_new( imx->nx , 1 , MRI_float ) ;
1053
1054 jj = first_principal_vectors( imx->nx , imx->ny , MRI_FLOAT_PTR(imx) ,
1055 1 , NULL , MRI_FLOAT_PTR(imu) ) ;
1056
1057 if( jj <= 0 ){ mri_free(imu) ; imu = NULL ; }
1058
1059 RETURN(imu) ;
1060 }
1061
1062 /*----------------------------------------------------------------------------*/
1063 /*! Square root of a mat44 struct (3x3 matrix + 3 vect). */
1064
THD_mat44_sqrt(mat44 A)1065 mat44 THD_mat44_sqrt( mat44 A ) /* 30 Jul 2007 */
1066 {
1067 MRI_IMAGE *imc, *imx ; float *far ; mat44 X ;
1068
1069 imc = mri_new(4,4,MRI_float) ; far = MRI_FLOAT_PTR(imc) ;
1070 UNLOAD_MAT44_AR( A , far ) ;
1071 far[12] = far[13] = far[14] = 0.0f ; far[15] = 1.0f ;
1072 imx = mri_matrix_sqrt(imc) ; mri_free(imc) ;
1073 if( imx == NULL ){
1074 LOAD_DIAG_MAT44(X,0,0,0) ; /* error! */
1075 } else {
1076 far = MRI_FLOAT_PTR(imx) ;
1077 LOAD_MAT44_AR(X,far) ; mri_free(imx) ;
1078 }
1079 return X ;
1080 }
1081
1082 /*----------------------------------------------------------------------------*/
1083 /*! Legendre polynomial of non-negative order m evaluated at x;
1084 x may be any real number, but the main use is for -1 <= x <= 1 (duuuh).
1085 ------------------------------------------------------------------------------*/
1086
Plegendre(double x,int m)1087 double Plegendre( double x , int m )
1088 {
1089 double pk=0.0, pkm1, pkm2 ; int k ; /* for the recurrence, when m > 20 */
1090
1091 if( m < 0 ) return 1.0 ; /* bad input */
1092
1093 switch( m ){ /*** direct formulas for P_m(x) for m=0..20 ***/
1094 case 0: return 1.0 ; /* that was easy */
1095 case 1: return x ; /* also pretty easy */
1096 case 2: return (3.0*x*x-1.0)/2.0 ; /* now it gets harder */
1097 case 3: return (5.0*x*x-3.0)*x/2.0 ;
1098 case 4: return ((35.0*x*x-30.0)*x*x+3.0)/8.0 ;
1099 case 5: return ((63.0*x*x-70.0)*x*x+15.0)*x/8.0 ;
1100 case 6: return (((231.0*x*x-315.0)*x*x+105.0)*x*x-5.0)/16.0 ;
1101 case 7: return (((429.0*x*x-693.0)*x*x+315.0)*x*x-35.0)*x/16.0 ;
1102 case 8: return ((((6435.0*x*x-12012.0)*x*x+6930.0)*x*x-1260.0)*x*x+35.0)/128.0;
1103
1104 /** 07 Feb 2005: this part generated by Maple, then hand massaged **/
1105
1106 case 9:
1107 return (0.24609375e1 + (-0.3609375e2 + (0.140765625e3 + (-0.20109375e3
1108 + 0.949609375e2 * x * x) * x * x) * x * x) * x * x) * x;
1109
1110 case 10:
1111 return -0.24609375e0 + (0.1353515625e2 + (-0.1173046875e3 +
1112 (0.3519140625e3 + (-0.42732421875e3 + 0.18042578125e3 * x * x)
1113 * x * x) * x * x) * x * x) * x * x;
1114
1115 case 11:
1116 return (-0.270703125e1 + (0.5865234375e2 + (-0.3519140625e3 +
1117 (0.8546484375e3 + (-0.90212890625e3 + 0.34444921875e3 * x * x)
1118 * x * x) * x * x) * x * x) * x * x) * x;
1119
1120 case 12:
1121 return 0.2255859375e0 + (-0.17595703125e2 + (0.2199462890625e3 +
1122 (-0.99708984375e3 + (0.20297900390625e4 + (-0.1894470703125e4
1123 + 0.6601943359375e3 * x * x) * x * x) * x * x) * x * x) * x * x)
1124 * x * x;
1125
1126 case 13:
1127 return (0.29326171875e1 + (-0.87978515625e2 + (0.7478173828125e3 +
1128 (-0.270638671875e4 + (0.47361767578125e4 + (-0.3961166015625e4
1129 + 0.12696044921875e4 * x * x) * x * x) * x * x) * x * x) * x * x)
1130 * x * x) * x;
1131
1132 case 14:
1133 return -0.20947265625e0 + (0.2199462890625e2 + (-0.37390869140625e3 +
1134 (0.236808837890625e4 + (-0.710426513671875e4 +
1135 (0.1089320654296875e5 + (-0.825242919921875e4 +
1136 0.244852294921875e4 * x * x) * x * x) * x * x) * x * x) * x * x)
1137 * x * x) * x * x;
1138
1139 case 15:
1140 return (-0.314208984375e1 + (0.12463623046875e3 + (-0.142085302734375e4
1141 + (0.710426513671875e4 + (-0.1815534423828125e5 +
1142 (0.2475728759765625e5 + (-0.1713966064453125e5 +
1143 0.473381103515625e4 * x * x) * x * x) * x * x) * x * x)
1144 * x * x) * x * x) * x * x) * x;
1145
1146 case 16:
1147 return 0.196380615234375e0 + (-0.26707763671875e2 + (0.5920220947265625e3
1148 + (-0.4972985595703125e4 + (0.2042476226806641e5 +
1149 (-0.4538836059570312e5 + (0.5570389709472656e5 +
1150 (-0.3550358276367188e5 + 0.9171758880615234e4 * x * x) * x * x)
1151 * x * x) * x * x) * x * x) * x * x) * x * x) * x * x;
1152
1153 case 17:
1154 return (0.3338470458984375e1 + (-0.1691491699218750e3 +
1155 (0.2486492797851562e4 + (-0.1633980981445312e5 +
1156 (0.5673545074462891e5 + (-0.1114077941894531e6 +
1157 (0.1242625396728516e6 + (-0.7337407104492188e5 +
1158 0.1780400253295898e5 * x * x) * x * x) * x * x) * x * x)
1159 * x * x) * x * x) * x * x) * x * x) * x;
1160
1161 case 18:
1162 return -0.1854705810546875e0 + (0.3171546936035156e2 +
1163 (-0.8880331420898438e3 + (0.9531555725097656e4 +
1164 (-0.5106190567016602e5 + (0.1531857170104980e6 +
1165 (-0.2692355026245117e6 + (0.2751527664184570e6 +
1166 (-0.1513340215301514e6 + 0.34618893814086910e5 * x * x) * x * x)
1167 * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x * x;
1168
1169 case 19:
1170 return (-0.3523941040039062e1 + (0.2220082855224609e3 +
1171 (-0.4084952453613281e4 + (0.3404127044677734e5 +
1172 (-0.1531857170104980e6 + (0.4038532539367676e6 +
1173 (-0.6420231216430664e6 + (0.6053360861206055e6 +
1174 (-0.3115700443267822e6 + 0.67415740585327150e5 * x * x) * x * x)
1175 * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x;
1176
1177 case 20:
1178 return 0.1761970520019531e0 + (-0.3700138092041016e2 +
1179 (0.1276547641754150e4 + (-0.1702063522338867e5 +
1180 (0.1148892877578735e6 + (-0.4442385793304443e6 +
1181 (0.1043287572669983e7 + (-0.1513340215301514e7 +
1182 (0.1324172688388824e7 + (-0.6404495355606079e6 +
1183 0.1314606941413879e6 * x * x) * x * x) * x * x) * x * x) * x * x)
1184 * x * x) * x * x) * x * x) * x * x) * x * x;
1185 }
1186
1187 /**----- if here, m > 20 ==> use recurrence relation ----**/
1188
1189 pkm2 = Plegendre( x , 19 ) ; /* recursion to start things off! */
1190 pkm1 = Plegendre( x , 20 ) ;
1191 for( k=21 ; k <= m ; k++ , pkm2=pkm1 , pkm1=pk )
1192 pk = ((2.0*k-1.0)*x*pkm1 - (k-1.0)*pkm2)/k ;
1193 return pk ;
1194 }
1195
1196 /******************************************************************************/
1197 /**** Matrix RPN calculator ***************************************************/
1198 /******************************************************************************/
1199
1200 static MRI_IMARR *matar = NULL ; /* list of named matrices */
1201
1202 /*----------------------------------------------------------------------------*/
1203
matrix_name_lookup(char * nam)1204 static int matrix_name_lookup( char *nam )
1205 {
1206 int ii ;
1207 MRI_IMAGE *im ;
1208
1209 if( nam == NULL || matar == NULL ) return -1 ;
1210
1211 for( ii=0 ; ii < IMARR_COUNT(matar) ; ii++ ){
1212 im = IMARR_SUBIM(matar,ii) ;
1213 if( im != NULL && strcmp(nam,im->name) == 0 ) return ii;
1214 }
1215
1216 return -1 ;
1217 }
1218
1219 /*----------------------------------------------------------------------------*/
1220
matrix_name_assign(char * nam,MRI_IMAGE * ima)1221 static void matrix_name_assign( char *nam , MRI_IMAGE *ima )
1222 {
1223 int ii ; MRI_IMAGE *imb ;
1224
1225 if( nam == NULL || *nam == '\0' || ima == NULL ) return ;
1226
1227 if( matar == NULL ) INIT_IMARR(matar) ;
1228
1229 ii = matrix_name_lookup( nam ) ;
1230 imb = mri_to_float(ima) ; mri_add_name(nam,imb) ;
1231 #pragma omp critical (MATAR)
1232 { if( ii < 0 ){
1233 ADDTO_IMARR(matar,imb) ;
1234 } else {
1235 mri_free( IMARR_SUBIM(matar,ii) ) ;
1236 IMARR_SUBIM(matar,ii) = imb ;
1237 }
1238 }
1239 return ;
1240 }
1241
1242 /*----------------------------------------------------------------------------*/
1243
char_check(char ccc)1244 static int char_check( char ccc )
1245 {
1246 return ( ccc == '&' || ccc == '@' || ccc == '%' ) ;
1247 }
1248
command_check(char * str,char * cmd)1249 static int command_check( char *str , char *cmd )
1250 {
1251 if( str == NULL || cmd == NULL || *str == '\0' || *cmd == '\0' ) return 0 ;
1252 if( char_check(str[0]) == 0 ) return 0 ;
1253 return (strncasecmp(str+1,cmd,strlen(cmd)) == 0) ;
1254 }
1255
1256 /*----------------------------------------------------------------------------*/
1257
1258 #undef ERREX
1259 #define ERREX(sss) \
1260 do{ ERROR_message("mri_matrix_evalrpn('%s') at '%s': %s" , expr,cmd,sss ); \
1261 DESTROY_IMARR(imstk); NI_delete_str_array(sar); RETURN(NULL); \
1262 } while(0)
1263
1264 #undef DECODE_VALUE
1265 #define DECODE_VALUE(ccc,vvv) \
1266 do{ float val=-666.0f ; \
1267 if( isdigit(*(ccc)) ){ \
1268 sscanf((ccc),"%f",&val) ; \
1269 } else if( char_check(*(ccc)) && toupper(*((ccc)+1)) == 'R' ){ \
1270 if( nstk > 0 ) val = IMARR_SUBIM(imstk,nstk-1)->nx ; \
1271 else ERROR_message("mri_matrix_evalrpn: bad (&R)") ; \
1272 } else if( char_check(*(ccc)) && toupper(*((ccc)+1)) == 'C' ){ \
1273 if( nstk > 0 ) val = IMARR_SUBIM(imstk,nstk-1)->ny ; \
1274 else ERROR_message("mri_matrix_evalrpn: bad (&C)") ; \
1275 } else if( *(ccc) == '=' && isalpha(*((ccc)+1)) ){ \
1276 char *xx=strdup((ccc)+1) , *pp=strchr(xx,')') ; \
1277 int qq = matrix_name_lookup(xx) ; /* swap these lines */ \
1278 if( pp != NULL ) *pp = '\0' ; /* 18 Apr 2006 [rickr] */ \
1279 if( qq >= 0 ){ \
1280 MRI_IMAGE *iq = IMARR_SUBIM(matar,qq) ; \
1281 float *qar = MRI_FLOAT_PTR(iq) ; \
1282 val = qar[0] ; \
1283 } \
1284 else ERROR_message("mri_matrix_evalrpn: bad (=%s)",xx); \
1285 } \
1286 (vvv) = val ; \
1287 } while(0)
1288
1289 /*----------------------------------------------------------------------------*/
1290 /*! Evaluate a space delimited RPN expression to evaluate matrices:
1291 - The operations allowed are described in the output of
1292 function mri_matrix_evalrpn_help().
1293 - The return value is the top matrix at the end of the expression.
1294 - The rest of the stack is discarded, but named matrices are preserved
1295 in a static array between calls.
1296 - If NULL is returned, some error occurred.
1297 ------------------------------------------------------------------------------*/
1298
mri_matrix_evalrpn(char * expr)1299 MRI_IMAGE * mri_matrix_evalrpn( char *expr )
1300 {
1301 NI_str_array *sar ;
1302 char *cmd , mess[512] ;
1303 MRI_IMARR *imstk ;
1304 int nstk , ii , ss ;
1305 MRI_IMAGE *ima, *imb, *imc ;
1306 float *car ;
1307
1308 ENTRY("mri_matrix_evalrpn") ;
1309
1310 /** break string into sub-strings, delimited by whitespace **/
1311
1312 sar = NI_decode_string_list( expr , "~" ) ;
1313 if( sar == NULL ) RETURN(NULL) ;
1314
1315 INIT_IMARR(imstk) ; /* initialize stack */
1316 mri_matrix_psinv_svd(1) ; /* always use SVD for inversion */
1317
1318 /** loop thru and process commands **/
1319
1320 if(verb_rpn)fprintf(stderr,"++ mri_matrix_evalrpn('%s')\n",expr) ;
1321
1322 for( ss=0 ; ss < sar->num ; ss++ ){
1323
1324 cmd = sar->str[ss] ;
1325 nstk = IMARR_COUNT(imstk) ;
1326
1327 if(verb_rpn){
1328 int ss , na,nb ; float *amat ;
1329 fprintf(stderr," + nstk=%d cmd='%s'\n",nstk,cmd) ;
1330 for( ss=nstk-1 ; ss >= 0 ; ss-- ){
1331 ima = IMARR_SUBIM(imstk,ss) ;
1332 fprintf(stderr," [%d] %d X %d",ss,ima->nx,ima->ny) ;
1333 if( ima->nx == 1 && ima->ny == 1 ){
1334 amat = MRI_FLOAT_PTR(ima) ;
1335 fprintf(stderr," = %g",amat[0]) ;
1336 }
1337 fprintf(stderr,"\n") ;
1338 }
1339 }
1340
1341 if( *cmd == '\0' ) continue ; /* WTF? */
1342
1343 /** push a scalar value (1x1 matrix) onto the stack **/
1344
1345 else if( isdigit(cmd[0]) || (cmd[0]=='-' && isdigit(cmd[1])) ){
1346 imc = mri_new(1,1,MRI_float) ; car = MRI_FLOAT_PTR(imc) ;
1347 car[0] = (float)strtod(cmd,NULL) ;
1348 ADDTO_IMARR(imstk,imc) ;
1349 }
1350
1351 /** =NAME --> save top of stack as the given NAME;
1352 stack itself is unchanged **/
1353
1354 else if( cmd[0] == '=' && isalpha(cmd[1]) ){
1355 if( nstk < 1 ) ERREX("no matrix") ;
1356 ima = IMARR_SUBIM(imstk,nstk-1) ;
1357 matrix_name_assign(cmd+1,ima) ;
1358 }
1359
1360 /** clear all named matrices; stack is unchanged **/
1361
1362 else if( command_check(cmd,"clear") ){
1363 DESTROY_IMARR(matar) ;
1364 }
1365
1366 /** purge the stack; named matrices are unchanged **/
1367
1368 else if( command_check(cmd,"purge") ){ /* 05 May 2020 */
1369 TRUNCATE_IMARR(imstk,0) ;
1370 }
1371
1372 /** transpose matrix on top of stack, replacing it **/
1373
1374 else if( strcmp(cmd,"'") == 0 || command_check(cmd,"trans") || command_check(cmd,"'") ){
1375 if( nstk < 1 ) ERREX("no matrix") ;
1376 ima = IMARR_SUBIM(imstk,nstk-1) ;
1377 imc = mri_matrix_transpose(ima) ;
1378 TRUNCATE_IMARR(imstk,nstk-1) ;
1379 ADDTO_IMARR(imstk,imc) ;
1380 }
1381
1382 /** &read(filename)=NAME --> read from ASCII file; optional NAME-ing;
1383 new matrix is pushed onto top of stack **/
1384
1385 else if( command_check(cmd,"read(") ){
1386 char *buf = strdup(cmd+6) , *bp ; int heq ;
1387 for( bp=buf ; *bp != '\0' && *bp != ')' ; bp++ ) ; /*nada*/
1388 heq = (*bp == ')' && *(bp+1) == '=' && isalpha(*(bp+2)) ) ;
1389 *bp = '\0' ;
1390 imc = mri_read_1D(buf) ;
1391 if( imc == NULL ){
1392 sprintf(mess,"can't read file '%s'",buf); free(buf); ERREX(mess);
1393 }
1394 ADDTO_IMARR(imstk,imc) ;
1395 if( heq ) matrix_name_assign( bp+2 , imc ) ;
1396 free(buf) ;
1397 }
1398
1399 else if( command_check(cmd,"read4x4Xform(") ){
1400 char *buf = strdup(cmd+14) , *bp ; int heq ;
1401 for( bp=buf ; *bp != '\0' && *bp != ')' ; bp++ ) ; /*nada*/
1402 heq = (*bp == ')' && *(bp+1) == '=' && isalpha(*(bp+2)) ) ;
1403 *bp = '\0' ;
1404 imc = mri_read_4x4AffXfrm_1D(buf) ;
1405 if( imc == NULL ){
1406 sprintf(mess,"can't read file '%s'",buf); free(buf); ERREX(mess);
1407 }
1408 ADDTO_IMARR(imstk,imc) ;
1409 if( heq ) matrix_name_assign( bp+2 , imc ) ;
1410 free(buf) ;
1411 }
1412
1413 /** &write(filename) --> write top of stack to ASCII file;
1414 top of stack is unchanged; 'filename'=='-' is stdout **/
1415
1416 else if( command_check(cmd,"write(") ){
1417 char *buf = strdup(cmd+7) , *bp ;
1418 if( nstk < 1 ){ free(buf); ERREX("no matrix"); }
1419 for( bp=buf ; *bp != '\0' && *bp != ')' ; bp++ ) ; /*nada*/
1420 *bp = '\0' ;
1421 ima = IMARR_SUBIM(imstk,nstk-1) ;
1422 mri_write_1D( buf , ima ) ;
1423 free(buf) ;
1424 }
1425
1426 /** &ident(N) --> create an NxN identity matrix **/
1427
1428 else if( command_check(cmd,"ident(") ){
1429 int nn ; char *cv=cmd+7 ;
1430 DECODE_VALUE(cv,nn) ;
1431 if( nn <= 0 ) ERREX("illegal dimension") ;
1432 imc = mri_new(nn,nn,MRI_float) ; car = MRI_FLOAT_PTR(imc) ;
1433 for( ii=0 ; ii < nn ; ii++ ) car[ii+ii*nn] = 1.0f ;
1434 ADDTO_IMARR(imstk,imc) ;
1435 }
1436
1437 /** log product of singular values [05 May 2020] **/
1438
1439 else if( command_check(cmd,"svprodlog") ){
1440 if( nstk < 1 ) ERREX("no matrix") ;
1441 ima = IMARR_SUBIM(imstk,nstk-1) ;
1442 imc = mri_matrix_svprod( ima , 1 ) ;
1443 if( imc == NULL ) ERREX("can't compute") ;
1444 TRUNCATE_IMARR(imstk,nstk-1) ;
1445 ADDTO_IMARR(imstk,imc) ;
1446 }
1447
1448 /** product of singular values [05 May 2020] **/
1449
1450 else if( command_check(cmd,"svprod") ){
1451 if( nstk < 1 ) ERREX("no matrix") ;
1452 ima = IMARR_SUBIM(imstk,nstk-1) ;
1453 imc = mri_matrix_svprod( ima , 0 ) ;
1454 if( imc == NULL ) ERREX("can't compute") ;
1455 TRUNCATE_IMARR(imstk,nstk-1) ;
1456 ADDTO_IMARR(imstk,imc) ;
1457 }
1458
1459 /** vector of singular values [05 May 2020] **/
1460
1461 else if( command_check(cmd,"svals") ){
1462 if( nstk < 1 ) ERREX("no matrix") ;
1463 ima = IMARR_SUBIM(imstk,nstk-1) ;
1464 imc = mri_matrix_svals( ima ) ;
1465 if( imc == NULL ) ERREX("can't compute") ;
1466 TRUNCATE_IMARR(imstk,nstk-1) ;
1467 ADDTO_IMARR(imstk,imc) ;
1468 }
1469
1470 /** pseudo-inverse **/
1471
1472 else if( command_check(cmd,"Psinv") ){
1473 if( nstk < 1 ) ERREX("no matrix") ;
1474 ima = IMARR_SUBIM(imstk,nstk-1) ;
1475 imc = mri_matrix_psinv( ima , NULL , 0.00001f ) ;
1476 if( imc == NULL ) ERREX("can't compute") ;
1477 TRUNCATE_IMARR(imstk,nstk-1) ;
1478 ADDTO_IMARR(imstk,imc) ;
1479 }
1480
1481 /** matrix square root **/
1482
1483 else if( command_check(cmd,"Sqrt") ){
1484 if( nstk < 1 ) ERREX("no matrix") ;
1485 ima = IMARR_SUBIM(imstk,nstk-1) ;
1486 imc = mri_matrix_sqrt( ima ) ;
1487 if( imc == NULL ) ERREX("can't compute") ;
1488 TRUNCATE_IMARR(imstk,nstk-1) ;
1489 ADDTO_IMARR(imstk,imc) ;
1490 }
1491
1492 /** orthogonal projection onto column space **/
1493
1494 else if( command_check(cmd,"Pproj") ){
1495 if( nstk < 1 ) ERREX("no matrix") ;
1496 ima = IMARR_SUBIM(imstk,nstk-1) ;
1497 imc = mri_matrix_ortproj( ima , 0 ) ;
1498 if( imc == NULL ) ERREX("can't compute") ;
1499 TRUNCATE_IMARR(imstk,nstk-1) ;
1500 ADDTO_IMARR(imstk,imc) ;
1501 }
1502
1503 /** orthogonal projection onto complement of column space **/
1504
1505 else if( command_check(cmd,"Qproj") ){
1506 if( nstk < 1 ) ERREX("no matrix") ;
1507 ima = IMARR_SUBIM(imstk,nstk-1) ;
1508 imc = mri_matrix_ortproj( ima , 1 ) ;
1509 if( imc == NULL ) ERREX("can't compute") ;
1510 TRUNCATE_IMARR(imstk,nstk-1) ;
1511 ADDTO_IMARR(imstk,imc) ;
1512 }
1513
1514 /** matrix multiply **/
1515
1516 else if( strcmp(cmd,"*") == 0 || command_check(cmd,"mult") ){
1517 if( nstk < 2 ) ERREX("no matrices") ;
1518 ima = IMARR_SUBIM(imstk,nstk-2) ;
1519 imb = IMARR_SUBIM(imstk,nstk-1) ;
1520 imc = mri_matrix_mult( ima , imb ) ;
1521 if( imc == NULL ) ERREX("illegal multiply") ;
1522 TRUNCATE_IMARR(imstk,nstk-2) ;
1523 ADDTO_IMARR(imstk,imc) ;
1524 }
1525
1526 /** matrix add **/
1527
1528 else if( strcmp(cmd,"+") == 0 || command_check(cmd,"add") ){
1529 if( nstk < 2 ) ERREX("no matrices") ;
1530 ima = IMARR_SUBIM(imstk,nstk-2) ;
1531 imb = IMARR_SUBIM(imstk,nstk-1) ;
1532 imc = mri_matrix_sadd( 1.0f,ima , 1.0f,imb ) ;
1533 if( imc == NULL ) ERREX("illegal add") ;
1534 TRUNCATE_IMARR(imstk,nstk-2) ;
1535 ADDTO_IMARR(imstk,imc) ;
1536 }
1537
1538 /** matrix subtract **/
1539
1540 else if( strcmp(cmd,"-") == 0 || command_check(cmd,"sub") ){
1541 if( nstk < 2 ) ERREX("no matrices") ;
1542 ima = IMARR_SUBIM(imstk,nstk-2) ;
1543 imb = IMARR_SUBIM(imstk,nstk-1) ;
1544 imc = mri_matrix_sadd( 1.0f,ima , -1.0f,imb ) ;
1545 if( imc == NULL ) ERREX("illegal subtract") ;
1546 TRUNCATE_IMARR(imstk,nstk-2) ;
1547 ADDTO_IMARR(imstk,imc) ;
1548 }
1549
1550 /** matrix Hglue [04 Jan 2012] **/
1551
1552 else if( command_check(cmd,"Hglue") ){
1553 MRI_IMARR *qar ;
1554 if( nstk < 2 ) ERREX("no matrices") ;
1555 ima = IMARR_SUBIM(imstk,nstk-2) ;
1556 imb = IMARR_SUBIM(imstk,nstk-1) ;
1557 INIT_IMARR(qar) ; ADDTO_IMARR(qar,ima) ; ADDTO_IMARR(qar,imb) ;
1558 imc = mri_cat2D( 2 , 1 , 0 , NULL , qar ) ; FREE_IMARR(qar) ;
1559 if( imc == NULL ) ERREX("illegal Hglue") ;
1560 TRUNCATE_IMARR(imstk,nstk-2) ;
1561 ADDTO_IMARR(imstk,imc) ;
1562 }
1563
1564 /** matrix Vglue [04 Jan 2012] **/
1565
1566 else if( command_check(cmd,"Vglue") ){
1567 MRI_IMARR *qar ;
1568 if( nstk < 2 ) ERREX("no matrices") ;
1569 ima = IMARR_SUBIM(imstk,nstk-2) ;
1570 imb = IMARR_SUBIM(imstk,nstk-1) ;
1571 INIT_IMARR(qar) ; ADDTO_IMARR(qar,ima) ; ADDTO_IMARR(qar,imb) ;
1572 imc = mri_cat2D( 1 , 2 , 0 , NULL , qar ) ; FREE_IMARR(qar) ;
1573 if( imc == NULL ) ERREX("illegal Vglue") ;
1574 TRUNCATE_IMARR(imstk,nstk-2) ;
1575 ADDTO_IMARR(imstk,imc) ;
1576 }
1577
1578 /** recall named matrix to top of stack **/
1579
1580 else if( isalpha(cmd[0]) ){
1581 ii = matrix_name_lookup(cmd) ;
1582 if( ii >= 0 ){
1583 imc = mri_to_float(IMARR_SUBIM(matar,ii)) ;
1584 ADDTO_IMARR(imstk,imc) ;
1585 } else {
1586 sprintf(mess,"can't lookup name '%s'",cmd); ERREX(mess);
1587 }
1588 }
1589
1590 /** matrix duplicate on top of stack **/
1591
1592 else if( command_check(cmd,"dup") ){
1593 if( nstk < 1 ) ERREX("no matrix") ;
1594 imc = mri_to_float( IMARR_SUBIM(imstk,nstk-1) ) ;
1595 ADDTO_IMARR(imstk,imc) ;
1596 }
1597
1598 /** pop top of stack to oblivion **/
1599
1600 else if( command_check(cmd,"pop") ){
1601 if( nstk < 1 ) ERREX("no matrix") ;
1602 TRUNCATE_IMARR(imstk,nstk-1) ;
1603 }
1604
1605 /** swap top 2 matrices on stack **/
1606
1607 else if( command_check(cmd,"swap") ){
1608 if( nstk < 2 ) ERREX("no matrices") ;
1609 ima = IMARR_SUBIM(imstk,nstk-2) ;
1610 imb = IMARR_SUBIM(imstk,nstk-1) ;
1611 IMARR_SUBIM(imstk,nstk-2) = imb ;
1612 IMARR_SUBIM(imstk,nstk-1) = ima ;
1613 }
1614
1615 #if 0
1616 /** store a string (not really used anywhere) **/
1617
1618 else if( cmd[0] == '\'' || cmd[0] == '"' ){
1619 char *buf = strdup(cmd+1) ;
1620 imc = mri_new(1,1,MRI_float) ;
1621 for( ii=0 ; buf[ii] != cmd[0] && buf[ii] != '\0' ; ii++ ) ; /*nada*/
1622 buf[ii] = '\0' ; mri_add_name(buf,imc) ; free(buf) ;
1623 ADDTO_IMARR(imstk,imc) ;
1624 }
1625 #endif
1626
1627 /** stupid user **/
1628
1629 else {
1630 ERREX("unknown operation!?") ;
1631 }
1632 }
1633
1634 if(verb_rpn)fprintf(stderr,"++ exit mri_matrix_evalrpn\n") ;
1635
1636 NI_delete_str_array(sar) ;
1637 nstk = IMARR_COUNT(imstk) ;
1638 if( nstk > 0 ) ima = mri_to_float( IMARR_SUBIM(imstk,nstk-1) ) ;
1639 else ima = NULL ;
1640 DESTROY_IMARR(imstk) ;
1641 RETURN(ima) ;
1642 }
1643
1644 /*----------------------------------------------------------------------------*/
1645
mri_matrix_evalrpn_help(void)1646 char * mri_matrix_evalrpn_help(void)
1647 {
1648 static char *hs =
1649 "Evaluate a space delimited RPN matrix-valued expression:\n"
1650 "\n"
1651 " * The operations are on a stack, each element of which is a\n"
1652 " real-valued matrix.\n"
1653 " * N.B.: This is a computer-science stack of separate matrices.\n"
1654 " If you want to join two matrices in separate files\n"
1655 " into one 'stacked' matrix, then you must use program\n"
1656 " 1dcat to join them as columns, or the system program\n"
1657 " cat to join them as rows.\n"
1658 " * You can also save matrices by name in an internal buffer\n"
1659 " using the '=NAME' operation and then retrieve them later\n"
1660 " using just the same NAME.\n"
1661 " * You can read and write matrices from files stored in ASCII\n"
1662 " columns (.1D format) using the &read and &write operations.\n"
1663 " * The following 5 operations, input as a single string,\n"
1664 " \'&read(V.1D) &read(U.1D) &transp * &write(VUT.1D)\'\n"
1665 " - reads matrices V and U from disk (separately),\n"
1666 " - transposes U (on top of the stack) into U',\n"
1667 " - multiplies V and U' (the two matrices on top of the stack),\n"
1668 " - and writes matrix VU' out (the matrix left on the stack by '*').\n"
1669 " * Calculations are carried out in single precision ('float').\n"
1670 " * Operations mostly contain characters such as '&' and '*' that\n"
1671 " are special to Unix shells, so you'll probably need to put\n"
1672 " the arguments to this program in 'single quotes'.\n"
1673 " * You can use '%' or '@' in place of the '&' character, if you wish.\n"
1674 "\n"
1675 " STACK OPERATIONS\n"
1676 " -----------------\n"
1677 " number == push scalar value (1x1 matrix) on stack;\n"
1678 " a number starts with a digit or a minus sign\n"
1679 " =NAME == save a copy matrix on top of stack as 'NAME'\n"
1680 " NAME == push a copy of NAME-ed matrix onto top of stack;\n"
1681 " names start with an alphabetic character\n"
1682 " &clear == erase all named matrices (to save memory);\n"
1683 " does not affect the stack at all\n"
1684 " &purge == erase the stack;\n"
1685 " does not affect named matrices\n"
1686 " &read(FF) == read ASCII (.1D) file onto top of stack from file 'FF'\n"
1687 " &read4x4Xform(FF)\n"
1688 " == Similar to &read(FF), except that it expects data\n"
1689 " for a 12-parameter spatial affine transform.\n"
1690 " FF can contain 12x1, 1x12, 16x1, 1x16, 3x4, or\n"
1691 " 4x4 values. \n"
1692 " The read operation loads the data into a 4x4 matrix\n"
1693 " r11 r12 r13 r14\n"
1694 " r21 r22 r23 r24\n"
1695 " r31 r32 r33 r34\n"
1696 " 0.0 0.0 0.0 1.0\n"
1697 " This option was added to simplify the combination of \n"
1698 " linear spatial transformations. However, you are better \n"
1699 " off using cat_matvec for that purpose.\n"
1700 "\n"
1701 " &write(FF) == write top matrix to ASCII file to file 'FF';\n"
1702 " if 'FF' == '-', writes to stdout\n"
1703 " &transp == replace top matrix with its transpose\n"
1704 " &ident(N) == push square identity matrix of order N onto stack\n"
1705 " N is an fixed integer, OR\n"
1706 " &R to indicate the row dimension of the\n"
1707 " current top matrix, OR\n"
1708 " &C to indicate the column dimension of the\n"
1709 " current top matrix, OR\n"
1710 " =X to indicate the (1,1) element of the\n"
1711 " matrix named X\n"
1712 " &Psinv == replace top matrix with its pseudo-inverse\n"
1713 " [computed via SVD, not via inv(A'*A)*A']\n"
1714 " &Sqrt == replace top matrix with its square root\n"
1715 " [computed via Denman & Beavers iteration]\n"
1716 " N.B.: not all real matrices have real square\n"
1717 " roots, and &Sqrt will fail if you push it\n"
1718 " N.B.: the matrix must be square!\n"
1719 " &Pproj == replace top matrix with the projection onto\n"
1720 " its column space; Input=A; Output = A*Psinv(A)\n"
1721 " N.B.: result P is symmetric and P*P=P\n"
1722 " &Qproj == replace top matrix with the projection onto\n"
1723 " the orthogonal complement of its column space\n"
1724 " Input=A; Output=I-Pproj(A)\n"
1725 " * == replace top 2 matrices with their product;\n"
1726 " OR stack = [ ... C A B ] (where B = top) goes to\n"
1727 " &mult stack = [ ... C AB ]\n"
1728 " if either of the top matrices is a 1x1 scalar,\n"
1729 " then the result is the scalar multiplication of\n"
1730 " the other matrix; otherwise, matrices must conform\n"
1731 " + OR &add == replace top 2 matrices with sum A+B\n"
1732 " - OR &sub == replace top 2 matrices with difference A-B\n"
1733 " &dup == push duplicate of top matrix onto stack\n"
1734 " &pop == discard top matrix\n"
1735 " &swap == swap top two matrices (A <-> B)\n"
1736 " &Hglue == glue top two matrices together horizontally:\n"
1737 " stack = [ ... C A B ] goes to\n"
1738 " stack = [ ... C A|B ]\n"
1739 " this is like what program 1dcat does.\n"
1740 " &Vglue == glue top two matrices together vertically:\n"
1741 " stack = [ ... C A B ] goes to\n"
1742 "\n"
1743 " A\n"
1744 " stack = [ ... C - ]\n"
1745 " B\n"
1746 "\n"
1747 " this is like what program cat does.\n"
1748 ;
1749
1750 return hs ;
1751 }
1752