00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008
00009
00010
00011 int equal_bitvector_piece( MRI_IMAGE * b , MRI_IMAGE * c , int aa , int bb )
00012 {
00013 int ii ; byte * bv=MRI_BYTE_PTR(b) , * cv=MRI_BYTE_PTR(c) ;
00014
00015 if( aa < 0 ) aa = 0 ;
00016 if( bb >= b->nx ) bb = b->nx - 1 ;
00017 for( ii=aa ; ii <= bb ; ii++ ) if( bv[ii] != cv[ii] ) return 0 ;
00018 return 1;
00019 }
00020
00021 int equal_bitvector( MRI_IMAGE * b , MRI_IMAGE * c )
00022 {
00023 return equal_bitvector_piece( b , c , 0 , b->nx - 1 ) ;
00024 }
00025
00026
00027
00028 void randomize_bitvector_piece( MRI_IMAGE * b , int aa , int bb )
00029 {
00030 int ii ; byte * bv=MRI_BYTE_PTR(b) ;
00031
00032 if( aa < 0 ) aa = 0 ;
00033 if( bb >= b->nx ) bb = b->nx - 1 ;
00034 for( ii=aa ; ii <= bb ; ii++ ) bv[ii] = ( drand48() > 0.5 ) ;
00035 return ;
00036 }
00037
00038 void randomize_bitvector( MRI_IMAGE * b )
00039 {
00040 randomize_bitvector_piece( b , 0 , b->nx - 1 ) ;
00041 return ;
00042 }
00043
00044
00045
00046 void zero_bitvector_piece( MRI_IMAGE * b , int aa , int bb )
00047 {
00048 int ii ; byte * bv=MRI_BYTE_PTR(b) ;
00049
00050 if( aa < 0 ) aa = 0 ;
00051 if( bb >= b->nx ) bb = b->nx - 1 ;
00052 for( ii=aa ; ii <= bb ; ii++ ) bv[ii] = 0 ;
00053 return ;
00054 }
00055
00056 void zero_bitvector( MRI_IMAGE * b )
00057 {
00058 zero_bitvector_piece( b , 0 , b->nx - 1 ) ;
00059 return ;
00060 }
00061
00062
00063
00064 void one_bitvector_piece( MRI_IMAGE * b , int aa , int bb )
00065 {
00066 int ii ; byte * bv=MRI_BYTE_PTR(b) ;
00067
00068 if( aa < 0 ) aa = 0 ;
00069 if( bb >= b->nx ) bb = b->nx - 1 ;
00070 for( ii=aa ; ii <= bb ; ii++ ) bv[ii] = 1 ;
00071 return ;
00072 }
00073
00074 void one_bitvector( MRI_IMAGE * b )
00075 {
00076 one_bitvector_piece( b , 0 , b->nx - 1 ) ;
00077 return ;
00078 }
00079
00080
00081
00082 void invert_bitvector_piece( MRI_IMAGE * b , int aa , int bb )
00083 {
00084 int ii ; byte * bv=MRI_BYTE_PTR(b) ;
00085
00086 if( aa < 0 ) aa = 0 ;
00087 if( bb >= b->nx ) bb = b->nx - 1 ;
00088 for( ii=aa ; ii <= bb ; ii++ ) bv[ii] = 1 - bv[ii] ;
00089 return ;
00090 }
00091
00092 void invert_bitvector( MRI_IMAGE * b )
00093 {
00094 invert_bitvector_piece( b , 0 , b->nx - 1 ) ;
00095 return ;
00096 }
00097
00098
00099
00100 MRI_IMAGE * new_bitvector( int n )
00101 {
00102 MRI_IMAGE * b ;
00103 b = mri_new( n , 1 , MRI_byte ) ;
00104 return b ;
00105 }
00106
00107 MRI_IMAGE * copy_bitvector( MRI_IMAGE * b )
00108 {
00109 int ii ;
00110 MRI_IMAGE * c = mri_new( b->nx , 1 , MRI_byte ) ;
00111 memcpy( MRI_BYTE_PTR(c) , MRI_BYTE_PTR(b) , b->nx ) ;
00112 return c ;
00113 }
00114
00115
00116
00117 int count_bitvector( MRI_IMAGE * b )
00118 {
00119 int ii,ss , n=b->nx ;
00120 byte * bv = MRI_BYTE_PTR(b) ;
00121 for( ii=ss=0 ; ii < n ; ii++ ) if( bv[ii] ) ss++ ;
00122 return ss ;
00123 }
00124
00125
00126
00127 void normalize_floatvector( MRI_IMAGE * fim )
00128 {
00129 int ii , n=fim->nx ;
00130 float ff,gg , * far=MRI_FLOAT_PTR(fim) ;
00131
00132 for( ff=0.0,ii=0 ; ii < n ; ii++ ) ff += far[ii] ;
00133 ff /= n ;
00134 for( gg=0.0,ii=0 ; ii < n ; ii++ ) gg += SQR( (far[ii]-ff) ) ;
00135 if( gg <= 0.0 ) return ;
00136 gg = 1.0 / sqrt(gg) ;
00137 for( ii=0 ; ii < n ; ii++ ) far[ii] = (far[ii]-ff)*gg ;
00138 return ;
00139 }
00140
00141
00142
00143 float corr_floatbit( MRI_IMAGE * fim , MRI_IMAGE * bim )
00144 {
00145 int ii , n=fim->nx , ns ;
00146 float * far=MRI_FLOAT_PTR(fim) , ss ;
00147 byte * bar=MRI_BYTE_PTR(bim) ;
00148
00149 for( ss=0.0,ns=ii=0 ; ii < n ; ii++ )
00150 if( bar[ii] ){ ns++ ; ss += far[ii] ; }
00151
00152 if( ns == 0 || ns == n ) return 0.0 ;
00153
00154 ss *= sqrt( ((float) n) / (float)(ns*(n-ns)) ) ;
00155 return ss ;
00156 }
00157
00158
00159
00160 MRI_IMARR * init_bitvector_array( int nbv , MRI_IMAGE * fim )
00161 {
00162 MRI_IMARR * imar ;
00163 MRI_IMAGE * bim ;
00164 int ii , n=fim->nx ;
00165 byte * bar ; float * far=MRI_FLOAT_PTR(fim) ;
00166
00167 INIT_IMARR(imar) ;
00168 normalize_floatvector( fim ) ;
00169
00170 for( ii=0 ; ii < nbv ; ii++ ){
00171 bim = new_bitvector( n ) ;
00172 randomize_bitvector( bim ) ;
00173 bim->dx = corr_floatbit( fim , bim ) ;
00174 if( bim->dx < 0.0 ){ bim->dx = -bim->dx; invert_bitvector( bim ); }
00175 ADDTO_IMARR(imar,bim) ;
00176 }
00177
00178 return imar ;
00179 }
00180
00181
00182
00183 #define IRAN(j) (lrand48() % (j))
00184
00185 void evolve_bitvector_array( MRI_IMARR * bvar , MRI_IMAGE * fim )
00186 {
00187 static int nrvec=-1 ;
00188 static float * rvec=NULL ;
00189
00190 int ii,nn=fim->nx , nbv=IMARR_COUNT(bvar) , aa,bb , qbv ;
00191 MRI_IMAGE * bim , * qim ;
00192
00193
00194
00195 for( ii=0 ; ii < nbv ; ii++ ){
00196
00197 bim = IMARR_SUBIMAGE(bvar,ii) ;
00198
00199 aa = IRAN(nn) ; bb = aa + IRAN(5) ; if( bb >= nn ) bb = nn-1 ;
00200
00201 qim = copy_bitvector(bim) ;
00202 zero_bitvector_piece( qim , aa , bb ) ;
00203 if( equal_bitvector_piece(bim,qim,aa,bb) ){
00204 mri_free(qim) ;
00205 } else {
00206 qim->dx = corr_floatbit( fim , qim ) ;
00207 if( qim->dx < 0.0 ){ qim->dx = -qim->dx; invert_bitvector( qim ); }
00208 ADDTO_IMARR(bvar,qim) ;
00209 }
00210
00211 qim = copy_bitvector(bim) ;
00212 one_bitvector_piece( qim , aa , bb ) ;
00213 if( equal_bitvector_piece(bim,qim,aa,bb) ){
00214 mri_free(qim) ;
00215 } else {
00216 qim->dx = corr_floatbit( fim , qim ) ;
00217 if( qim->dx < 0.0 ){ qim->dx = -qim->dx; invert_bitvector( qim ); }
00218 ADDTO_IMARR(bvar,qim) ;
00219 }
00220
00221 qim = copy_bitvector(bim) ;
00222 randomize_bitvector_piece( qim , aa , bb ) ;
00223 if( equal_bitvector_piece(bim,qim,aa,bb) ){
00224 mri_free(qim) ;
00225 } else {
00226 qim->dx = corr_floatbit( fim , qim ) ;
00227 if( qim->dx < 0.0 ){ qim->dx = -qim->dx; invert_bitvector( qim ); }
00228 ADDTO_IMARR(bvar,qim) ;
00229 }
00230
00231 qim = copy_bitvector(bim) ;
00232 invert_bitvector_piece( qim , aa , bb ) ;
00233 if( equal_bitvector_piece(bim,qim,aa,bb) ){
00234 mri_free(qim) ;
00235 } else {
00236 qim->dx = corr_floatbit( fim , qim ) ;
00237 if( qim->dx < 0.0 ){ qim->dx = -qim->dx; invert_bitvector( qim ); }
00238 ADDTO_IMARR(bvar,qim) ;
00239 }
00240 }
00241
00242
00243
00244 qbv = IMARR_COUNT(bvar) ;
00245 if( nrvec < qbv ){
00246 if( rvec != NULL ) free(rvec) ;
00247 rvec = (float *) malloc(sizeof(float)*qbv) ;
00248 nrvec = qbv ;
00249 }
00250 for( ii=0 ; ii < qbv ; ii++ ) rvec[ii] = -fabs(IMARR_SUBIM(bvar,ii)->dx) ;
00251
00252 qsort_floatstuff( qbv , rvec , (void **) bvar->imarr ) ;
00253
00254 TRUNCATE_IMARR( bvar , nbv ) ;
00255 return ;
00256 }
00257
00258
00259
00260 int main( int argc , char * argv[] )
00261 {
00262 MRI_IMAGE * fim , * bim ;
00263 MRI_IMARR * bvar ;
00264 int ii , nv , nite=0 , neq=0 ;
00265 float fold , fnew ;
00266
00267 if( argc < 2 ){printf("Usage: bitvec fname.1D > bname.1D\n");exit(0);}
00268
00269 fim = mri_read_1D( argv[1] ) ;
00270 if( fim == NULL ){fprintf(stderr,"Can't read %s\n",argv[1]);exit(1);}
00271 fprintf(stderr,"vector length = %d\n",fim->nx) ;
00272
00273 srand48((long)time(NULL)) ;
00274
00275 bvar = init_bitvector_array( 256 , fim ) ;
00276 fold = fabs(IMARR_SUBIM(bvar,0)->dx) ;
00277 fprintf(stderr,"fold = %7.4f\n",fold) ;
00278
00279 while(1){
00280 evolve_bitvector_array( bvar , fim ) ;
00281 nv = IMARR_COUNT(bvar) ;
00282 nite++ ;
00283 #if 0
00284 fprintf(stderr,"nite=%d\n",nite) ;
00285 for( ii=0 ; ii < nv ; ii++ )
00286 fprintf(stderr," %7.4f",IMARR_SUBIM(bvar,ii)->dx) ;
00287 fprintf(stderr,"\n\n") ;
00288 #endif
00289
00290 fnew = fabs(IMARR_SUBIM(bvar,0)->dx) ;
00291 if( fnew <= fold ){
00292 neq++ ; if( neq == 10 ) break ;
00293 } else {
00294 neq = 0 ;
00295 fold = fnew ;
00296 fprintf(stderr,"%d: %7.4f\n",nite,fold) ;
00297 }
00298 }
00299
00300 bim = IMARR_SUBIM(bvar,0) ;
00301 for( ii=0 ; ii < fim->nx ; ii++ )
00302 printf(" %d\n",(int)MRI_BYTE_PTR(bim)[ii]) ;
00303
00304 exit(0) ;
00305 }
