Doxygen Source Code Documentation
3dAnhist.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Defines | |
#define | pmodel_bin(x1, x2, pk, wd) (pmodel_cdf((x2-pk)/wd)-pmodel_cdf((x1-pk)/wd)) |
#define | pmodel_mean(pk, wd) ((pk)+0.078125*apar*(wd)) |
#define | pmodel_sigma(wd) (0.002232142857*(wd)*sqrt(28672.-1225.*apar*apar)) |
Functions | |
void | print_results (char *, int, float *, float *) |
double | pmodel_pdf (double x) |
double | pmodel_cdf (double x) |
int | main (int argc, char *argv[]) |
Variables | |
float | threshold = 0.0 |
float | cer = 0.0 |
float | cjv = 0.0 |
float | gmcount = 0.0 |
float | wmcount = 0.0 |
double | apar = 0.0 |
Define Documentation
|
Definition at line 31 of file 3dAnhist.c. Referenced by main(). |
|
Definition at line 33 of file 3dAnhist.c. Referenced by main(). |
|
Definition at line 35 of file 3dAnhist.c. Referenced by main(). |
Function Documentation
|
compute the overall minimum and maximum voxel values for a dataset Definition at line 39 of file 3dAnhist.c. References abs, AFNI_logger(), apar, argc, cer, cjv, cl1_solve_res(), THD_3dim_dataset::daxes, DSET_ARRAY, DSET_BRICK_TYPE, DSET_BRIKNAME, DSET_delete, DSET_DX, DSET_DY, DSET_DZ, DSET_load, DSET_LOADED, DSET_NVOX, DSET_NX, DSET_NY, DSET_NZ, fit, free, gmcount, ISVALID_DSET, machdep(), mainENTRY, malloc, MAX, MIN, nz, pmodel_bin, pmodel_mean, pmodel_sigma, print_results(), PRINT_VERSION, THD_automask(), THD_countmask(), THD_filename_ok(), THD_mask_clust(), THD_mask_dilate(), THD_mask_fillin_completely(), THD_open_dataset(), threshold, win, WINsorize(), wmcount, THD_dataxes::xxorient, THD_dataxes::yyorient, and THD_dataxes::zzorient.
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 /*-- options --*/ 00121 00122 while( iarg < argc && argv[iarg][0] == '-' ){ 00123 00124 if( strcmp(argv[iarg],"-label") == 0 ){ /* 14 Mar 2003 */ 00125 label = argv[++iarg] ; 00126 iarg++ ; continue ; 00127 } 00128 00129 if( strcmp(argv[iarg],"-fname") == 0 ){ /* 14 Mar 2003 */ 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 /*-- read dataset --*/ 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 ; /* 14 Mar 2003 */ 00170 00171 /*-----------------------------*/ 00172 /*** FIND THE BRAIN-ish MASK ***/ 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 /* do fillin, etc, after dilation */ 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 /* find most superior point, clip off all more than SIhh (130) mm below it */ 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 /* cut off stuff below SIhh mm from most Superior point */ 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 /* Winsorize? */ 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 /* copy data in mask region to private array, then expunge dataset */ 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 /*** Mask is finished. Now build histogram of data in mask region. ***/ 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 /* build histogram */ 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 /* find cliplevel */ 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] ; /* number >= cut */ 00344 nhalf = npos/2 ; 00345 for( kk=0,ii=ncut ; ii <= stop && kk < nhalf ; ii++ ) /* find median */ 00346 kk += hist[ii] ; 00347 nold = ncut ; 00348 ncut = 0.5 * ii ; /* new cut */ 00349 qq++ ; 00350 } while( qq < 20 && ncut != nold ) ; 00351 00352 /* make smoothed histogram */ 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] && /* very special case */ 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] && /* super special case */ 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 ){ /* find largest two peaks and keep only them */ 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 /* make sure peaks are increasing in pval */ 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 ){ /* fit histogram? */ 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 ){ /* keep previous best estimates */ 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++ ){ /* initial estimates */ 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++ ){ /* random search nw times */ 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++ ){ /* create basis vectors */ 00520 for( jj=0 ; jj < nvec ; jj++ ){ /* basis for jj-th peak pdf */ 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 ; /* constraints */ 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 } /* end of nw searches */ 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 } /* end of fitting */ 00564 00565 /* save some results */ 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 /* compute individual fits into Gmat[jj] curves, 00583 and total number in each fit into wt[jj] values */ 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 /* 03 Nov 2004: calculate some measures of the overlap */ 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 /* statistics of variation and location */ 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 /* scan for crossover between fitted densities */ 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 { /* 15 Nov 2004 */ 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 } |
|
Definition at line 18 of file 3dAnhist.c.
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 } |
|
skewness parameter [-1..1] * Definition at line 10 of file 3dAnhist.c.
|
|
Definition at line 715 of file 3dAnhist.c. References cer, cjv, gmcount, threshold, and wmcount. Referenced by main().
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 } |
Variable Documentation
|
Definition at line 8 of file 3dAnhist.c. Referenced by main(), pmodel_cdf(), and pmodel_pdf(). |
|
Definition at line 5 of file 3dAnhist.c. Referenced by main(), and print_results(). |
|
Definition at line 5 of file 3dAnhist.c. Referenced by main(), and print_results(). |
|
Definition at line 5 of file 3dAnhist.c. Referenced by main(), and print_results(). |
|
Definition at line 5 of file 3dAnhist.c. Referenced by main(), and print_results(). |
|
Definition at line 5 of file 3dAnhist.c. Referenced by main(), and print_results(). |