00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #include "mrilib.h"
00008 
00009 /*--------------------------------------------------------------------
00010   routine to edit an input dataset in place according to inputs
00011   in "edopt" (see editvol.h).
00012 
00013   Feb 1996: This routine is much more complex now due to the need to deal
00014             with byte, short, float, or complex data in sub-bricks.
00015 
00016   30 Nov 1997: added ability to edit a given sub-brick, using the
00017                edopt->iv_fim entry
00018 
00019   17 June 1998:  Modifications for erosion and dilation of clusters.
00020 ----------------------------------------------------------------------*/
00021 
00022 void EDIT_one_dataset( THD_3dim_dataset * dset , EDIT_options * edopt )
00023 {
00024    int   edit_thtoin   = edopt->thtoin ;       /* copy into local variables */
00025    int   edit_noneg    = edopt->noneg ;        /* for historical reasons    */
00026    int   edit_abs      = edopt->abss ;
00027    float edit_clip_bot = edopt->clip_bot ;     /* Nov 1995: changed to floats */
00028    float edit_clip_top = edopt->clip_top ;
00029    float edit_thresh   = edopt->thresh ;
00030    int   edit_clust    = edopt->edit_clust ;     /* 10 Sept 1996 */
00031    float clust_rmm     = edopt->clust_rmm ;
00032    float clust_vmul    = edopt->clust_vmul ;
00033    float erode_pv      = edopt->erode_pv;        /* 17 June 1998 */
00034    int   dilate        = edopt->dilate;          /* 17 June 1998 */
00035    int   filter_opt    = edopt->filter_opt;      /* 11 Sept 1996 */
00036    float filter_rmm    = edopt->filter_rmm;      /* 11 Sept 1996 */
00037    int   thrfilter_opt = edopt->thrfilter_opt;   /* 1 Oct 1996 */
00038    float thrfilter_rmm = edopt->thrfilter_rmm;   /* 1 Oct 1996 */
00039    float edit_blur     = edopt->blur ;
00040    float edit_thrblur  = edopt->thrblur;         /* 4 Oct 1996 */
00041    int   edit_scale    = edopt->scale ;
00042    float edit_mult     = edopt->mult ;
00043    int   edit_zvol     = edopt->do_zvol ;
00044    int   edit_ivfim    = edopt->iv_fim ;         /* 30 Nov 1997 */
00045    int   edit_ivthr    = edopt->iv_thr ;         /* 30 Nov 1997 */
00046    int   verbose       = edopt->verbose ;        /* 01 Nov 1999 */
00047    int   fake_dxyz     = edopt->fake_dxyz ;      /* 11 Sep 2000 */
00048 
00049    int   edit_clip_unscaled = edopt->clip_unscaled ;  /* 09 Aug 1996 */
00050 
00051    THD_dataxes   * daxes ;
00052    short   * sfim = NULL , * sthr = NULL ;
00053    float   * ffim = NULL , * fthr = NULL ;
00054    complex * cfim = NULL ;
00055    byte    * bfim = NULL , * bthr = NULL ;
00056    void    * vfim = NULL , * vthr = NULL ;
00057    int nx,ny,nz,nxy,nxyz , jj,kk , ptmin , iclu,nclu , fim_max ;
00058    int iv_fim , iv_thr , fim_type , thr_type ;
00059    register int ii ;
00060    float dx,dy,dz , dxyz , rmm,vmul , val , vvv ;
00061    MCW_cluster_array * clar ;
00062    MCW_cluster       * blur=NULL ;
00063    int fimtype , thrtype ;
00064    float fimfac , thrfac ;
00065 
00066    /** get the data from this dataset **/
00067 
00068 ENTRY("EDIT_one_dataset") ;
00069 
00070    THD_force_malloc_type( dset->dblk , DATABLOCK_MEM_MALLOC ) ;
00071    THD_load_datablock( dset->dblk ) ;
00072 
00073    if( !DSET_LOADED(dset) ){
00074       fprintf(stderr,
00075               "\n*** Cannot read data brick for dataset %s\a\n",
00076               DSET_BRIKNAME(dset) ) ;
00077       EXRETURN ;
00078    }
00079 
00080    /** load the data sub-brick indexes (iv_*) and check types for legality **/
00081 
00082    if( ISANAT(dset) ){
00083       if( edit_ivfim >= 0 && edit_ivfim < DSET_NVALS(dset) )  /* 30 Nov 1997 */
00084          iv_fim = edit_ivfim ;
00085       else
00086          iv_fim = ANAT_ival_zero[dset->func_type] ;
00087 
00088       if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: fim index = %d\n",iv_fim) ;
00089 
00090       fim_type = DSET_BRICK_TYPE(dset,iv_fim) ;
00091       fimfac   = DSET_BRICK_FACTOR(dset,iv_fim) ;
00092       iv_thr   = -1 ;
00093       thr_type = ILLEGAL_TYPE ;
00094 
00095       if( !AFNI_GOOD_DTYPE(fim_type) || fim_type == MRI_rgb ){
00096          fprintf(stderr,"\n*** Illegal anatomy data type %s in dataset %s\a\n" ,
00097                     MRI_type_name[fim_type] ,
00098                     dset->dblk->diskptr->brick_name ) ;
00099          EXRETURN ;
00100       }
00101 
00102 #ifdef AFNI_DEBUG
00103 { char str[256] ;
00104   sprintf(str,"Anat dset: iv=%d type=%s fac=%g",iv_fim,MRI_TYPE_name[fim_type],fimfac) ;
00105   STATUS(str) ; }
00106 #endif
00107 
00108    }
00109 
00110    if( ISFUNC(dset) ){
00111       if( edit_ivfim >= 0 && edit_ivfim < DSET_NVALS(dset) )  /* 30 Nov 1997 */
00112          iv_fim = edit_ivfim ;
00113       else
00114          iv_fim = FUNC_ival_fim[dset->func_type] ;
00115 
00116       if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: fim index = %d\n",iv_fim) ;
00117 
00118       fim_type = DSET_BRICK_TYPE(dset,iv_fim) ;
00119       fimfac   = DSET_BRICK_FACTOR(dset,iv_fim) ;
00120 
00121       if( edit_ivthr >= 0 && edit_ivthr < DSET_NVALS(dset) )  /* 30 Nov 1997 */
00122          iv_thr = edit_ivthr ;
00123       else
00124          iv_thr = FUNC_ival_thr[dset->func_type] ;
00125 
00126       if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: thr index = %d\n",iv_thr) ;
00127 
00128       if( iv_thr < 0 ){
00129          thr_type = ILLEGAL_TYPE ;
00130          thrfac   = 0.0 ;
00131       } else {
00132          thr_type = DSET_BRICK_TYPE(dset,iv_thr) ;
00133          thrfac   = DSET_BRICK_FACTOR(dset,iv_thr) ;
00134          if( thrfac == 0.0 ){
00135             switch( thr_type ){
00136                case MRI_short: thrfac = 1.0/FUNC_scale_short[dset->func_type]; break;
00137                case MRI_byte : thrfac = 1.0/FUNC_scale_byte [dset->func_type]; break;
00138             }
00139          }
00140       }
00141 
00142       if( !AFNI_GOOD_FUNC_DTYPE(fim_type) || fim_type == MRI_rgb ){
00143          fprintf(stderr,"\n*** Illegal functional data type %s in dataset %s\a\n" ,
00144                    MRI_type_name[fim_type], dset->dblk->diskptr->brick_name ) ;
00145          EXRETURN ;
00146       }
00147 
00148       if( thr_type >= 0 && (!AFNI_GOOD_FUNC_DTYPE(thr_type) || fim_type == MRI_rgb) ){
00149          fprintf(stderr,"\n*** Illegal threshold data type %s in dataset %s\a\n" ,
00150                     MRI_type_name[fim_type] , dset->dblk->diskptr->brick_name ) ;
00151          EXRETURN ;
00152       }
00153 
00154 #ifdef AFNI_DEBUG
00155 { char str[256] ;
00156   sprintf(str,"Func dset: iv_fim=%d type=%s fac=%g",iv_fim,MRI_TYPE_name[fim_type],fimfac) ;
00157   STATUS(str) ;
00158   if( iv_thr >= 0 ){
00159   sprintf(str,"Func dset: iv_thr=%d type=%s fac=%g",iv_thr,MRI_TYPE_name[thr_type],thrfac) ;
00160   STATUS(str) ; } }
00161 #endif
00162 
00163    }
00164 
00165    /** load the pointers to the sub-bricks **/
00166 
00167    vfim = DSET_ARRAY(dset,iv_fim) ;
00168    switch( fim_type ){
00169       default:
00170          fprintf(stderr,"\n*** Illegal data type in dataset %s\a\n",
00171                  dset->dblk->diskptr->brick_name ) ;
00172       EXRETURN ;
00173 
00174       case MRI_short:   sfim = (short *)   vfim ; break ;
00175       case MRI_float:   ffim = (float *)   vfim ; break ;
00176       case MRI_byte:    bfim = (byte *)    vfim ; break ;
00177       case MRI_complex: cfim = (complex *) vfim ; break ;
00178    }
00179 
00180    if( iv_thr >= 0 ){
00181       vthr = DSET_ARRAY(dset,iv_thr) ;
00182       switch( thr_type ){
00183          default:
00184             fprintf(stderr,"\n*** Illegal thresh data type in dataset %s\a\n",
00185                     dset->dblk->diskptr->brick_name ) ;
00186          EXRETURN ;
00187 
00188          case MRI_short:   sthr = (short *) vthr ; break ;
00189          case MRI_float:   fthr = (float *) vthr ; break ;
00190          case MRI_byte:    bthr = (byte *)  vthr ; break ;
00191       }
00192    }
00193 
00194    /** load the grid parameters **/
00195 
00196    daxes = dset->daxes ;
00197    nx    = daxes->nxx ; dx = fabs(daxes->xxdel) ;
00198    ny    = daxes->nyy ; dy = fabs(daxes->yydel) ;
00199    nz    = daxes->nzz ; dz = fabs(daxes->zzdel) ;
00200 
00201    if( fake_dxyz ) dx = dy = dz = 1.0 ;  /* 11 Sep 2000 */
00202 
00203    nxy = nx * ny ; nxyz = nxy * nz ; dxyz = dx*dy*dz ;
00204 
00205    /*----- copy threshold over intensity? -----*/
00206 
00207 STATUS("dataset loaded") ;
00208 
00209    if( edit_thtoin && iv_thr >= 0 ){
00210       float new_fimfac , scaling ;
00211 
00212       if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: copy thr over fim\n") ;
00213 
00214       /****
00215             Find scaling factors for various conversions (0 --> no scaling)
00216             scaling    = factor to actually scale data by when copying to new brick
00217             new_fimfac = factor to later scale data by when converting to floats
00218       ****/
00219 
00220       if( edit_thtoin == 2 ){
00221          new_fimfac = scaling = 0.0 ;  /** -2thtoin --> no scaling **/
00222       } else {
00223          switch( thr_type ){
00224 
00225          /** threshold datum is shorts **/
00226 
00227            case MRI_short:{
00228               switch( fim_type ){
00229                  case MRI_short:   /* fim datum is shorts --> no new scaling needed */
00230                     new_fimfac = thrfac ;
00231                     scaling    = 0.0 ;
00232                  break ;
00233 
00234                  case MRI_float:   /* fim datum is floats --> will be scaled properly */
00235                     new_fimfac = 0.0 ;
00236                     scaling    = thrfac ;
00237                  break ;
00238 
00239                  case MRI_byte:    /* fim datum is bytes */
00240                     new_fimfac = 1.0 / FUNC_scale_byte[dset->func_type] ;
00241                     scaling    = thrfac * FUNC_scale_byte[dset->func_type] ;
00242                  break ;
00243               }
00244            }
00245            break ;
00246 
00247            /** threshold datum is bytes **/
00248 
00249            case MRI_byte:{
00250               switch( fim_type ){
00251                  case MRI_short:   /* fim datum is shorts */
00252                     new_fimfac = 1.0 / FUNC_scale_short[dset->func_type] ;
00253                     scaling    = thrfac * FUNC_scale_short[dset->func_type] ;
00254                  break ;
00255 
00256                  case MRI_float:   /* fim datum is floats */
00257                     new_fimfac = 0.0 ;
00258                     scaling    = thrfac ;
00259                  break ;
00260 
00261                  case MRI_byte:    /* fim datum is bytes */
00262                     new_fimfac = thrfac ;
00263                     scaling    = 0.0 ;
00264                  break ;
00265               }
00266            }
00267            break ;
00268 
00269            /** threshold datum is floats **/
00270 
00271            case MRI_float:{
00272               switch( fim_type ){
00273                  case MRI_short:  /* fim datum is shorts */
00274                     new_fimfac = 1.0 / FUNC_scale_short[dset->func_type] ;
00275                     scaling    = FUNC_scale_short[dset->func_type] ;
00276                  break ;
00277 
00278                  case MRI_float:  /* fim datum is floats --> no scaling needed */
00279                     new_fimfac = 0.0 ;
00280                     scaling    = 0.0 ;
00281                  break ;
00282 
00283                  case MRI_byte:   /* fim datum is bytes */
00284                     new_fimfac = 1.0 / FUNC_scale_byte[dset->func_type] ;
00285                     scaling    = FUNC_scale_byte[dset->func_type] ;
00286                  break ;
00287               }
00288            }
00289            break ;
00290         }
00291       }
00292 
00293 #ifdef AFNI_DEBUG
00294 { char str[256] ;
00295   sprintf(str,"thtoin: scaling=%f new_fimfac=%f input=%s output=%s",
00296           scaling,new_fimfac,MRI_TYPE_name[thr_type],MRI_TYPE_name[fim_type]) ;
00297   STATUS(str) ; }
00298 #endif
00299 
00300       /** have scaling factors, so use them **/
00301 
00302       EDIT_coerce_scale_type( nxyz , scaling ,
00303                               thr_type , vthr , fim_type , vfim ) ;
00304 
00305       DSET_BRICK_FACTOR(dset,iv_fim) = fimfac = new_fimfac ;
00306    } /* end -1thtoin */
00307 
00308    /*----- non-negative? -----*/
00309 
00310    if( edit_noneg ){   /* meaningless for byte and complex */
00311 STATUS("noneg") ;
00312 
00313       if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: remove negative values\n") ;
00314 
00315       switch( fim_type ){
00316          case MRI_short:
00317             for( ii=0 ; ii < nxyz ; ii++ ) if( sfim[ii] < 0 ) sfim[ii] = 0 ;
00318          break ;
00319 
00320          case MRI_float:
00321             for( ii=0 ; ii < nxyz ; ii++ ) if( ffim[ii] < 0 ) ffim[ii] = 0 ;
00322          break ;
00323 
00324          default:
00325 STATUS("noneg applied to meaningless type: will be ignored") ;
00326       }
00327    }
00328 
00329    /*----- absolute? -----*/
00330 
00331    if( edit_abs ){   /* meaningless for byte */
00332 STATUS("abs") ;
00333 
00334       if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: take absolute value\n") ;
00335 
00336       switch( fim_type ){
00337          case MRI_short:
00338             for( ii=0 ; ii < nxyz ; ii++ ) if( sfim[ii] < 0 ) sfim[ii] = -sfim[ii] ;
00339          break ;
00340 
00341          case MRI_float:
00342             for( ii=0 ; ii < nxyz ; ii++ ) if( ffim[ii] < 0 ) ffim[ii] = -ffim[ii] ;
00343          break ;
00344 
00345          case MRI_complex:
00346             for( ii=0 ; ii < nxyz ; ii++ ){
00347                cfim[ii].r = CABS(cfim[ii]) ; cfim[ii].i = 0.0 ;
00348             }
00349          break ;
00350 
00351          default:
00352 STATUS("abs applied to meaningless type: will be ignored") ;
00353       }
00354    }
00355 
00356    /*----- clip? -----*/
00357 
00358    if( edit_clip_bot < edit_clip_top ){
00359 
00360       if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: clip fim values\n") ;
00361 
00362       switch( fim_type ){
00363          case MRI_short:{
00364             int top , bot ;
00365             float ftop,fbot ;
00366             if( fimfac > 0.0 && ! edit_clip_unscaled ){
00367                ftop = edit_clip_top / fimfac ;
00368                fbot = edit_clip_bot / fimfac ;
00369             } else {
00370                ftop = edit_clip_top ;
00371                fbot = edit_clip_bot ;
00372             }
00373 
00374             top = rint(ftop) ;  /* this code was modifed 28 Sep 1998 */
00375             if( top >=  MRI_maxshort ) top =   MRI_maxshort + 1  ;
00376             if( top <= -MRI_maxshort ) top = -(MRI_maxshort + 1) ;
00377 
00378             bot = rint(fbot) ;
00379             if( bot >=  MRI_maxshort ) bot =   MRI_maxshort + 1  ;
00380             if( bot <= -MRI_maxshort ) bot = -(MRI_maxshort + 1) ;
00381 
00382 #ifdef AFNI_DEBUG
00383 { char str[256] ;
00384   sprintf(str,"clipping short from %d to %d",bot,top) ;
00385   STATUS(str) ; }
00386 #endif
00387             for( ii=0 ; ii < nxyz ; ii++ )
00388                if( sfim[ii] > bot && sfim[ii] < top ) sfim[ii] = 0 ;
00389          }
00390          break ;
00391 
00392          case MRI_byte:{
00393             int top , bot ;
00394             float ftop,fbot ;
00395             if( fimfac > 0.0 && ! edit_clip_unscaled ){
00396                ftop = edit_clip_top / fimfac ;
00397                fbot = edit_clip_bot / fimfac ;
00398             } else {
00399                ftop = edit_clip_top ;
00400                fbot = edit_clip_bot ;
00401             }
00402 
00403             top = rint(ftop) ;
00404             if( top >=  MRI_maxbyte ) top =   MRI_maxbyte + 1  ;
00405             if( top <= -MRI_maxbyte ) top = -(MRI_maxbyte + 1) ;
00406 
00407             bot = rint(fbot) ;
00408             if( bot >=  MRI_maxbyte ) bot =   MRI_maxbyte + 1  ;
00409             if( bot <= -MRI_maxbyte ) bot = -(MRI_maxbyte + 1) ;
00410 
00411             if( bot < 0 )   bot = 0 ;
00412             if( top < bot ) top = bot ;
00413 #ifdef AFNI_DEBUG
00414 { char str[256] ;
00415   sprintf(str,"clipping byte from %d to %d",bot,top) ;
00416   STATUS(str) ; }
00417 #endif
00418             for( ii=0 ; ii < nxyz ; ii++ )
00419                if( bfim[ii] > bot && bfim[ii] < top ) bfim[ii] = 0 ;
00420          }
00421          break ;
00422 
00423          case MRI_float:{
00424             float top , bot ;
00425             if( fimfac > 0.0 && ! edit_clip_unscaled ){
00426                top = edit_clip_top / fimfac ;
00427                bot = edit_clip_bot / fimfac ;
00428             } else {
00429                top = edit_clip_top ;
00430                bot = edit_clip_bot ;
00431             }
00432 #ifdef AFNI_DEBUG
00433 { char str[256] ;
00434   sprintf(str,"clipping float from %g to %g",bot,top) ;
00435   STATUS(str) ; }
00436 #endif
00437             for( ii=0 ; ii < nxyz ; ii++ )
00438                if( ffim[ii] > bot && ffim[ii] < top ) ffim[ii] = 0.0 ;
00439          }
00440          break ;
00441 
00442          case MRI_complex:{
00443             float val ;
00444             float top , bot ;
00445             if( fimfac > 0.0 && ! edit_clip_unscaled ){
00446                top = edit_clip_top / fimfac ;
00447                bot = edit_clip_bot / fimfac ;
00448             } else {
00449                top = edit_clip_top ;
00450                bot = edit_clip_bot ;
00451             }
00452 #ifdef AFNI_DEBUG
00453 { char str[256] ;
00454   sprintf(str,"clipping complex from %g to %g",bot,top) ;
00455   STATUS(str) ; }
00456 #endif
00457             for( ii=0 ; ii < nxyz ; ii++ ){
00458                val = CABS(cfim[ii]) ;
00459                if( val > bot && val < top ) cfim[ii].r = cfim[ii].i = 0.0 ;
00460             }
00461          }
00462          break ;
00463       }
00464    }
00465 
00466    /*----- apply threshold? -----*/
00467 
00468    if( edit_thresh > 0.0 && iv_thr >= 0 ){
00469 #ifdef AFNI_DEBUG
00470    int nthresh = 0 ;
00471 #  define THADD (nthresh++)
00472 #else
00473 #  define THADD /* nada */
00474 #endif
00475 
00476       if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: apply threshold\n") ;
00477 
00478       switch( thr_type ){
00479 
00480          /** threshold datum is shorts **/
00481 
00482          case MRI_short:{
00483             short thrplu , thrmin ;
00484             float fplu = edit_thresh / thrfac ;
00485             if( fplu > 32767.0 ){
00486                fprintf(stderr,"\n*** -1thresh out of range: reset to %g\n",
00487                                32767.0 * thrfac ) ;
00488                fplu = 32767.0 ;
00489             }
00490             thrplu = (short) fplu ;
00491             thrmin = -thrplu ;
00492 #ifdef AFNI_DEBUG
00493 { char str[256] ;
00494   sprintf(str,"short threshold = %d\n",(int)thrplu) ; STATUS(str) ; }
00495 #endif
00496             switch( fim_type ){
00497                case MRI_short:   /* fim datum is shorts */
00498                   for( ii=0 ; ii < nxyz ; ii++ )
00499                      if( sthr[ii] < thrplu && sthr[ii] > thrmin ){ sfim[ii] = 0 ; THADD ; }
00500                break ;
00501 
00502                case MRI_byte:    /* fim datum is bytes */
00503                   for( ii=0 ; ii < nxyz ; ii++ )
00504                      if( sthr[ii] < thrplu && sthr[ii] > thrmin ){ bfim[ii] = 0 ; THADD ; }
00505                break ;
00506 
00507                case MRI_float:   /* fim datum is floats */
00508                   for( ii=0 ; ii < nxyz ; ii++ )
00509                      if( sthr[ii] < thrplu && sthr[ii] > thrmin ){ ffim[ii] = 0.0 ; THADD ; }
00510                break ;
00511 
00512                case MRI_complex: /* fim datum is complex */
00513                   for( ii=0 ; ii < nxyz ; ii++ )
00514                      if( sthr[ii] < thrplu && sthr[ii] > thrmin ){
00515                        cfim[ii].r = cfim[ii].i = 0.0 ; THADD ;
00516                      }
00517                break ;
00518             }
00519          }
00520          break ;
00521 
00522          /** threshold datum is bytes **/
00523 
00524          case MRI_byte:{
00525             byte thrplu ;
00526             float fplu = edit_thresh / thrfac ;
00527             if( fplu > 255.0 ){
00528                fprintf(stderr,"\n*** -1thresh out of range: reset to %g\n",
00529                                255.0 * thrfac ) ;
00530                fplu = 255.0 ;
00531             }
00532             thrplu = (byte) fplu ;
00533 #ifdef AFNI_DEBUG
00534 { char str[256] ;
00535   sprintf(str,"byte threshold = %d\n",(int)thrplu) ; STATUS(str) ; }
00536 #endif
00537             switch( fim_type ){
00538                case MRI_short:   /* fim datum is shorts */
00539                   for( ii=0 ; ii < nxyz ; ii++ ) if( bthr[ii] < thrplu ){ sfim[ii] = 0 ; THADD ; }
00540                break ;
00541 
00542                case MRI_byte:    /* fim datum is bytes */
00543                   for( ii=0 ; ii < nxyz ; ii++ ) if( bthr[ii] < thrplu ){ bfim[ii] = 0 ; THADD ; }
00544                break ;
00545 
00546                case MRI_float:   /* fim datum is floats */
00547                   for( ii=0 ; ii < nxyz ; ii++ ) if( bthr[ii] < thrplu ){ ffim[ii] = 0.0 ; THADD ; }
00548                break ;
00549 
00550                case MRI_complex: /* fim datum is complex */
00551                   for( ii=0 ; ii < nxyz ; ii++ )
00552                      if( bthr[ii] < thrplu ){
00553                        cfim[ii].r = cfim[ii].i = 0.0 ; THADD ;
00554                      }
00555                break ;
00556             }
00557          }
00558          break ;
00559 
00560          /** threshold datum is floats **/
00561 
00562          case MRI_float:{
00563             float thrplu , thrmin ;
00564             thrplu = edit_thresh ; if( thrfac > 0.0 ) thrplu /= thrfac ;
00565             thrmin = -thrplu ;
00566 #ifdef AFNI_DEBUG
00567 { char str[256] ;
00568   sprintf(str,"float threshold = %g\n",thrplu) ; STATUS(str) ; }
00569 #endif
00570             switch( fim_type ){
00571                case MRI_short:   /* fim datum is shorts */
00572                   for( ii=0 ; ii < nxyz ; ii++ )
00573                      if( fthr[ii] < thrplu && fthr[ii] > thrmin ){ sfim[ii] = 0 ; THADD ; }
00574                break ;
00575 
00576                case MRI_byte:    /* fim datum is bytes */
00577                   for( ii=0 ; ii < nxyz ; ii++ )
00578                      if( fthr[ii] < thrplu && fthr[ii] > thrmin ){ bfim[ii] = 0 ; THADD ; }
00579                break ;
00580 
00581                case MRI_float:   /* fim datum is floats */
00582                   for( ii=0 ; ii < nxyz ; ii++ )
00583                      if( fthr[ii] < thrplu && fthr[ii] > thrmin ){ ffim[ii] = 0.0 ; THADD ; }
00584                break ;
00585 
00586                case MRI_complex: /* fim datum is complex */
00587                   for( ii=0 ; ii < nxyz ; ii++ )
00588                      if( fthr[ii] < thrplu && fthr[ii] > thrmin ){
00589                        cfim[ii].r = cfim[ii].i = 0.0 ; THADD ;
00590                      }
00591                break ;
00592             }
00593          }
00594          break ;
00595       }
00596 #ifdef AFNI_DEBUG
00597 { char str[256] ;
00598   sprintf(str,"number thresholded to zero = %d",nthresh) ;
00599   STATUS(str) ; }
00600 #endif
00601    }
00602 
00603    /*----- blur? -----*/
00604 
00605    if( edit_blur > 0.0 ){
00606 
00607       if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: blurring fim\n") ;
00608 
00609       EDIT_blur_volume( nx,ny,nz, dx,dy,dz , fim_type,vfim , edit_blur ) ;
00610    }
00611 
00612    /*----- threshold blur? -----*/   /* 4 Oct 1996 */
00613    if(( edit_thrblur > 0.0) && (vthr != NULL) ){
00614 
00615       if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: blurring threshold\n") ;
00616 
00617       EDIT_blur_volume( nx,ny,nz, dx,dy,dz , thr_type,vthr , edit_thrblur ) ;
00618    }
00619 
00620 
00621    /*----- zvol? -----*/
00622 
00623    if( edit_zvol ){
00624       THD_ivec3 iv1 , iv2 ;
00625       int ix1,ix2 , jy1,jy2 , kz1,kz2 , jj,kk ;
00626 
00627       iv1 = THD_3dmm_to_3dind(dset,TEMP_FVEC3(edopt->zv_x1,edopt->zv_y1,edopt->zv_z1));
00628       iv2 = THD_3dmm_to_3dind(dset,TEMP_FVEC3(edopt->zv_x2,edopt->zv_y2,edopt->zv_z2));
00629 
00630       ix1 = iv1.ijk[0] ; ix2 = iv2.ijk[0] ;
00631       jy1 = iv1.ijk[1] ; jy2 = iv2.ijk[1] ;
00632       kz1 = iv1.ijk[2] ; kz2 = iv2.ijk[2] ;
00633 
00634       if( ix1 > ix2 ){ ii=ix1 ; ix1=ix2 ; ix2=ii ; }
00635       if( jy1 > jy2 ){ ii=jy1 ; jy1=jy2 ; jy2=ii ; }
00636       if( kz1 > kz2 ){ ii=kz1 ; kz1=kz2 ; kz2=ii ; }
00637 
00638 #ifdef AFNI_DEBUG
00639 { char str[256] ;
00640   sprintf(str,"edit_zvol: x1=%g x2=%g y1=%g y2=%g z1=%g z2=%g",
00641           edopt->zv_x1,edopt->zv_x2,edopt->zv_y1,edopt->zv_y2,edopt->zv_z1,edopt->zv_z2) ;
00642   STATUS(str) ;
00643   sprintf(str,"         : ix1=%d ix2=%d jy1=%d jy2=%d kz1=%d kz2=%d",
00644           ix1,ix2,jy1,jy2,kz1,kz2) ;
00645   STATUS(str) ; }
00646 #endif
00647 
00648       if( verbose )
00649          fprintf(stderr,"--- EDIT_one_dataset: zeroing indexes [%d,%d]x[%d,%d]x[%d,%d]\n",
00650                  ix1,ix2,jy1,jy2,kz1,kz2 ) ;
00651 
00652       for( kk=kz1 ; kk <= kz2 ; kk++ ){
00653          for( jj=jy1 ; jj <= jy2 ; jj++ ){
00654             switch( fim_type ){
00655                case MRI_short:
00656                   for( ii=ix1 ; ii <= ix2 ; ii++ ) sfim[ii+jj*nx+kk*nxy] = 0 ;
00657                break ;
00658 
00659                case MRI_byte:
00660                   for( ii=ix1 ; ii <= ix2 ; ii++ ) bfim[ii+jj*nx+kk*nxy] = 0 ;
00661                break ;
00662 
00663                case MRI_float:
00664                   for( ii=ix1 ; ii <= ix2 ; ii++ ) ffim[ii+jj*nx+kk*nxy] = 0 ;
00665                break ;
00666 
00667                case MRI_complex:
00668                   for( ii=ix1 ; ii <= ix2 ; ii++ )
00669                      cfim[ii+jj*nx+kk*nxy].r = cfim[ii+jj*nx+kk*nxy].i = 0 ;
00670                break ;
00671             }
00672          }
00673       }
00674    }
00675 
00676    /*----- form clusters? -----*/
00677 
00678    rmm  = clust_rmm ;
00679    vmul = clust_vmul ;
00680 
00681    if( rmm >= 0.0 ){       /* do clustering? */
00682 
00683       MCW_cluster_array * clbig ;
00684       MCW_cluster * cl ;
00685 
00686       if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: clustering with rmm=%g vmul=%g\n",
00687                             rmm,vmul ) ;
00688 
00689      /*----- Erosion and dilation of clusters -----*/   /* 17 June 1998 */
00690      if (erode_pv > 0.0)
00691        MCW_erode_clusters (nx, ny, nz, dx, dy, dz, fim_type, vfim, rmm,
00692                            erode_pv, dilate);
00693 
00694 
00695 STATUS("clustering") ;
00696 
00697       if( vmul >= 0.0 )
00698         ptmin = (int)( vmul / dxyz + 0.99 ) ;
00699       else
00700         ptmin = (int) fabs(vmul) ;  /* 30 Apr 2002 */
00701 
00702       vmul = MAX(1,ptmin) * dxyz ;  /* for use below */
00703 
00704       clar  = MCW_find_clusters( nx,ny,nz , dx,dy,dz , fim_type,vfim , rmm ) ;
00705       nclu  = 0 ;
00706 
00707       if( clar != NULL ){
00708          INIT_CLARR(clbig) ;
00709          for( iclu=0 ; iclu < clar->num_clu ; iclu++ ){
00710             cl = clar->clar[iclu] ;
00711             if( cl->num_pt >= ptmin ){ /* big enough */
00712                ADDTO_CLARR(clbig,cl) ;    /* copy pointer */
00713                clar->clar[iclu] = NULL ;  /* null out original */
00714                nclu++ ;
00715             }
00716          }
00717          DESTROY_CLARR(clar) ;
00718          clar = clbig ;
00719          if( nclu == 0 || clar == NULL || clar->num_clu == 0 ){
00720             printf("*** NO CLUSTERS FOUND ***\n") ;
00721             if( clar != NULL ) DESTROY_CLARR(clar) ;
00722             EXRETURN ;
00723          }
00724          SORT_CLARR(clar) ;
00725       }
00726 
00727       if( nclu == 0 ){  /* no data left */
00728 STATUS("no data left after cluster edit!") ;
00729          DESTROY_CLARR(clar) ;
00730          EXRETURN ;
00731       }
00732 
00733 #ifdef AFNI_DEBUG
00734 { char str[256] ;
00735   sprintf(str,"number clusters = %d",nclu) ; STATUS(str) ; }
00736 #endif
00737 
00738       /*----- edit clusters? -----*/   /* 10 Sept 1996 */
00739       if (edit_clust > ECFLAG_SAME)
00740          EDIT_cluster_array (clar, edit_clust, dxyz, vmul);
00741       if (edit_clust == ECFLAG_SIZE || edit_clust == ECFLAG_ORDER)
00742          DSET_BRICK_FACTOR(dset,iv_fim) = 0.0;
00743 
00744       for( iclu=0 ; iclu < clar->num_clu ; iclu++ )
00745          if( clar->clar[iclu] != NULL && clar->clar[iclu]->num_pt > 0 ){
00746             MCW_cluster_to_vol( nx,ny,nz , fim_type,vfim , clar->clar[iclu] ) ;
00747          } else {
00748          }
00749 
00750       DESTROY_CLARR(clar) ;
00751    }
00752 
00753 
00754    /*----- filter? -----*/   /* 11 Sept 1996 */
00755    if (filter_opt > FCFLAG_NONE){
00756 
00757       if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: filtering fim\n") ;
00758 
00759       EDIT_filter_volume (nx, ny, nz, dx, dy, dz, fim_type, vfim,
00760                           filter_opt, filter_rmm , edopt->fmask , edopt->fexpr );
00761    }
00762 
00763 
00764    /*----- threshold filter? -----*/   /* 1 Oct 1996 */
00765    if ((thrfilter_opt > FCFLAG_NONE) && (vthr != NULL)){
00766 
00767       if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: filtering thr\n") ;
00768 
00769       EDIT_filter_volume (nx, ny, nz, dx, dy, dz, thr_type, vthr,
00770                           thrfilter_opt, thrfilter_rmm , edopt->fmask , edopt->fexpr );
00771    }
00772 
00773 
00774    /*----- scale? -----*/
00775 
00776 #ifdef ALLOW_SCALE_TO_MAX
00777    if( edit_scale ){
00778 STATUS("scale") ;
00779 
00780       if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: scaling fim to max\n") ;
00781 
00782       MCW_scale_to_max( nx,ny,nz , fim_type , vfim ) ;
00783    }
00784 #endif
00785 
00786    /*----- mult? -----*/
00787    /*--- correction for scaling of short and byte bricks (13 Sept. 1996) ---*/
00788 
00789    if( edit_mult != 0.0 ){
00790 STATUS("mult") ;
00791 
00792     if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: multiplying fim\n") ;
00793 
00794      switch( fim_type ){
00795         case MRI_short:
00796            if (fimfac > 0)
00797               DSET_BRICK_FACTOR(dset,iv_fim) =
00798                  DSET_BRICK_FACTOR(dset,iv_fim) * edit_mult ;
00799            else
00800               for( ii=0 ; ii < nxyz ; ii++ ) sfim[ii] *= edit_mult ;
00801         break ;
00802 
00803         case MRI_byte :
00804            if (fimfac > 0)
00805               DSET_BRICK_FACTOR(dset,iv_fim) =
00806                  DSET_BRICK_FACTOR(dset,iv_fim) * edit_mult ;
00807            else
00808               for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] *= edit_mult ;
00809         break ;
00810 
00811         case MRI_float: for( ii=0 ; ii < nxyz ; ii++ ) ffim[ii] *= edit_mult ;
00812         break ;
00813 
00814         case MRI_complex: for( ii=0 ; ii < nxyz ; ii++ )
00815                              cfim[ii].r *= edit_mult , cfim[ii].i *= edit_mult ;
00816         break ;
00817       }
00818    }
00819 
00820    /*----- 17 Sep 1998: conversion to z-score? -----*/
00821 
00822    if( edopt->zscore ){                          /* 07 Jun 1999! How did this get lost? */
00823      int kv = DSET_BRICK_STATCODE(dset,iv_fim) ;
00824      float par[2] ;
00825 
00826      if( FUNC_IS_STAT(kv) && kv != FUNC_ZT_TYPE ){
00827 
00828 #if 0
00829 fprintf(stderr," -1zscore: converting\n") ;
00830 #endif
00831 
00832        if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: converting to zscore\n") ;
00833 
00834         EDIT_zscore_vol( nxyz , fim_type , fimfac , vfim ,
00835                          kv , DSET_BRICK_STATAUX(dset,iv_fim) ) ;
00836 
00837         if( ISBUCKET(dset) ){
00838 
00839 #if 0
00840 fprintf(stderr," -1zscore: bucketing\n") ;
00841 #endif
00842 
00843            par[0] = FUNC_ZT_TYPE ;
00844            par[1] = 0 ;
00845            EDIT_dset_items( dset , ADN_brick_stataux_one+iv_fim,par , ADN_none ) ;
00846 
00847         } else if( ISFUNC(dset)                  &&
00848                    FUNC_IS_STAT(dset->func_type) &&
00849                    iv_fim == FUNC_ival_thr[dset->func_type]  ){
00850 
00851 #if 0
00852 fprintf(stderr," -1zscore: retyping\n") ;
00853 #endif
00854 
00855            dset->func_type   = FUNC_ZT_TYPE ;
00856            dset->stat_aux[0] = 0.0 ;
00857 
00858         } else {
00859            fprintf(stderr,"*** -1zscore error: non-bucket & non-func!\n") ;
00860         }
00861 
00862         if( fim_type == MRI_short )
00863            DSET_BRICK_FACTOR(dset,iv_fim) = 1.0 / FUNC_ZT_SCALE_SHORT ;
00864       }
00865    }
00866 
00867    /*------ DONE! -----*/
00868 
00869    EXRETURN ;
00870 }