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 #include <stdio.h>
8 #include <math.h>
9 #include <stdlib.h>
10
11 #include "pcor.h"
12
13 /************************************************************************/
14 /************************************************************************/
15
16 /***
17 recursive calculation of partial correlation coefficients for
18 a lot of voxels at once. -- RWCox, Feb 1994
19 ***/
20
21 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
22
23 /***
24 create a new references data structure:
25 input: numref = number of reference vectors to allow for
26 output: pointer to the data structure
27 ***/
28
new_references(numref)29 references * new_references(numref)
30 int numref;
31 {
32 references *ref ;
33 int ii,jj ;
34
35 /*** check input for reasonableness ***/
36
37 if( numref < 1 ){
38 fprintf( stderr , "new_references called with numref=%d\n" , numref ) ;
39 exit(-1) ;
40 }
41
42 /*** allocate storage for top level data ***/
43
44 /* malloc->calloc 13 Feb 2009 [lesstif patrol] */
45 ref = (references *) calloc( 1, sizeof(references) ) ;
46 if( ref == NULL ){
47 fprintf( stderr , "new_references: malloc error for base\n" ) ;
48 exit(-1) ;
49 }
50 ref->nref = numref ;
51 ref->nupdate = 0 ; /* June 1995: not updated at all yet */
52
53 /*** allocate storage for Cholesky factor
54 (an array of rows, row #ii is length ii+1, for ii=0..numref-1) ***/
55
56 ref->chol = (ref_float **) malloc( sizeof(ref_float *) * numref ) ;
57 if( ref->chol == NULL ){
58 fprintf( stderr , "new_references: malloc error for chol\n" ) ;
59 exit(-1) ;
60 }
61
62 for( ii=0 ; ii < numref ; ii++ ){
63 ref->chol[ii] = (ref_float *) malloc( sizeof(ref_float) * (ii+1) ) ;
64 if( ref->chol[ii] == NULL ){
65 fprintf( stderr , "new_references: malloc error for chol[ii]\n" ) ;
66 exit(-1) ;
67 }
68 }
69
70 /*** allocate storage for vectors of alpha, f, g ***/
71
72 ref->alp = (ref_float *) malloc( sizeof(ref_float) * numref ) ;
73 ref->ff = (ref_float *) malloc( sizeof(ref_float) * numref ) ;
74 ref->gg = (ref_float *) malloc( sizeof(ref_float) * numref ) ;
75
76 if( ref->alp == NULL || ref->ff == NULL || ref->gg == NULL ){
77 fprintf( stderr , "new_references: malloc error for data\n" ) ;
78 exit(-1) ;
79 }
80
81 /*** initialize Cholesky factor ***/
82
83 for( ii=0 ; ii < numref ; ii++ ){
84 for( jj=0 ; jj < ii ; jj++ ) RCH(ref,ii,jj) = 0.0 ;
85 RCH(ref,ii,ii) = REF_EPS ;
86 #ifdef OV_DEBUG2
87 ref->alp[ii] = ref->ff[ii] = ref->gg[ii] = 0.0 ;
88 #endif
89 }
90
91 #ifdef OV_DEBUG2
92 ref->betasq = 0.0 ;
93 #endif
94
95 return ref ;
96 }
97
98 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
99
100 /***
101 update a references structure:
102 input: vec = pointer to nref values of the reference vectors
103 at the new time point
104 ref = references data structure
105 output: ref is updated;
106 the Cholesky factor is modified via the Carlson algorithm,
107 the alpha, f, and g factors are saved for later use
108 ***/
109
update_references(vec,ref)110 void update_references(vec, ref)
111 float *vec;
112 references *ref;
113 {
114 int nr = ref->nref , jj,kk ;
115 ref_float bold , bnew , aaa,fff,ggg ;
116 static ref_float * zz = NULL ;
117 static int zz_size = -1 ;
118
119 #ifdef OV_DEBUG2
120 static ref_float qinput[50] ; /* to hold the sums of squares of inputs */
121 #endif
122
123 /*** copy vector data into local storage ***/
124
125 if( zz_size < nr ){ /* get new space, if not enough is present */
126
127 if( zz != NULL ) free( zz ) ;
128 zz = (ref_float *) malloc( sizeof(ref_float) * nr ) ;
129 zz_size = nr ;
130 if( zz == NULL ){
131 fprintf( stderr , "update_references: cannot malloc!\n" ) ;
132 exit(-1) ;
133 }
134 }
135
136 for( jj=0 ; jj < nr ; jj++) zz[jj] = (ref_float) vec[jj] ;
137
138 #ifdef OV_DEBUG2
139 for( jj=0 ; jj < nr ; jj++) qinput[jj] += SQR(zz[jj]) ; /* for later */
140
141 REF_DUMP(ref,"before update") ;
142 fprintf(stderr," input vec= ") ;
143 for( jj=0 ; jj < nr ; jj++ ) fprintf(stderr,"%11.4e ",zz[jj]) ;
144 fprintf(stderr,"\n") ;
145 #endif
146
147 /*** Carlson algorithm ***/
148
149 bold = 1.0 ;
150
151 for( jj=0 ; jj < nr ; jj++ ){
152
153 aaa = zz[jj] / RCH(ref,jj,jj) ; /* alpha */
154 bnew = sqrt( bold*bold + aaa*aaa ) ; /* new beta */
155 fff = bnew / bold ; /* f factor */
156 ggg = aaa / (bnew*bold) ; /* g factor */
157 bold = bnew ; /* new beta becomes old beta */
158
159 ref->alp[jj] = aaa ; /* save these for later use */
160 ref->ff[jj] = fff ;
161 ref->gg[jj] = ggg ;
162
163 for( kk=jj ; kk < nr ; kk++ ){
164 zz[kk] -= aaa * RCH(ref,kk,jj) ;
165 RCH(ref,kk,jj) = fff * RCH(ref,kk,jj) + ggg * zz[kk] ;
166 }
167 }
168
169 ref->betasq = 1.0 / ( bold * bold ) ; /* and save this too! */
170
171 #ifdef OV_DEBUG2
172 REF_DUMP(ref,"after update") ;
173 fprintf(stderr," qsum of input vecs= ") ;
174 for( jj=0 ; jj < nr ; jj++ ) fprintf(stderr,"%11.4e ",qinput[jj]) ;
175 fprintf(stderr,"\n") ;
176 #endif
177
178 (ref->nupdate)++ ; /* June 1995: another update! */
179 return ;
180 }
181
182 /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
183
184 /***
185 create a new voxel partial correlation data structure
186 inputs: numvox = number of voxels in image
187 numref = number of reference vectors to allow for
188 output: pointer to voxel_corr data structure
189 ***/
190
new_voxel_corr(numvox,numref)191 voxel_corr * new_voxel_corr(numvox, numref)
192 int numvox;
193 int numref;
194 {
195 int vox , jj ;
196 voxel_corr *vc ;
197
198 /*** check input for OK-osity ***/
199
200 if( numvox < 1 ){
201 fprintf( stderr , "new_voxel_corr: numvox=%d\n" , numvox ) ;
202 exit(-1) ;
203 }
204
205 /*** get the base storage ***/
206
207 vc = (voxel_corr *) malloc( sizeof(voxel_corr) ) ;
208 if( vc == NULL ){
209 fprintf( stderr , "new_voxel_corr: cannot malloc base\n" ) ;
210 exit(-1) ;
211 }
212
213 /*** setup the references common to all voxels ***/
214
215 vc->nvox = numvox ;
216 vc->nref = numref ;
217 vc->nupdate = 0 ; /* June 1995: not updated at all yet */
218
219 /*** setup the storage of the last row for each voxel ***/
220
221 vc->chrow = (ref_float *)malloc( sizeof(ref_float) * numvox*(numref+1) );
222 if( vc->chrow == NULL ){
223 fprintf( stderr , "new_voxel_corr: cannot malloc last rows\n" ) ;
224 exit(-1) ;
225 }
226
227 /*** initialize each voxel ***/
228
229 for( vox=0 ; vox < numvox ; vox++ ){
230 for( jj=0 ; jj < numref ; jj++ ) VCH(vc,vox,jj) = 0 ;
231 VCH(vc,vox,numref) = REF_EPS ;
232 }
233
234 return vc ;
235 }
236
237 /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
238
239 /*** de-allocate a references data structure ***/
240
free_references(ref)241 void free_references(ref)
242 references *ref;
243 {
244 int ii , nr ;
245
246 if( ref == NULL ) return ;
247
248 nr = ref->nref ; if( nr <= 0 ) return ;
249
250 free(ref->alp) ; free(ref->ff) ; free(ref->gg) ;
251
252 for( ii=0 ; ii < nr ; ii++ ) free( ref->chol[ii] ) ;
253
254 free(ref->chol) ; free(ref) ;
255
256 return ;
257 }
258
259 /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
260
free_voxel_corr(vc)261 void free_voxel_corr(vc)
262 voxel_corr *vc;
263 {
264 if( vc != NULL ){
265 free( vc->chrow ) ;
266 free( vc ) ;
267 }
268 return ;
269 }
270
271 /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
272
273 /*** update all voxels with a new array of data
274 inputs: vdata = array[nvox] of new data for each voxel
275 ref = pointer to references structure to use
276 vc = pointer to correlation data structure
277 output: updated vc
278 ***/
279
update_voxel_corr(vdata,ref,vc)280 void update_voxel_corr(vdata, ref, vc)
281 vox_data *vdata;
282 references *ref;
283 voxel_corr *vc;
284 {
285 int vox , jj , /* loop indices */
286 nv = vc->nvox , /* number of voxels */
287 nr = vc->nref ; /* number of references */
288
289 ref_float *aaa = ref->alp ,
290 *fff = ref->ff ,
291 *ggg = ref->gg ;
292
293 ref_float zz ,
294 bq = ref->betasq ;
295
296 #ifdef OV_DEBUG2
297 static ref_float qvox = 0.0 ;
298 #endif
299
300 /*** check inputs for OK-ness ***/
301
302 if( vc->nref != ref->nref ){
303 fprintf( stderr , "update_voxel_corr: reference size mismatch!\n" ) ;
304 exit(-1) ;
305 }
306
307 #ifdef OV_DEBUG2
308 VOX_DUMP(vc,VD,"before update") ;
309 #ifdef VOX_SHORT
310 fprintf(stderr," integer input data = %d\n" , vdata[VD] ) ;
311 #else
312 fprintf(stderr," float input data = %11.4e\n" , vdata[VD] ) ;
313 #endif
314 qvox += SQR(vdata[VD]) ;
315 #endif
316
317 /** innermost loop expansion is for speedup if nref is small, if enabled **/
318
319 #ifdef EXPAND_UPDATE
320 #define UPZZ(j) zz -= aaa[j] * VCH(vc,vox,j)
321 #define UPCH(j) VCH(vc,vox,j) = fff[j] * VCH(vc,vox,j) + ggg[j] * zz
322 #define UPLL(j) VCH(vc,vox,j) += bq * zz * zz
323
324 switch( nr ){
325 default:
326 #endif
327
328 /*** for each voxel ***/
329
330 for( vox=0 ; vox < nv ; vox++ ){
331
332 /*** update last row of each Cholesky factor ***/
333
334 zz = (ref_float) vdata[vox] ;
335 for( jj=0 ; jj < nr ; jj++ ){
336 zz -= aaa[jj] * VCH(vc,vox,jj) ;
337 VCH(vc,vox,jj) = fff[jj] * VCH(vc,vox,jj) + ggg[jj] * zz ;
338 }
339 VCH(vc,vox,nr) += bq * zz * zz ; /* square of true Cholesky diagonal */
340 }
341
342 #ifdef EXPAND_UPDATE
343 break ;
344
345 case 1:
346 for( vox=0 ; vox < nv ; vox++ ){
347 zz = (ref_float) vdata[vox] ;
348 UPZZ(0) ; UPCH(0) ; UPLL(1) ;
349 }
350 break ;
351
352 case 2:
353 for( vox=0 ; vox < nv ; vox++ ){
354 zz = (ref_float) vdata[vox] ;
355 UPZZ(0) ; UPCH(0) ;
356 UPZZ(1) ; UPCH(1) ; UPLL(2) ;
357 }
358 break ;
359
360 case 3:
361 for( vox=0 ; vox < nv ; vox++ ){
362 zz = (ref_float) vdata[vox] ;
363 UPZZ(0) ; UPCH(0) ;
364 UPZZ(1) ; UPCH(1) ;
365 UPZZ(2) ; UPCH(2) ; UPLL(3) ;
366 }
367 break ;
368
369 case 4:
370 for( vox=0 ; vox < nv ; vox++ ){
371 zz = (ref_float) vdata[vox] ;
372 UPZZ(0) ; UPCH(0) ;
373 UPZZ(1) ; UPCH(1) ;
374 UPZZ(2) ; UPCH(2) ;
375 UPZZ(3) ; UPCH(3) ; UPLL(4) ;
376 }
377 break ;
378
379 case 5:
380 for( vox=0 ; vox < nv ; vox++ ){
381 zz = (ref_float) vdata[vox] ;
382 UPZZ(0) ; UPCH(0) ;
383 UPZZ(1) ; UPCH(1) ;
384 UPZZ(2) ; UPCH(2) ;
385 UPZZ(3) ; UPCH(3) ;
386 UPZZ(4) ; UPCH(4) ; UPLL(5) ;
387 }
388 break ;
389
390 case 6:
391 for( vox=0 ; vox < nv ; vox++ ){
392 zz = (ref_float) vdata[vox] ;
393 UPZZ(0) ; UPCH(0) ;
394 UPZZ(1) ; UPCH(1) ;
395 UPZZ(2) ; UPCH(2) ;
396 UPZZ(3) ; UPCH(3) ;
397 UPZZ(4) ; UPCH(4) ;
398 UPZZ(5) ; UPCH(5) ; UPLL(6) ;
399 }
400 break ;
401
402 case 7:
403 for( vox=0 ; vox < nv ; vox++ ){
404 zz = (ref_float) vdata[vox] ;
405 UPZZ(0) ; UPCH(0) ;
406 UPZZ(1) ; UPCH(1) ;
407 UPZZ(2) ; UPCH(2) ;
408 UPZZ(3) ; UPCH(3) ;
409 UPZZ(4) ; UPCH(4) ;
410 UPZZ(5) ; UPCH(5) ;
411 UPZZ(6) ; UPCH(6) ; UPLL(7) ;
412 }
413 break ;
414
415 case 8:
416 for( vox=0 ; vox < nv ; vox++ ){
417 zz = (ref_float) vdata[vox] ;
418 UPZZ(0) ; UPCH(0) ;
419 UPZZ(1) ; UPCH(1) ;
420 UPZZ(2) ; UPCH(2) ;
421 UPZZ(3) ; UPCH(3) ;
422 UPZZ(4) ; UPCH(4) ;
423 UPZZ(5) ; UPCH(5) ;
424 UPZZ(6) ; UPCH(6) ;
425 UPZZ(7) ; UPCH(7) ; UPLL(8) ;
426 }
427 break ;
428
429 case 9:
430 for( vox=0 ; vox < nv ; vox++ ){
431 zz = (ref_float) vdata[vox] ;
432 UPZZ(0) ; UPCH(0) ;
433 UPZZ(1) ; UPCH(1) ;
434 UPZZ(2) ; UPCH(2) ;
435 UPZZ(3) ; UPCH(3) ;
436 UPZZ(4) ; UPCH(4) ;
437 UPZZ(5) ; UPCH(5) ;
438 UPZZ(6) ; UPCH(6) ;
439 UPZZ(7) ; UPCH(7) ;
440 UPZZ(8) ; UPCH(8) ; UPLL(9) ;
441 }
442 break ;
443
444 }
445 #endif /* EXPAND_UPDATE */
446
447 #ifdef OV_DEBUG2
448 VOX_DUMP(vc,VD,"after update") ;
449 fprintf(stderr," qsum of vox[VD]=%11.4e\n",qvox) ;
450 #endif
451
452 (vc->nupdate)++ ; /* June 1995: another update */
453 return ;
454 }
455
456 /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
457
458 /*** compute the partial correlation coefficients ***/
459
get_pcor(ref,vc,pcor)460 void get_pcor(ref, vc, pcor)
461 references *ref;
462 voxel_corr *vc;
463 float *pcor;
464 {
465 int vox , nv = vc->nvox , nr = vc->nref ;
466 ref_float den ;
467 #define DENEPS 1.e-5
468
469 /*** check inputs for OK-ness ***/
470
471 if( vc->nref != ref->nref ){
472 fprintf( stderr , "get_pcor: reference size mismatch!\n" ) ;
473 exit(-1) ;
474 }
475
476 /*** Work ***/
477
478 for( vox=0 ; vox < nv ; vox++ ){
479
480 den = VCH(vc,vox,nr) ;
481 if( den > DENEPS ){
482 pcor[vox] = VCH(vc,vox,nr-1)
483 / sqrt( den + SQR(VCH(vc,vox,nr-1)) ) ;
484 } else {
485 pcor[vox] = 0.0 ;
486 }
487
488 }
489
490 #ifdef OV_DEBUG2
491 fprintf(stderr,"get_pcor: pcor[VD]=%11.4e\n",pcor[VD]) ;
492 #endif
493
494 return ;
495 }
496
497 /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
498
499 /*** get activation coefficient ***/
500
get_coef(ref,vc,coef)501 void get_coef(ref, vc, coef)
502 references *ref;
503 voxel_corr *vc;
504 float *coef;
505 {
506 int vox , nv = vc->nvox , nr = vc->nref ;
507 ref_float scale ;
508
509 /*** check inputs for OK-ness ***/
510
511 if( vc->nref != ref->nref ){
512 fprintf( stderr , "get_coef: reference size mismatch!\n" ) ;
513 exit(-1) ;
514 }
515
516 /*** Work ***/
517
518 scale = 1.0 / RCH(ref,nr-1,nr-1) ;
519
520 for( vox=0 ; vox < nv ; vox++ ){
521 coef[vox] = scale * VCH(vc,vox,nr-1) ;
522 }
523
524 #ifdef OV_DEBUG2
525 fprintf(stderr,"get_coef: coef[VD]=%11.4e\n",coef[VD]) ;
526 #endif
527
528 return ;
529 }
530
531 /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
532
533 /*** get variance estimate (June 1995) ***/
534
get_variance(vc,var)535 void get_variance(vc, var)
536 voxel_corr *vc;
537 float *var;
538 {
539 int vox , nv = vc->nvox , nr = vc->nref , nup = vc->nupdate ;
540 ref_float scale ;
541
542 /*** check inputs for OK-ness ***/
543
544 if( nup <= nr ){
545 fprintf(stderr,"get_variance: not enough data to compute!\n") ;
546 for( vox=0 ; vox < nv ; vox++ ) var[vox] = 0.0 ;
547 return ;
548 }
549
550 /*** Work ***/
551
552 scale = 1.0 / ( nup - nr ) ;
553
554 for( vox=0 ; vox < nv ; vox++ ){
555 var[vox] = scale * VCH(vc,vox,nr) ;
556 }
557
558 #ifdef OV_DEBUG2
559 fprintf(stderr,"get_variance: var[VD]=%11.4e\n",var[VD]) ;
560 #endif
561
562 return ;
563 }
564
565 /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
566
567 /*** Get least squares fit coefficients (all of them, Frank).
568 "fit" is an array of pointers to floats, of length nref.
569 If fit[j] != NULL, then it points to an array of size nvox that
570 will get the coefficient for reference #j (j=0..nref-1). [June 1995] ***/
571
get_lsqfit(ref,vc,fit)572 void get_lsqfit(ref, vc, fit)
573 references *ref;
574 voxel_corr *vc;
575 float *fit[] ;
576 {
577 int vox,jj,kk , nv = vc->nvox , nr = vc->nref ;
578 ref_float sum ;
579 ref_float * ff ;
580
581 /*** check inputs for OK-ness ***/
582
583 if( vc->nref != ref->nref ){
584 fprintf( stderr , "get_lsqfit: reference size mismatch!\n" ) ;
585 exit(-1) ;
586 }
587
588 kk = 0 ;
589 for( jj=0 ; jj < nr ; jj++ ) kk += (fit[jj] != NULL) ;
590 if( kk == 0 ){
591 fprintf(stderr,"get_lsqfit: NO OUTPUT REQUESTED!\n") ;
592 return ;
593 }
594
595 ff = (ref_float *) malloc( sizeof(ref_float) * nr ) ;
596 if( ff == NULL ){
597 fprintf( stderr, "get_lsqfit: cannot malloc workspace!\n") ;
598 exit(-1) ;
599 }
600
601 /*** for each voxel, compute the nr fit coefficients (backwards) ***/
602
603 for( vox=0 ; vox < nv ; vox++ ){
604
605 for( jj=nr-1 ; jj >=0 ; jj-- ){
606 sum = VCH(vc,vox,jj) ;
607 for( kk=jj+1 ; kk < nr ; kk++ ) sum -= ff[kk] * RCH(ref,kk,jj) ;
608 ff[jj] = sum / RCH(ref,jj,jj) ;
609 }
610
611 for( jj=0 ; jj < nr ; jj++ )
612 if( fit[jj] != NULL ) fit[jj][vox] = ff[jj] ;
613 }
614
615 free( ff ) ;
616 return ;
617 }
618
619
620 /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
621
622 /*** get correlation and thresholded alpha:
623
624 only |pcor| >= pcthresh will be computed;
625 only those voxels will have coef computed;
626 if cothresh > 0, then voxels whose |coef| is less than
627 cothresh * max|coef| will be also be set to zero
628 ***/
629
get_pcor_thresh_coef(ref,vc,pcthresh,cothresh,pcor,coef,thr)630 void get_pcor_thresh_coef(ref, vc, pcthresh, cothresh, pcor, coef, thr)
631 references *ref;
632 voxel_corr *vc;
633 float pcthresh;
634 float cothresh;
635 float *pcor;
636 float *coef;
637 thresh_result *thr;
638 {
639 int vox , nv = vc->nvox , nr = vc->nref ;
640 ref_float den , num , scale ;
641 #define DENEPS 1.e-5
642
643 float pc , co , thfac ;
644 int npc_pos=0 , npc_neg=0 , nco_pos=0 , nco_neg=0 ;
645 float mpc_pos=0., mpc_neg=0., mco_pos=0., mco_neg=0.;
646
647 int do_pcth , do_coth ;
648
649 /*** check inputs for OK-ness ***/
650
651 if( vc->nref != ref->nref ){
652 fprintf( stderr , "get_pcor: reference size mismatch!\n" ) ;
653 exit(-1) ;
654 }
655
656 scale = 1.0 / RCH(ref,nr-1,nr-1) ; /* for coef calculation */
657 thfac = SQR(pcthresh)/(1.0-SQR(pcthresh)) ;
658 do_pcth = pcthresh <= 0.0 ; /* whether to do these tests */
659 do_coth = cothresh > 0.0 ;
660
661 /*** Compute pcor and coef, thresholded on pcthresh ***/
662
663 for( vox=0 ; vox < nv ; vox++ ){
664
665 den = VCH(vc,vox,nr) ;
666 num = VCH(vc,vox,nr-1) ;
667 if( do_pcth || SQR(num) > thfac*den ){ /* fancy threshold test */
668
669 pc = pcor[vox] = num / sqrt(den+SQR(num)) ;
670 co = coef[vox] = scale * num ;
671
672 if( pc > 0 ){
673 npc_pos++ ;
674 if( mpc_pos < pc ) mpc_pos = pc ;
675 if( mco_pos < co ) mco_pos = co ;
676 } else {
677 npc_neg++ ;
678 if( mpc_neg > pc ) mpc_neg = pc ;
679 if( mco_neg > co ) mco_neg = co ;
680 }
681
682 } else { /* fails pcor thresh */
683 pcor[vox] = coef[vox] = 0.0 ;
684 }
685
686 }
687
688 nco_pos = npc_pos ;
689 nco_neg = npc_neg ;
690
691 /*** threshold coef on cothresh as well ***/
692
693 if( do_coth && nco_pos+nco_neg > 0 ){
694
695 thfac = cothresh * MAX(mco_pos,-mco_neg) ;
696
697 for( vox=0 ; vox < nv ; vox++ ){
698 if( coef[vox] > 0.0 && coef[vox] < thfac ){
699 coef[vox] = 0.0 ;
700 nco_pos-- ;
701 } else if( coef[vox] < 0.0 && coef[vox] > -thfac ){
702 coef[vox] = 0.0 ;
703 nco_neg-- ;
704 }
705 }
706 }
707
708 /*** load threshold output report ***/
709
710 thr->num_pcor_pos = npc_pos ;
711 thr->num_pcor_neg = npc_neg ;
712 thr->max_pcor_pos = mpc_pos ;
713 thr->max_pcor_neg = mpc_neg ;
714
715 thr->num_coef_pos = nco_pos ;
716 thr->num_coef_neg = nco_neg ;
717 thr->max_coef_pos = mco_pos ;
718 thr->max_coef_neg = mco_neg ;
719
720 return ;
721 }
722