00001 #include "mrilib.h"
00002 
00003 #define SKIPVOX(i) ( mask != NULL && !mask[i] )
00004 
00005 /*-----------------------------------------------------------------------*/
00006 /*! Compute the median filter of an input image.
00007     Output image is always in float format.
00008 -------------------------------------------------------------------------*/
00009 
00010 MRI_IMAGE *mri_medianfilter( MRI_IMAGE *imin, float irad, byte *mask, int verb )
00011 {
00012    MRI_IMAGE *imout ;
00013    float *fin , *fout , *tmp ;
00014    short *sin ; byte *bin ; void *vin ;
00015    short *di , *dj , *dk ;
00016    int nd, ii,jj,kk, ip,jp,kp, nx,ny,nz, nxy, ijk, dd,nt,pjk, kd  ;
00017    MCW_cluster *cl ;
00018 
00019 ENTRY("mri_medianfilter") ;
00020 
00021    if( imin == NULL || irad < 1.0 ) RETURN(NULL) ;
00022 
00023    /** if not a good input data type, floatize and try again **/
00024 
00025    if( imin->kind != MRI_float &&
00026        imin->kind != MRI_short &&
00027        imin->kind != MRI_byte    ){
00028 
00029      MRI_IMAGE *qim ;
00030      qim = mri_to_float( imin ) ;
00031      imout = mri_medianfilter( qim , irad , mask , verb ) ;
00032      mri_free( qim ) ;
00033      RETURN(imout) ;
00034    }
00035 
00036    /** build cluster of points for median-izing **/
00037 
00038    if( irad < 1.01 ) irad = 1.01 ;
00039    cl = MCW_build_mask( 0,0,0 , 1.0,1.0,1.0 , irad ) ;
00040 
00041    if( cl == NULL || cl->num_pt < 6 ){ KILL_CLUSTER(cl); RETURN(NULL); }
00042 
00043    ADDTO_CLUSTER(cl,0,0,0,0) ;
00044 
00045    di = cl->i    ; dj = cl->j    ; dk = cl->k    ; nd  = cl->num_pt ;
00046    nx = imin->nx ; ny = imin->ny ; nz = imin->nz ; nxy = nx*ny ;
00047 
00048    if( verb ){
00049      fprintf(stderr,"++ Median mask=%d",nd) ;
00050      if( mask != NULL )
00051        fprintf(stderr," Data mask=%d",THD_countmask(nxy*nz,mask)) ;
00052    }
00053 
00054    imout = mri_new_conforming( imin , MRI_float ) ;
00055    fout  = MRI_FLOAT_PTR( imout ) ;
00056    vin   = mri_data_pointer( imin ) ;
00057    tmp   = (float *) malloc(sizeof(float)*nd) ;
00058    switch( imin->kind ){
00059      case MRI_float:  fin = (float *)vin ; break ;
00060      case MRI_short:  sin = (short *)vin ; break ;
00061      case MRI_byte :  bin = (byte  *)vin ; break ;
00062    }
00063 
00064    if( verb ){ kd = (int)rint(0.03*nz); if( kd < 1 ) kd = 1; }
00065 
00066    for( kk=0 ; kk < nz ; kk++ ){
00067     if( verb && kk%kd == 0 ) fprintf(stderr,".") ;
00068     for( jj=0 ; jj < ny ; jj++ ){
00069       for( ii=0 ; ii < nx ; ii++ ){
00070         ijk = ii + jj*nx + kk*nxy ;
00071         if( SKIPVOX(ijk) ){ fout[ijk] = 0.0; continue; }
00072 
00073         /* extract neighborhood values */
00074 
00075         switch( imin->kind ){
00076           case MRI_float:
00077             for( nt=dd=0 ; dd < nd ; dd++ ){
00078               ip = ii+di[dd] ; if( ip < 0 || ip >= nx ) continue ;
00079               jp = jj+dj[dd] ; if( jp < 0 || jp >= ny ) continue ;
00080               kp = kk+dk[dd] ; if( kp < 0 || kp >= nz ) continue ;
00081               pjk = ip+jp*nx+kp*nxy ;
00082               if( SKIPVOX(pjk) ) continue ;
00083               tmp[nt++] = fin[pjk] ;
00084             }
00085           break ;
00086           case MRI_short:
00087             for( nt=dd=0 ; dd < nd ; dd++ ){
00088               ip = ii+di[dd] ; if( ip < 0 || ip >= nx ) continue ;
00089               jp = jj+dj[dd] ; if( jp < 0 || jp >= ny ) continue ;
00090               kp = kk+dk[dd] ; if( kp < 0 || kp >= nz ) continue ;
00091               pjk = ip+jp*nx+kp*nxy ;
00092               if( SKIPVOX(pjk) ) continue ;
00093               tmp[nt++] = sin[pjk] ;
00094             }
00095           break ;
00096           case MRI_byte:
00097             for( nt=dd=0 ; dd < nd ; dd++ ){
00098               ip = ii+di[dd] ; if( ip < 0 || ip >= nx ) continue ;
00099               jp = jj+dj[dd] ; if( jp < 0 || jp >= ny ) continue ;
00100               kp = kk+dk[dd] ; if( kp < 0 || kp >= nz ) continue ;
00101               pjk = ip+jp*nx+kp*nxy ;
00102               if( SKIPVOX(pjk) ) continue ;
00103               tmp[nt++] = bin[pjk] ;
00104             }
00105           break ;
00106         }
00107 
00108         fout[ijk] = qmed_float( nt , tmp ) ;  /* the actual median-izing */
00109    }}}
00110    if( verb ) fprintf(stderr,"\n") ;
00111 
00112    KILL_CLUSTER(cl); free((void *)tmp);  /* toss the trash */
00113    RETURN(imout) ;
00114 }