1 #include "mrilib.h"
2
3 /*----------------------------------------------------------------------*/
4 /*! Return a measure of the difference between 2 images bim and nim,
5 possibly with a mask to indicate which voxels to include.
6
7 The difference is defined as dd below:
8
9 ss = min sum [ bim[i] - a * nim[i] ]**2
10 a i
11
12 dd = sqrt( ss / #voxels ) i.e., the RMS difference.
13
14 That is, we least squares fit image nim to bim by a scale factor,
15 and the residual determins the difference. The minimizing value of
16 'a' is given by
17
18 a = BN / NN where BN = sum [ bim[i]*nim[i] ]
19 NN = sum [ nim[i]**2 ]
20
21 so the difference is given by
22
23 diff = BB - (BN)**2/NN where BB = sum[ bim[i]**2 ]
24
25 - A negative return value indicates an error.
26 - If bim=0, then the return value will be 0.
27 - If bim!=0 and nim=0, then the return value is computed with a=0
28 (that is, it will be the RMS value of bim)
29 ------------------------------------------------------------------------*/
30
mri_scaled_diff(MRI_IMAGE * bim,MRI_IMAGE * nim,MRI_IMAGE * msk)31 float mri_scaled_diff( MRI_IMAGE *bim , MRI_IMAGE *nim , MRI_IMAGE *msk )
32 {
33 int nvox , ii , nmask=0 ;
34 MRI_IMAGE *fim , *gim ;
35 float *far , *gar , sdif , ff,gg,fg ;
36 byte *mar=NULL ;
37
38 ENTRY("mri_scaled_diff") ;
39
40 if( bim == NULL || nim == NULL ) RETURN(-1.0f) ;
41
42 nvox = bim->nvox ; if( nim->nvox != nvox ) RETURN(-1.0f) ;
43
44 if( msk != NULL && msk->kind == MRI_byte && msk->nvox == nvox ){
45 mar = MRI_BYTE_PTR(msk) ;
46 nmask = THD_countmask( nvox , mar ) ;
47 if( nmask < 3 ){ nmask = 0 ; mar = NULL ; }
48 }
49
50 fim = (bim->kind != MRI_float) ? mri_to_float(bim) : bim ;
51 gim = (nim->kind != MRI_float) ? mri_to_float(nim) : nim ;
52
53 far = MRI_FLOAT_PTR(fim) ; gar = MRI_FLOAT_PTR(gim) ;
54
55 ff = gg = fg = 0.0f ;
56 for( ii=0 ; ii < nvox ; ii++ ){
57 if( nmask == 0 || mar[ii] != 0 ){
58 ff += far[ii] * far[ii] ;
59 gg += gar[ii] * gar[ii] ;
60 fg += far[ii] * gar[ii] ;
61 }
62 }
63 if( gg > 0.0f ){ /* normal case */
64 sdif = ff - fg*fg/gg ;
65 if( sdif > 0.0f )
66 sdif = sqrt( sdif / ((nmask > 0) ? nmask : nvox) ) ;
67
68 } else if( ff <= 0.0f ) /* far=gar=0 */
69 sdif = 0.0f ;
70 else /* gar=0 far!=0 */
71 sdif = sqrt( ff / ((nmask > 0) ? nmask : nvox) ) ;
72
73 /** done **/
74
75 if( fim != bim ) mri_free(fim) ;
76 if( gim != nim ) mri_free(gim) ;
77 RETURN(sdif) ;
78 }
79