Doxygen Source Code Documentation
mri_medianfilter.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Defines | |
#define | SKIPVOX(i) ( mask != NULL && !mask[i] ) |
Functions | |
MRI_IMAGE * | mri_medianfilter (MRI_IMAGE *imin, float irad, byte *mask, int verb) |
Define Documentation
|
Definition at line 3 of file mri_medianfilter.c. Referenced by mri_medianfilter(). |
Function Documentation
|
Compute the median filter of an input image. Output image is always in float format. ------------------------------------------------------------------------- Definition at line 10 of file mri_medianfilter.c. References ADDTO_CLUSTER, ENTRY, fout, free, MCW_cluster::i, MCW_cluster::j, MCW_cluster::k, KILL_CLUSTER, MRI_IMAGE::kind, malloc, MCW_build_mask(), mri_data_pointer(), MRI_FLOAT_PTR, mri_free(), mri_new_conforming, mri_to_float(), MCW_cluster::num_pt, MRI_IMAGE::nx, MRI_IMAGE::ny, nz, MRI_IMAGE::nz, qmed_float(), RETURN, SKIPVOX, and THD_countmask(). Referenced by main().
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 } |