00001 #include "mrilib.h"
00002
00003 void print_results( char * , int , float *, float * ) ;
00004
00005 static float threshold=0.0 , cer=0.0 , cjv=0.0 , gmcount=0.0 , wmcount=0.0 ;
00006
00007
00008 static double apar=0.0 ;
00009
00010 static double pmodel_pdf( double x )
00011 {
00012 double q = 1.0-x*x ;
00013 return (q <= 0.0) ? 0.0
00014 : 0.9375*q*q
00015 * (1.0 + apar*x*fabs(x)) ;
00016 }
00017
00018 static double pmodel_cdf( double x )
00019 {
00020 double q , val , sss ;
00021 if( x <= -1.0 ) return 0.0 ;
00022 if( x >= 1.0 ) return 1.0 ;
00023 q = x*x ;
00024 val = 0.5 + x*(15.0-10.0*q+3.0*q*q)/16.0 - apar/14.0 ;
00025 sss = 0.9375 * (apar*x*q)*(1.0/3.0-0.4*q+q*q/7.0) ;
00026 if( x < 0.0 ) val -= sss ;
00027 else val += sss ;
00028 return val ;
00029 }
00030
00031 #define pmodel_bin(x1,x2,pk,wd) (pmodel_cdf((x2-pk)/wd)-pmodel_cdf((x1-pk)/wd))
00032
00033 #define pmodel_mean(pk,wd) ((pk)+0.078125*apar*(wd))
00034
00035 #define pmodel_sigma(wd) (0.002232142857*(wd)*sqrt(28672.-1225.*apar*apar))
00036
00037
00038
00039 int main( int argc , char * argv[] )
00040 {
00041 THD_3dim_dataset *dset , *mset ;
00042 byte *mask ;
00043 int iarg=1 , nmask,nfill , dilate=1 , dd,nx,ny,nz,nxy,nxyz , ii,jj,kk ;
00044 float SIhh=130.0 ;
00045 int SIax=0 , SIbot,SItop ;
00046 short *sar , *mar ;
00047 float pval[128] , wval[128] ;
00048 int npk , verb=1 , win=0 , his=0 , fit=1 , do2=0 ;
00049 char *dname , cmd[2222] , *label=NULL , *fname="Anhist" ;
00050
00051 static int hist[32768] , gist[32768] ;
00052 int qq,ncut,ib,nold , sbot,stop , npos,nhalf , cbot,ctop,nwid ;
00053 double dsum ;
00054 float *wt , ws ;
00055 int ibot,itop ;
00056 FILE *hf ;
00057 float **Gmat, *Hvec, *lam, *rez, sum,wtm, wbot,wtop, ebest=-1.0;
00058 float *ap,*pk,*ww ; int nregtry , nbetter ;
00059 float *pkbest,*wwbest,*apbest,*lambest , pplm,aplm,wplm ;
00060 float *pklast,*wwlast,*aplast ;
00061
00062 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00063 printf("Usage: 3dAnhist [options] dataset\n"
00064 "Input dataset is a T1-weighted high-res of the brain (shorts only).\n"
00065 "Output is a list of peaks in the histogram, to stdout, in the form\n"
00066 " ( datasetname #peaks peak1 peak2 ... )\n"
00067 "In the C-shell, for example, you could do\n"
00068 " set anhist = `3dAnhist -q -w1 dset+orig`\n"
00069 "Then the number of peaks found is in the shell variable $anhist[2].\n"
00070 "\n"
00071 "Options:\n"
00072 " -q = be quiet (don't print progress reports)\n"
00073 " -h = dump histogram data to Anhist.1D and plot to Anhist.ps\n"
00074 " -F = DON'T fit histogram with stupid curves.\n"
00075 " -w = apply a Winsorizing filter prior to histogram scan\n"
00076 " (or -w7 to Winsorize 7 times, etc.)\n"
00077 " -2 = Analyze top 2 peaks only, for overlap etc.\n"
00078 "\n"
00079 " -label xxx = Use 'xxx' for a label on the Anhist.ps plot file\n"
00080 " instead of the input dataset filename.\n"
00081 " -fname fff = Use 'fff' for the filename instead of 'Anhist'.\n"
00082 "\n"
00083 "If the '-2' option is used, AND if 2 peaks are detected, AND if\n"
00084 "the -h option is also given, then stdout will be of the form\n"
00085 " ( datasetname 2 peak1 peak2 thresh CER CJV count1 count2 count1/count2)\n"
00086 "where 2 = number of peaks\n"
00087 " thresh = threshold between peak1 and peak2 for decision-making\n"
00088 " CER = classification error rate of thresh\n"
00089 " CJV = coefficient of joint variation\n"
00090 " count1 = area under fitted PDF for peak1\n"
00091 " count2 = area under fitted PDF for peak2\n"
00092 " count1/count2 = ratio of the above quantities\n"
00093 "NOTA BENE\n"
00094 "---------\n"
00095 "* If the input is a T1-weighted MRI dataset (the usual case), then\n"
00096 " peak 1 should be the gray matter (GM) peak and peak 2 the white\n"
00097 " matter (WM) peak.\n"
00098 "* For the definitions of CER and CJV, see the paper\n"
00099 " Method for Bias Field Correction of Brain T1-Weighted Magnetic\n"
00100 " Resonance Images Minimizing Segmentation Error\n"
00101 " JD Gispert, S Reig, J Pascau, JJ Vaquero, P Garcia-Barreno,\n"
00102 " and M Desco, Human Brain Mapping 22:133-144 (2004).\n"
00103 "* Roughly speaking, CER is the ratio of the overlapping area of the\n"
00104 " 2 peak fitted PDFs to the total area of the fitted PDFS. CJV is\n"
00105 " (sigma_GM+sigma_WM)/(mean_WM-mean_GM), and is a different, ad hoc,\n"
00106 " measurement of how much the two PDF overlap.\n"
00107 "* The fitted PDFs are NOT Gaussians. They are of the form\n"
00108 " f(x) = b((x-p)/w,a), where p=location of peak, w=width, 'a' is\n"
00109 " a skewness parameter between -1 and 1; the basic distribution\n"
00110 " is defined by b(x)=(1-x^2)^2*(1+a*x*abs(x)) for -1 < x < 1.\n"
00111 "\n"
00112 "-- RWCox - November 2004\n"
00113 ) ;
00114 exit(0) ;
00115 }
00116
00117 mainENTRY("3dAnhist main"); machdep(); AFNI_logger("3dAnhist",argc,argv);
00118 PRINT_VERSION("3dAnhist") ;
00119
00120
00121
00122 while( iarg < argc && argv[iarg][0] == '-' ){
00123
00124 if( strcmp(argv[iarg],"-label") == 0 ){
00125 label = argv[++iarg] ;
00126 iarg++ ; continue ;
00127 }
00128
00129 if( strcmp(argv[iarg],"-fname") == 0 ){
00130 fname = argv[++iarg] ;
00131 if( !THD_filename_ok(fname) ){
00132 fprintf(stderr,"** Bad name after -fname!\n"); exit(1);
00133 }
00134 iarg++ ; continue ;
00135 }
00136
00137 if( argv[iarg][1] == 'q' ){ verb=0 ; iarg++ ; continue ; }
00138 if( argv[iarg][1] == 'v' ){ verb++ ; iarg++ ; continue ; }
00139 if( argv[iarg][1] == 'h' ){ his =1 ; iarg++ ; continue ; }
00140 if( argv[iarg][1] == 'F' ){ fit =0 ; iarg++ ; continue ; }
00141 if( argv[iarg][1] == '2' ){ do2 =1 ; iarg++ ; continue ; }
00142 if( argv[iarg][1] == 'w' ){
00143 win = -1 ;
00144 if( isdigit(argv[iarg][2]) )
00145 sscanf(argv[iarg],"-w%d",&win) ;
00146 if( win <= 0 ) win = 1 ;
00147 iarg++ ; continue ;
00148 }
00149
00150 fprintf(stderr,"** ILLEGAL option: %s\n",argv[iarg]) ; exit(1) ;
00151 }
00152
00153
00154
00155
00156 dset = THD_open_dataset(argv[iarg]) ; dname = argv[iarg] ;
00157 if( !ISVALID_DSET(dset) ){ fprintf(stderr,"** CAN'T open dataset\n");exit(1); }
00158 if( DSET_BRICK_TYPE(dset,0) != MRI_short ){
00159 fprintf(stderr,"** ILLEGAL non-short dataset type\n"); exit(1);
00160 }
00161 nx = DSET_NX(dset); ny = DSET_NY(dset); nz = DSET_NZ(dset); nxy = nx*ny; nxyz = nxy*nz;
00162 if( nx < 32 || ny < 32 || nz < 32 ){
00163 fprintf(stderr,"** Dataset dimensions are less than 32x32x32?!\n"); exit(1);
00164 }
00165 if( verb ) fprintf(stderr,"++ Loading dataset %s\n", DSET_BRIKNAME(dset) ) ;
00166 DSET_load(dset) ;
00167 if( !DSET_LOADED(dset) ){ fprintf(stderr,"** CAN'T load dataset\n");exit(1); }
00168
00169 if( label == NULL ) label = dname ;
00170
00171
00172
00173
00174 if( verb ) fprintf(stderr,"++ Forming automask of connected high-intensity voxels\n") ;
00175 mask = THD_automask( dset ) ;
00176 if( mask == NULL ){
00177 fprintf(stderr,"** Mask creation fails for unknown reasons!\n"); exit(1);
00178 }
00179
00180
00181
00182 if( dilate > 0 ){
00183 int nmm=1 ;
00184 if( verb ) fprintf(stderr,"++ Dilating automask %d time%c\n",dilate,(dilate==1)?'.':'s') ;
00185 ii = rint(0.032*nx) ; nmm = MAX(nmm,ii) ;
00186 ii = rint(0.032*ny) ; nmm = MAX(nmm,ii) ;
00187 ii = rint(0.032*nz) ; nmm = MAX(nmm,ii) ;
00188 for( dd=0 ; dd < dilate ; dd++ ){
00189 THD_mask_dilate ( nx,ny,nz , mask, 3 ) ;
00190 THD_mask_fillin_completely( nx,ny,nz , mask, nmm ) ;
00191 }
00192 for( ii=0 ; ii < nxyz ; ii++ ) mask[ii] = !mask[ii] ;
00193 THD_mask_clust( nx,ny,nz, mask ) ;
00194 for( ii=0 ; ii < nxyz ; ii++ ) mask[ii] = !mask[ii] ;
00195 }
00196
00197 nmask = THD_countmask( DSET_NVOX(dset) , mask ) ;
00198 if( nmask == 0 ){
00199 fprintf(stderr,"** No voxels in the automask?!\n");
00200 print_results( label,0,NULL,NULL ) ; exit(1);
00201 }
00202
00203
00204
00205 if( SIhh > 0.0 ){
00206
00207 if( ORIENT_tinystr[dset->daxes->xxorient][0] == 'S' ){
00208 for( ii=0 ; ii < nx ; ii++ )
00209 for( kk=0 ; kk < nz ; kk++ )
00210 for( jj=0 ; jj < ny ; jj++ )
00211 if( mask[ii+jj*nx+kk*nxy] ) goto CP5 ;
00212 CP5:
00213 SIax = 1 ; SIbot = ii + (int)(SIhh/fabs(DSET_DX(dset))+0.5) ; SItop = nx-1 ;
00214 }
00215
00216 else if( ORIENT_tinystr[dset->daxes->xxorient][1] == 'S' ){
00217 for( ii=nx-1 ; ii >= 0 ; ii-- )
00218 for( kk=0 ; kk < nz ; kk++ )
00219 for( jj=0 ; jj < ny ; jj++ )
00220 if( mask[ii+jj*nx+kk*nxy] ) goto CP6 ;
00221 CP6:
00222 SIax = 1 ; SIbot = 0 ; SItop = ii - (int)(SIhh/fabs(DSET_DX(dset))+0.5) ;
00223 }
00224
00225 else if( ORIENT_tinystr[dset->daxes->yyorient][0] == 'S' ){
00226 for( jj=0 ; jj < ny ; jj++ )
00227 for( kk=0 ; kk < nz ; kk++ )
00228 for( ii=0 ; ii < nx ; ii++ )
00229 if( mask[ii+jj*nx+kk*nxy] ) goto CP3 ;
00230 CP3:
00231 SIax = 2 ; SIbot = jj + (int)(SIhh/fabs(DSET_DY(dset))+0.5) ; SItop = ny-1 ;
00232 }
00233
00234 else if( ORIENT_tinystr[dset->daxes->yyorient][1] == 'S' ){
00235 for( jj=ny-1 ; jj >= 0 ; jj-- )
00236 for( kk=0 ; kk < nz ; kk++ )
00237 for( ii=0 ; ii < nx ; ii++ )
00238 if( mask[ii+jj*nx+kk*nxy] ) goto CP4 ;
00239 CP4:
00240 SIax = 2 ; SIbot = 0 ; SItop = jj - (int)(SIhh/fabs(DSET_DY(dset))+0.5) ;
00241 }
00242
00243 else if( ORIENT_tinystr[dset->daxes->zzorient][0] == 'S' ){
00244 for( kk=0 ; kk < nz ; kk++ )
00245 for( jj=0 ; jj < ny ; jj++ )
00246 for( ii=0 ; ii < nx ; ii++ )
00247 if( mask[ii+jj*nx+kk*nxy] ) goto CP1 ;
00248 CP1:
00249 SIax = 3 ; SIbot = kk + (int)(SIhh/fabs(DSET_DZ(dset))+0.5) ; SItop = nz-1 ;
00250 }
00251
00252 else if( ORIENT_tinystr[dset->daxes->zzorient][1] == 'S' ){
00253 for( kk=nz-1 ; kk >= 0 ; kk-- )
00254 for( jj=0 ; jj < ny ; jj++ )
00255 for( ii=0 ; ii < nx ; ii++ )
00256 if( mask[ii+jj*nx+kk*nxy] ) goto CP2 ;
00257 CP2:
00258 SIax = 3 ; SIbot = 0 ; SItop = kk - (int)(SIhh/fabs(DSET_DZ(dset))+0.5) ;
00259 }
00260
00261
00262
00263 if( SIax > 0 && SIbot <= SItop ){
00264 char *cax="xyz" ;
00265 if( verb )
00266 fprintf(stderr,"++ SI clipping mask along axis %c: slices %d..%d\n" ,
00267 cax[SIax-1] , SIbot,SItop ) ;
00268 switch( SIax ){
00269 case 1:
00270 for( ii=SIbot ; ii <= SItop ; ii++ )
00271 for( kk=0 ; kk < nz ; kk++ )
00272 for( jj=0 ; jj < ny ; jj++ ) mask[ii+jj*nx+kk*nxy] = 0 ;
00273 break ;
00274 case 2:
00275 for( jj=SIbot ; jj <= SItop ; jj++ )
00276 for( kk=0 ; kk < nz ; kk++ )
00277 for( ii=0 ; ii < nx ; ii++ ) mask[ii+jj*nx+kk*nxy] = 0 ;
00278 break ;
00279 case 3:
00280 for( kk=SIbot ; kk <= SItop ; kk++ )
00281 for( jj=0 ; jj < ny ; jj++ )
00282 for( ii=0 ; ii < nx ; ii++ ) mask[ii+jj*nx+kk*nxy] = 0 ;
00283 break ;
00284 }
00285 }
00286 }
00287 nmask = THD_countmask( DSET_NVOX(dset) , mask ) ;
00288 if( nmask <= 999 ){
00289 fprintf(stderr,"** Only %d voxels in the automask?!\n",nmask);
00290 print_results( label,0,NULL,NULL); exit(1);
00291 }
00292
00293
00294
00295 if( win ){
00296 THD_3dim_dataset *wset ; float irad ;
00297 irad = (verb) ? 2.5 : -2.5 ;
00298 wset = WINsorize( dset , win , -1,-1 , irad , "winsor" , 0,0 , mask ) ;
00299 DSET_delete(dset) ;
00300 dset = wset ;
00301 }
00302
00303
00304
00305 sar = DSET_ARRAY(dset,0) ;
00306 mar = (short *) malloc( sizeof(short)*nmask ) ;
00307 for( jj=ii=0 ; ii < nxyz ; ii++ )
00308 if( mask[ii] && sar[ii] > 0 ) mar[jj++] = sar[ii] ;
00309 nmask = jj ;
00310 if( verb )
00311 fprintf(stderr,"++ %d voxels in mask [out of %d in dataset]\n",nmask,DSET_NVOX(dset)) ;
00312
00313 DSET_delete(dset) ; free(mask) ;
00314
00315
00316
00317
00318 sbot = stop = mar[0] ;
00319 for( ii=1 ; ii < nmask ; ii++ )
00320 if( mar[ii] < sbot ) sbot = mar[ii] ;
00321 else if( mar[ii] > stop ) stop = mar[ii] ;
00322 if( sbot == stop ){
00323 if( verb ) fprintf(stderr,"+ All voxels inside mask have value=%d ?!\n",sbot);
00324 pval[0] = sbot ; wval[0] = 0 ;
00325 print_results( label,1 , pval,wval ) ; exit(0) ;
00326 }
00327 if( verb > 1 ) fprintf(stderr,"++ Data range: %d .. %d\n",sbot,stop) ;
00328
00329
00330
00331 memset( hist , 0 , sizeof(int)*32768 ) ;
00332 for( ii=0 ; ii < nmask ; ii++ ) hist[mar[ii]]++ ;
00333 dsum = 0.0 ;
00334 for( ii=sbot ; ii <= stop ; ii++ )
00335 dsum += (double)(hist[ii])*(double)(hist[ii])*ii ;
00336
00337
00338
00339 qq = 0.65 * nmask ; ib = rint(0.5*sqrt(dsum/nmask)) ;
00340 for( kk=0,ii=stop ; ii >= ib && kk < qq ; ii-- ) kk += hist[ii] ;
00341 ncut = ii ; qq = 0 ;
00342 do{
00343 for( npos=0,ii=ncut ; ii <= stop ; ii++ ) npos += hist[ii] ;
00344 nhalf = npos/2 ;
00345 for( kk=0,ii=ncut ; ii <= stop && kk < nhalf ; ii++ )
00346 kk += hist[ii] ;
00347 nold = ncut ;
00348 ncut = 0.5 * ii ;
00349 qq++ ;
00350 } while( qq < 20 && ncut != nold ) ;
00351
00352
00353
00354 cbot = ncut ; ctop = 3*ncut ;
00355 if( verb ) fprintf(stderr,"++ Histogram range = %d .. %d\n",cbot,ctop) ;
00356 nwid = rint(0.1*ncut) ;
00357 if( nwid <= 0 ){
00358 if( verb ) fprintf(stderr,"++ Using unsmoothed histogram\n") ;
00359 memcpy( gist , hist , sizeof(int)*32768 ) ;
00360 } else {
00361 if( verb > 1 ) fprintf(stderr,"++ Computing smoothing weights:") ;
00362 ws = 0.0 ;
00363 wt = (float *)malloc(sizeof(float)*(2*nwid+1)) ;
00364 for( ii=0 ; ii <= 2*nwid ; ii++ ){
00365 wt[ii] = nwid-abs(nwid-ii) + 0.5f ;
00366 if( verb > 1 ) fprintf(stderr," (%d,%g)" , ii,wt[ii]) ;
00367 ws += wt[ii] ;
00368 }
00369 if( verb > 1 ) fprintf(stderr,"\n") ;
00370 for( ii=0 ; ii <= 2*nwid ; ii++ ) wt[ii] /= ws ;
00371
00372 if( verb )
00373 fprintf(stderr,"++ Smoothing histogram = %d .. %d around each value\n",
00374 -nwid,nwid) ;
00375
00376 for( jj=cbot ; jj <= ctop ; jj++ ){
00377 ibot = jj-nwid ; if( ibot < sbot ) ibot = sbot ;
00378 itop = jj+nwid ; if( itop > stop ) itop = stop ;
00379 ws = 0.0 ;
00380 for( ii=ibot ; ii <= itop ; ii++ )
00381 ws += wt[nwid-jj+ii] * hist[ii] ;
00382 gist[jj] = rint(ws) ;
00383 if( verb > 1 )
00384 fprintf(stderr," + %3d: unsmoothed=%d smoothed=%d\n",
00385 jj,hist[jj],gist[jj]) ;
00386 }
00387 free(wt) ;
00388 }
00389
00390 if( verb ) fprintf(stderr,"++ Scanning histogram for peaks:" ) ;
00391 npk = 0 ;
00392 for( ii=cbot+2 ; ii <= ctop-2 ; ii++ ){
00393 if( gist[ii] > gist[ii-1] &&
00394 gist[ii] > gist[ii-2] &&
00395 gist[ii] > gist[ii+1] &&
00396 gist[ii] > gist[ii+2] ){
00397 pval[npk]=ii; wval[npk++] = gist[ii];
00398 if( verb ) fprintf(stderr," %.1f",pval[npk-1]) ;
00399 }
00400
00401 else if( gist[ii] == gist[ii+1] &&
00402 gist[ii] > gist[ii-1] &&
00403 gist[ii] > gist[ii-2] &&
00404 gist[ii] > gist[ii+2] ){
00405 pval[npk]=ii+0.5; wval[npk++] = gist[ii];
00406 if( verb ) fprintf(stderr," %.1f",pval[npk-1]) ;
00407 }
00408
00409 else if( gist[ii] == gist[ii+1] &&
00410 gist[ii] == gist[ii-1] &&
00411 gist[ii] > gist[ii-2] &&
00412 gist[ii] > gist[ii+2] ){
00413 pval[npk]=ii; wval[npk++] = gist[ii];
00414 if( verb ) fprintf(stderr," %.1f",pval[npk-1]) ;
00415 }
00416 }
00417 if( verb ) fprintf(stderr,"\n") ;
00418
00419 if( do2 && npk > 2 ){
00420 float pval_top, pval_2nd, wval_top, wval_2nd , www; int iii,itop ;
00421 www = wval[0] ; iii = 0 ;
00422 for( ii=1 ; ii < npk ; ii++ ){
00423 if( wval[ii] > www ){ www = wval[ii] ; iii = ii ; }
00424 }
00425 pval_top = pval[iii] ; wval_top = www ; itop = iii ;
00426 www = -1 ; iii = -1 ;
00427 for( ii=0 ; ii < npk ; ii++ ){
00428 if( ii != itop && wval[ii] > www ){ www = wval[ii] ; iii = ii ; }
00429 }
00430 pval_2nd = pval[iii] ; wval_2nd = www ;
00431
00432
00433
00434 if( pval_top < pval_2nd ){
00435 pval[0] = pval_top ; wval[0] = wval_top ;
00436 pval[1] = pval_2nd ; wval[1] = wval_2nd ;
00437 } else {
00438 pval[0] = pval_2nd ; wval[0] = wval_2nd ;
00439 pval[1] = pval_top ; wval[1] = wval_top ;
00440 }
00441 npk = 2 ;
00442 if( verb )
00443 fprintf(stderr,"++ Keeping top 2 peaks: %.1f %.1f\n",pval[0],pval[1]) ;
00444 }
00445
00446 if( his ){
00447 npos = 0 ;
00448 if( npk > 0 && fit ){
00449 int ndim=ctop-cbot+1 , nvec=npk , iw,nw ;
00450 srand48((long)time(NULL)) ;
00451 Gmat = (float **)malloc(sizeof(float *)*nvec) ;
00452 for( jj=0 ; jj < nvec ; jj++ )
00453 Gmat[jj] = (float *)malloc(sizeof(float)*ndim) ;
00454 Hvec = (float *)malloc(sizeof(float)*ndim) ;
00455 rez = (float *)malloc(sizeof(float)*ndim) ;
00456 wt = (float *)malloc(sizeof(float)*ndim) ;
00457 lam = (float *)malloc(sizeof(float)*nvec) ;
00458 ww = (float *)malloc(sizeof(float)*nvec) ;
00459 pk = (float *)malloc(sizeof(float)*nvec) ;
00460 ap = (float *)malloc(sizeof(float)*nvec) ;
00461 apbest = (float *)malloc(sizeof(float)*nvec) ;
00462 pkbest = (float *)malloc(sizeof(float)*nvec) ;
00463 wwbest = (float *)malloc(sizeof(float)*nvec) ;
00464 lambest= (float *)malloc(sizeof(float)*nvec) ;
00465 aplast = (float *)malloc(sizeof(float)*nvec) ;
00466 pklast = (float *)malloc(sizeof(float)*nvec) ;
00467 wwlast = (float *)malloc(sizeof(float)*nvec) ;
00468
00469 if( verb ) fprintf(stderr,"++ Regressing histogram") ;
00470 wbot=0.1*cbot; wtop=0.9*cbot; pplm=0.05*cbot; aplm=0.75; wplm=0.4*cbot;
00471 nregtry = 0 ;
00472 RegTry:
00473 switch(nvec){
00474 case 1: nw = 30000 ; break ;
00475 case 2: nw = 400000 ; break ;
00476 default: nw = 666000 ; break ;
00477 }
00478 if( nregtry > 0 ){
00479 pplm *= 0.7 ; aplm *= 0.7 ; wplm *= 0.7 ; nw /= 2 ;
00480 memcpy(aplast,apbest,sizeof(float)*nvec) ;
00481 memcpy(pklast,pkbest,sizeof(float)*nvec) ;
00482 memcpy(wwlast,wwbest,sizeof(float)*nvec) ;
00483 } else {
00484 for( jj=0 ; jj < nvec ; jj++ ){
00485 wwlast[jj] = 0.5*cbot ;
00486 aplast[jj] = 0.0 ;
00487 pklast[jj] = pval[jj] ;
00488 }
00489 }
00490
00491 for( ii=0 ; ii < ndim ; ii++ ){
00492 wt[ii] = pow( (double)(1+gist[cbot+ii]) , 1.25 ) ;
00493 }
00494 for( ii=0,sum=0.0 ; ii < ndim ; ii++ ) sum += wt[ii] ;
00495 for( ii=0 ; ii < ndim ; ii++ ) wt[ii] /= sum ;
00496 wtm = 1.0 ;
00497 for( ii=0 ; ii < ndim ; ii++ )
00498 if( wt[ii] < wtm && wt[ii] > 0.0 ) wtm = wt[ii] ;
00499 for( ii=0 ; ii < ndim ; ii++ )
00500 if( wt[ii] == 0.0 ) wt[ii] = wtm ;
00501
00502 nbetter = 0 ;
00503 for( iw=0 ; iw < nw ; iw++ ){
00504 for( jj=0 ; jj < nvec ; jj++ ){
00505 ww[jj] = wwlast[jj] + (2.*drand48()-1.)*wplm ;
00506 pk[jj] = pklast[jj] + (2.*drand48()-1.)*pplm ;
00507 ap[jj] = aplast[jj] + (2.*drand48()-1.)*aplm ;
00508 if( pk[jj] > 0.9*ctop ) pk[jj] = 0.9*ctop ;
00509 else if( pk[jj] < 1.1*cbot ) pk[jj] = 1.1*cbot ;
00510 if( ap[jj] > 1.0 ) ap[jj] = 1.0 ;
00511 else if( ap[jj] < -1.0 ) ap[jj] = -1.0 ;
00512 if( ww[jj] < 0.1*cbot ) ww[jj] = 0.1*cbot ;
00513 else if( ww[jj] > 0.9*cbot ) ww[jj] = 0.9*cbot ;
00514
00515 if( pk[jj]+ww[jj] > ctop ) ww[jj] = ctop-pk[jj] ;
00516 else if( pk[jj]-ww[jj] < cbot ) ww[jj] = pk[jj]-cbot ;
00517 }
00518 sum = wtm = 0.0 ;
00519 for( ii=0 ; ii < ndim ; ii++ ){
00520 for( jj=0 ; jj < nvec ; jj++ ){
00521 apar = ap[jj] ;
00522 Gmat[jj][ii] = pmodel_bin(cbot+ii-0.5,cbot+ii+0.5,pk[jj],ww[jj]) ;
00523 if( Gmat[jj][ii] < 0.0 ) Gmat[jj][ii] = 0.0 ;
00524 }
00525 }
00526 for( ii=0 ; ii < ndim ; ii++ ){
00527 Hvec[ii] = gist[cbot+ii] * wt[ii] ;
00528 for( jj=0 ; jj < nvec ; jj++ ) Gmat[jj][ii] *= wt[ii] ;
00529 }
00530 for( jj=0 ; jj < nvec ; jj++ ) lam[jj] = 1.0 ;
00531 for( ii=0 ; ii < ndim ; ii++ ) rez[ii] = 1.0 ;
00532 sum = cl1_solve_res( ndim,nvec , Hvec,Gmat , lam,1 , rez,1 ) ;
00533 if( sum >= 0.0 ){
00534 if( ebest < 0.0 || sum < ebest ){
00535 ebest = sum ;
00536 for( ii=0 ; ii < ndim ; ii++ ){
00537 sum = 0.0 ;
00538 for( jj=0 ; jj < nvec ; jj++ ) sum += Gmat[jj][ii] * lam[jj] ;
00539 hist[cbot+ii] = (int)rint(sum/wt[ii]) ;
00540 }
00541 npos = 1 ; nbetter++ ;
00542 memcpy( apbest , ap , sizeof(float)*nvec ) ;
00543 memcpy( pkbest , pk , sizeof(float)*nvec ) ;
00544 memcpy( wwbest , ww , sizeof(float)*nvec ) ;
00545 memcpy( lambest, lam, sizeof(float)*nvec ) ;
00546 if( verb ) fprintf(stderr,"+") ;
00547 if( verb > 1 ){
00548 fprintf(stderr,"[") ;
00549 for( jj=0 ; jj < nvec ; jj++ ){
00550 fprintf(stderr,"p=%.1f w=%.2f a=%.3f",
00551 pkbest[jj],wwbest[jj],apbest[jj]) ;
00552 if( jj < nvec-1 ) fprintf(stderr,";") ;
00553 }
00554 fprintf(stderr,"]") ;
00555 }
00556 }
00557 }
00558 }
00559 if( verb ) fprintf(stderr,".") ;
00560 if( nregtry < 8 ){ nregtry++ ; goto RegTry ; }
00561 if( nregtry < 12 && nbetter > 0 ){ nregtry++ ; goto RegTry ; }
00562 if( verb ) fprintf(stderr,"!\n") ;
00563 }
00564
00565
00566
00567 sprintf(cmd,"%s.1D",fname) ;
00568 hf = fopen( cmd , "w" ) ;
00569 if( hf == NULL ){
00570 fprintf(stderr,"** Can't open file %s for output!\n",cmd) ;
00571 } else {
00572 static char pbuf[1024]="\0" ;
00573 fprintf(hf,"# 3dAnhist") ;
00574 for( ii=1 ; ii < argc ; ii++ ) fprintf(hf," %s",argv[ii]) ;
00575 fprintf(hf,"\n") ;
00576 for( jj=0 ; jj < npk ; jj++ ){
00577 fprintf(hf,"# Peak %d: location=%.1f histpeak=%.1f\n",jj+1,pval[jj],wval[jj]) ;
00578 }
00579 if( fit && npos > 0 ){
00580 float gval ; int gtot ;
00581
00582
00583
00584
00585 if( verb > 1 ) fprintf(stderr,"++ Computing individual fits:") ;
00586 for( jj=0 ; jj < npk ; jj++ ) wt[jj] = 0.0 ;
00587 gtot = 0 ;
00588 for( ii=cbot ; ii <= ctop ; ii++ ){
00589 gtot += gist[ii] ;
00590 for( jj=0 ; jj < npk ; jj++ ){
00591 if( verb > 1 ) fprintf(stderr,"(%d,%d)",jj,ii) ;
00592 apar = apbest[jj] ;
00593 gval = lambest[jj] *
00594 pmodel_bin(ii-0.5,ii+0.5,pkbest[jj],wwbest[jj]) ;
00595 Gmat[jj][ii-cbot] = rint(gval) ;
00596 wt[jj] += Gmat[jj][ii-cbot] ;
00597 }
00598 }
00599 if( verb > 1 ) fprintf(stderr,"\n") ;
00600
00601 for( jj=0 ; jj < npk ; jj++ ){
00602 fprintf(hf,"# Peak %d fit: location=%.1f width=%.2f skew=%.3f height=%.1f\n",
00603 jj+1,pkbest[jj],wwbest[jj],apbest[jj],lambest[jj] ) ;
00604 if( jj < 4 ){
00605 ii = strlen(pbuf) ;
00606 sprintf(pbuf+ii," #%d=%.1f\\pm%.1f",jj+1,pkbest[jj],wwbest[jj]) ;
00607 }
00608 }
00609 fprintf(hf,"#\n") ;
00610
00611
00612
00613 if( do2 && npk == 2 ){
00614 float mu_bot,sd_bot , mu_top,sd_top ;
00615 int ibot,itop , ix , gx0,gx1 ;
00616
00617
00618
00619 apar = apbest[0] ;
00620 mu_bot = pmodel_mean(pkbest[0],wwbest[0]) ;
00621 sd_bot = pmodel_sigma(wwbest[0]) ;
00622 fprintf(hf,"# mu_bot=%.1f sd_bot=%.2f\n",mu_bot,sd_bot) ;
00623 apar = apbest[1] ;
00624 mu_top = pmodel_mean(pkbest[1],wwbest[1]) ;
00625 sd_top = pmodel_sigma(wwbest[1]) ;
00626 fprintf(hf,"# mu_top=%.1f sd_top=%.2f\n",mu_top,sd_top) ;
00627 cjv = (sd_bot+sd_top)/fabs(mu_top-mu_bot) ;
00628 fprintf(hf,"# CJV = %.3f%%\n",100.0*cjv) ;
00629
00630 fprintf(hf,"# Histogram = %d voxels\n",gtot) ;
00631 fprintf(hf,"# Fit 1 = %d voxels\n",(int)wt[0]) ;
00632 fprintf(hf,"# Fit 2 = %d voxels\n",(int)wt[1]) ;
00633 fprintf(hf,"# Ratio 1/2 = %.3f\n" ,wt[0]/wt[1]) ;
00634
00635 gmcount = wt[0] ; wmcount = wt[1] ;
00636
00637
00638
00639 ibot = (int)mu_bot; ibot = MAX(ibot,cbot); ibot = MIN(ibot,ctop);
00640 itop = 1+(int)mu_top; itop = MAX(itop,cbot); itop = MIN(itop,ctop);
00641 for( ix=ibot ; ix < itop ; ix++ )
00642 if( Gmat[0][ix-cbot] < Gmat[1][ix-cbot] ) break ;
00643 if( ix < itop && ix > ibot ){
00644 threshold = ix-0.5 ;
00645 gx0 = gx1 = 0 ;
00646 for( ii=cbot ; ii < ix ; ii++ ) gx0 += Gmat[1][ii-cbot] ;
00647 for( ii=ix ; ii <= ctop ; ii++ ) gx1 += Gmat[0][ii-cbot] ;
00648 fprintf(hf,"# Overlaps: Thresh=%.1f Fit 1=%d Fit 2=%d\n",
00649 threshold,gx0,gx1) ;
00650 cer = ((float)(gx0+gx1))/(wt[0]+wt[1]) ;
00651 fprintf(hf,"# CER = %.3f%%\n",100.0*cer) ;
00652 #if 0
00653 } else {
00654 fprintf(hf,"# Overlaps: None\n") ;
00655 threshold = 0.5*(mu_bot+mu_top) ;
00656 cer = MIN(gmcount,wmcount)/(gmcount+wmcount) ;
00657 #endif
00658 }
00659 }
00660
00661 fprintf(hf,"# Histogram Region: min=%d max=%d\n",cbot,ctop) ;
00662
00663 fprintf(hf,"# Val Histog FitSum Hi-Fit") ;
00664 for( jj=0 ; jj < npk ; jj++ ) fprintf(hf," Fit#%1d ",jj+1) ;
00665 fprintf(hf,"\n") ;
00666
00667 fprintf(hf,"# --- ------ ------ ------") ;
00668 for( jj=0 ; jj < npk ; jj++ ) fprintf(hf," ------") ;
00669 fprintf(hf,"\n") ;
00670 for( ii=cbot ; ii <= ctop ; ii++ ){
00671 fprintf(hf,"%5d %6d %6d %6d",ii,gist[ii],hist[ii],gist[ii]-hist[ii]) ;
00672 for( jj=0 ; jj < npk ; jj++ )
00673 fprintf(hf," %6d",(int)Gmat[jj][ii-cbot]) ;
00674 fprintf(hf,"\n") ;
00675 }
00676
00677 #if 0
00678 for( jj=0; jj < npk ; jj++ ) free(Gmat[jj]);
00679 free(Gmat);free(pk);free(ww);free(ap);free(Hvec);free(lam);free(rez);free(wt);
00680 free(apbest);free(wwbest);free(pkbest);free(lambest);
00681 free(aplast);free(wwlast);free(pklast);
00682 #endif
00683 ii = (do2 && npk == 2) ? 5 : 3 ;
00684 sprintf(cmd,
00685 "1dplot -ps -nopush -one -xzero %d -xlabel '%s:%s' '%s.1D[1..%d]' > %s.ps" ,
00686 cbot , label , pbuf , fname,ii,fname ) ;
00687 } else {
00688 fprintf(hf,"# Histogram Region: min=%d max=%d\n",cbot,ctop) ;
00689 fprintf(hf,"# Val Histog\n") ;
00690 fprintf(hf,"# --- ------\n") ;
00691 for( ii=cbot ; ii <= ctop ; ii++ )
00692 fprintf(hf,"%5d %6d\n",ii,gist[ii]) ;
00693 for( jj=0 ; jj < npk && jj < 4 ; jj++ ){
00694 ii = strlen(pbuf) ;
00695 sprintf(pbuf+ii," #1%d=%.1f",jj+1,pval[jj]) ;
00696 }
00697 sprintf(cmd,"1dplot -ps -nopush -xzero %d -xlabel '%s:%s' '%s.1D[1]' > %s.ps",
00698 cbot , label , pbuf , fname,fname ) ;
00699 }
00700 fclose(hf) ;
00701 if( verb ){
00702 fprintf(stderr,"++ %s\n",cmd) ;
00703 fprintf(stderr,"++ To view plot: gv -landscape %s.ps\n",fname) ;
00704 }
00705 system( cmd ) ;
00706 }
00707 }
00708
00709 print_results( label,npk , pval , wval ) ;
00710 exit(0) ;
00711 }
00712
00713
00714
00715 void print_results( char *dname , int np , float *pk , float *wd )
00716 {
00717 int ii ;
00718 printf(" ( %s %d " , dname,np ) ;
00719 for( ii=0 ; ii < np ; ii++ ) printf(" %.2f ",pk[ii]) ;
00720
00721 if( threshold > 0.0 ){
00722 printf(" %.2f %.4f %.4f %.0f %.0f %.3f",
00723 threshold , cer , cjv , gmcount , wmcount , gmcount/wmcount ) ;
00724 }
00725
00726 printf(" )\n") ;
00727 }
