00001 #include "mrilib.h"
00002 
00003 /*-------------------------------------------------------------------
00004    Fill in the blanks (zeros) in a row of shorts that are at most
00005    maxgap in length, and are bounded by the same value.
00006    Return value is number of filled-in blanks.
00007 ---------------------------------------------------------------------*/
00008 
00009 static int THD_rowfillin_short( int nrow , short * row , int maxgap )
00010 {
00011    int ii , nfill=0 , jj ;
00012    short vbot ;
00013 
00014    /*-- skip zeros at start --*/
00015 
00016    for( ii=0 ; ii < nrow && row[ii] == 0 ; ii++ ) ; /* nada */
00017    if( ii == nrow ) return nfill ;                    /*** was all zeros ***/
00018 
00019    /*-- row[ii] is not zero here; now find a gap at or above this --*/
00020 
00021 #if 0
00022 { int kk ; fprintf(stderr,"maxgap=%d: row in:",maxgap);
00023   for( kk=0 ; kk < nrow ; kk++ ) fprintf(stderr," %d",row[kk]) ;
00024   fprintf(stderr,"\n") ; }
00025 #endif
00026 
00027    while( 1 ){
00028       /*-- find next zero value (start of blank region) --*/
00029 
00030       for( ; ii < nrow && row[ii] != 0 ; ii++ ) ; /* nada */
00031       if( ii == nrow ) return nfill ;          /*** didn't find any zero ***/
00032 
00033 #if 0
00034 fprintf(stderr,"  gap start at %d\n",ii) ;
00035 #endif
00036 
00037       vbot = row[ii-1] ;                          /* value just before gap */
00038 
00039       /*-- find next nonzero value above this zero --*/
00040 
00041       for( jj=ii+1 ; jj < nrow && row[jj] == 0 ; jj++ ) ; /* nada */
00042       if( jj == nrow ) return nfill ;      /*** was all zeros to the end ***/
00043 
00044 #if 0
00045 fprintf(stderr,"  gap end at %d: vbot=%d row[jj]=%d\n",jj,vbot,row[jj]) ;
00046 #endif
00047 
00048       if( row[jj] == vbot && jj-ii <= maxgap ){ /*** fill in this gap!!! ***/
00049          nfill += (jj-ii) ;
00050          for( ; ii < jj ; ii++ ) row[ii] = vbot ;
00051 #if 0
00052 { int kk ; fprintf(stderr,"filled:");
00053   for( kk=0 ; kk < nrow ; kk++ ) fprintf(stderr," %d",row[kk]) ;
00054   fprintf(stderr,"\n") ; }
00055 #endif
00056       }
00057 
00058       ii = jj ;                      /* loop back and look for another gap */
00059 
00060    } /* endless loop */
00061 }
00062 
00063 /*------------------------------------------------------------------------*/
00064 
00065 static int THD_rowfillin_byte( int nrow , byte * row , int maxgap )
00066 {
00067    int ii , nfill=0 , jj ;
00068    byte vbot ;
00069 
00070    /*-- skip zeros --*/
00071 
00072    for( ii=0 ; ii < nrow && row[ii] == 0 ; ii++ ) ; /* nada */
00073    if( ii == nrow ) return nfill ;                    /*** was all zeros ***/
00074 
00075    /*-- row[ii] is not zero here; now find a gap at or above this --*/
00076 
00077    while( 1 ){
00078       /*-- find next zero value --*/
00079 
00080       for( ; ii < nrow && row[ii] != 0 ; ii++ ) ; /* nada */
00081       if( ii == nrow ) return nfill ;          /*** didn't find any zero ***/
00082 
00083       vbot = row[ii-1] ;                          /* value at start of gap */
00084 
00085       /*-- find next nonzero value above this zero --*/
00086 
00087       for( jj=ii+1 ; jj < nrow && row[jj] == 0 ; jj++ ) ; /* nada */
00088       if( jj == nrow ) return nfill ;      /*** was all zeros to the end ***/
00089 
00090       if( row[jj] == vbot && jj-ii <= maxgap ){ /*** fill in this gap!!! ***/
00091          nfill += (jj-ii) ;
00092          for( ; ii < jj ; ii++ ) row[ii] = vbot ;
00093       }
00094 
00095       ii = jj ;                      /* loop back and look for another gap */
00096 
00097    } /* endless loop */
00098 }
00099 
00100 /*------------------------------------------------------------------------*/
00101 
00102 static int THD_rowfillin_float( int nrow , float * row , int maxgap )
00103 {
00104    int ii , nfill=0 , jj ;
00105    float vbot ;
00106 
00107    /*-- skip zeros --*/
00108 
00109    for( ii=0 ; ii < nrow && row[ii] == 0 ; ii++ ) ; /* nada */
00110    if( ii == nrow ) return nfill ;                    /*** was all zeros ***/
00111 
00112    /*-- row[ii] is not zero here; now find a gap at or above this --*/
00113 
00114    while( 1 ){
00115       /*-- find next zero value --*/
00116 
00117       for( ; ii < nrow && row[ii] != 0 ; ii++ ) ; /* nada */
00118       if( ii == nrow ) return nfill ;          /*** didn't find any zero ***/
00119 
00120       vbot = row[ii-1] ;                          /* value at start of gap */
00121 
00122       /*-- find next nonzero value above this zero --*/
00123 
00124       for( jj=ii+1 ; jj < nrow && row[jj] == 0 ; jj++ ) ; /* nada */
00125       if( jj == nrow ) return nfill ;      /*** was all zeros to the end ***/
00126 
00127       if( row[jj] == vbot && jj-ii <= maxgap ){ /*** fill in this gap!!! ***/
00128          nfill += (jj-ii) ;
00129          for( ; ii < jj ; ii++ ) row[ii] = vbot ;
00130       }
00131 
00132       ii = jj ;                      /* loop back and look for another gap */
00133 
00134    } /* endless loop */
00135 }
00136 
00137 /*---------------------------------------------------------------------------
00138    Do the fillin thing on each 1D row from a dataset sub-brick.
00139    Return value is number of voxels filled in (-1 if an error transpired).
00140 -----------------------------------------------------------------------------*/
00141 
00142 int THD_dataset_rowfillin( THD_3dim_dataset * dset, int ival, int dcode, int maxgap )
00143 {
00144    int kind , xx,yy,zz , nrow , nx,ny,nz ;
00145    int xtop,ytop,ztop , nff , nfftot=0 ;
00146 
00147 ENTRY("THD_dataset_rowfillin") ;
00148 
00149    if( !ISVALID_DSET(dset)      ||
00150        ival < 0                 ||
00151        ival >= DSET_NVALS(dset) ||
00152        maxgap < 1                 ) RETURN(-1) ;  /* bad things */
00153 
00154    kind = DSET_BRICK_TYPE(dset,ival) ;
00155    if( kind != MRI_short && kind != MRI_byte && kind != MRI_float ) RETURN(-1) ;  /* bad */
00156 
00157    nrow = THD_get_dset_rowcount( dset , dcode ) ;
00158    if( nrow < 1 ) RETURN(-1) ;  /* bad */
00159 
00160    nx = DSET_NX(dset) ;
00161    ny = DSET_NY(dset) ;
00162    nz = DSET_NZ(dset) ;
00163 
00164    xtop = ytop = ztop = 1 ;
00165    switch( dcode ){
00166       case 1: case -1: ytop=ny ; ztop=nz ; break ;
00167       case 2: case -2: xtop=nx ; ztop=nz ; break ;
00168       case 3: case -3: xtop=nx ; ytop=ny ; break ;
00169    }
00170 
00171    switch( kind ){
00172 
00173       case MRI_short:{
00174          short * row ;
00175          for( zz=0 ; zz < ztop ; zz++ )
00176           for( yy=0 ; yy < ytop ; yy++ )
00177             for( xx=0 ; xx < xtop ; xx++ ){
00178                row = THD_get_dset_row( dset,ival , dcode,xx,yy,zz ) ;
00179                nff = THD_rowfillin_short( nrow , row , maxgap ) ;
00180                if( nff > 0 ){
00181                   THD_put_dset_row( dset,ival , dcode,xx,yy,zz , row ) ;
00182                   nfftot += nff ;
00183                }
00184                free(row) ;
00185             }
00186       }
00187       break ;
00188 
00189       case MRI_byte:{
00190          byte * row ;
00191          for( zz=0 ; zz < ztop ; zz++ )
00192           for( yy=0 ; yy < ytop ; yy++ )
00193             for( xx=0 ; xx < xtop ; xx++ ){
00194                row = THD_get_dset_row( dset,ival , dcode,xx,yy,zz ) ;
00195                nff = THD_rowfillin_byte( nrow , row , maxgap ) ;
00196                if( nff > 0 ){
00197                   THD_put_dset_row( dset,ival , dcode,xx,yy,zz , row ) ;
00198                   nfftot += nff ;
00199                }
00200                free(row) ;
00201             }
00202       }
00203       break ;
00204 
00205       case MRI_float:{
00206          float * row ;
00207          for( zz=0 ; zz < ztop ; zz++ )
00208           for( yy=0 ; yy < ytop ; yy++ )
00209             for( xx=0 ; xx < xtop ; xx++ ){
00210                row = THD_get_dset_row( dset,ival , dcode,xx,yy,zz ) ;
00211                nff = THD_rowfillin_float( nrow , row , maxgap ) ;
00212                if( nff > 0 ){
00213                   THD_put_dset_row( dset,ival , dcode,xx,yy,zz , row ) ;
00214                   nfftot += nff ;
00215                }
00216                free(row) ;
00217             }
00218       }
00219       break ;
00220 
00221    }
00222 
00223    RETURN(nfftot) ;
00224 }