1 #include "mrilib.h"
2 
3 /*-------------------------------------------------------------------
4    Fill in the blanks (zeros) in a row of shorts that are at most
5    maxgap in length, and are bounded by the same value.
6    Return value is number of filled-in blanks.
7 ---------------------------------------------------------------------*/
8 
THD_rowfillin_short(int nrow,short * row,int maxgap)9 static int THD_rowfillin_short( int nrow , short * row , int maxgap )
10 {
11    int ii , nfill=0 , jj ;
12    short vbot ;
13 
14    /*-- skip zeros at start --*/
15 
16    for( ii=0 ; ii < nrow && row[ii] == 0 ; ii++ ) ; /* nada */
17    if( ii == nrow ) return nfill ;                    /*** was all zeros ***/
18 
19    /*-- row[ii] is not zero here; now find a gap at or above this --*/
20 
21 #if 0
22 { int kk ; fprintf(stderr,"maxgap=%d: row in:",maxgap);
23   for( kk=0 ; kk < nrow ; kk++ ) fprintf(stderr," %d",row[kk]) ;
24   fprintf(stderr,"\n") ; }
25 #endif
26 
27    while( 1 ){
28       /*-- find next zero value (start of blank region) --*/
29 
30       for( ; ii < nrow && row[ii] != 0 ; ii++ ) ; /* nada */
31       if( ii == nrow ) return nfill ;          /*** didn't find any zero ***/
32 
33 #if 0
34 fprintf(stderr,"  gap start at %d\n",ii) ;
35 #endif
36 
37       vbot = row[ii-1] ;                          /* value just before gap */
38 
39       /*-- find next nonzero value above this zero --*/
40 
41       for( jj=ii+1 ; jj < nrow && row[jj] == 0 ; jj++ ) ; /* nada */
42       if( jj == nrow ) return nfill ;      /*** was all zeros to the end ***/
43 
44 #if 0
45 fprintf(stderr,"  gap end at %d: vbot=%d row[jj]=%d\n",jj,vbot,row[jj]) ;
46 #endif
47 
48       if( row[jj] == vbot && jj-ii <= maxgap ){ /*** fill in this gap!!! ***/
49          nfill += (jj-ii) ;
50          for( ; ii < jj ; ii++ ) row[ii] = vbot ;
51 #if 0
52 { int kk ; fprintf(stderr,"filled:");
53   for( kk=0 ; kk < nrow ; kk++ ) fprintf(stderr," %d",row[kk]) ;
54   fprintf(stderr,"\n") ; }
55 #endif
56       }
57 
58       ii = jj ;                      /* loop back and look for another gap */
59 
60    } /* endless loop */
61 }
62 
63 /*------------------------------------------------------------------------*/
64 
THD_rowfillin_byte(int nrow,byte * row,int maxgap)65 static int THD_rowfillin_byte( int nrow , byte * row , int maxgap )
66 {
67    int ii , nfill=0 , jj ;
68    byte vbot ;
69 
70    /*-- skip zeros --*/
71 
72    for( ii=0 ; ii < nrow && row[ii] == 0 ; ii++ ) ; /* nada */
73    if( ii == nrow ) return nfill ;                    /*** was all zeros ***/
74 
75    /*-- row[ii] is not zero here; now find a gap at or above this --*/
76 
77    while( 1 ){
78       /*-- find next zero value --*/
79 
80       for( ; ii < nrow && row[ii] != 0 ; ii++ ) ; /* nada */
81       if( ii == nrow ) return nfill ;          /*** didn't find any zero ***/
82 
83       vbot = row[ii-1] ;                          /* value at start of gap */
84 
85       /*-- find next nonzero value above this zero --*/
86 
87       for( jj=ii+1 ; jj < nrow && row[jj] == 0 ; jj++ ) ; /* nada */
88       if( jj == nrow ) return nfill ;      /*** was all zeros to the end ***/
89 
90       if( row[jj] == vbot && jj-ii <= maxgap ){ /*** fill in this gap!!! ***/
91          nfill += (jj-ii) ;
92          for( ; ii < jj ; ii++ ) row[ii] = vbot ;
93       }
94 
95       ii = jj ;                      /* loop back and look for another gap */
96 
97    } /* endless loop */
98 }
99 
100 /*------------------------------------------------------------------------*/
101 
THD_rowfillin_float(int nrow,float * row,int maxgap)102 static int THD_rowfillin_float( int nrow , float * row , int maxgap )
103 {
104    int ii , nfill=0 , jj ;
105    float vbot ;
106 
107    /*-- skip zeros --*/
108 
109    for( ii=0 ; ii < nrow && row[ii] == 0 ; ii++ ) ; /* nada */
110    if( ii == nrow ) return nfill ;                    /*** was all zeros ***/
111 
112    /*-- row[ii] is not zero here; now find a gap at or above this --*/
113 
114    while( 1 ){
115       /*-- find next zero value --*/
116 
117       for( ; ii < nrow && row[ii] != 0 ; ii++ ) ; /* nada */
118       if( ii == nrow ) return nfill ;          /*** didn't find any zero ***/
119 
120       vbot = row[ii-1] ;                          /* value at start of gap */
121 
122       /*-- find next nonzero value above this zero --*/
123 
124       for( jj=ii+1 ; jj < nrow && row[jj] == 0 ; jj++ ) ; /* nada */
125       if( jj == nrow ) return nfill ;      /*** was all zeros to the end ***/
126 
127       if( row[jj] == vbot && jj-ii <= maxgap ){ /*** fill in this gap!!! ***/
128          nfill += (jj-ii) ;
129          for( ; ii < jj ; ii++ ) row[ii] = vbot ;
130       }
131 
132       ii = jj ;                      /* loop back and look for another gap */
133 
134    } /* endless loop */
135 }
136 
137 /*---------------------------------------------------------------------------
138    Do the fillin thing on each 1D row from a dataset sub-brick.
139    Return value is number of voxels filled in (-1 if an error transpired).
140 -----------------------------------------------------------------------------*/
141 
THD_dataset_rowfillin(THD_3dim_dataset * dset,int ival,int dcode,int maxgap)142 int THD_dataset_rowfillin( THD_3dim_dataset * dset, int ival, int dcode, int maxgap )
143 {
144    int kind , xx,yy,zz , nrow , nx,ny,nz ;
145    int xtop,ytop,ztop , nff , nfftot=0 ;
146 
147 ENTRY("THD_dataset_rowfillin") ;
148 
149    if( !ISVALID_DSET(dset)      ||
150        ival < 0                 ||
151        ival >= DSET_NVALS(dset) ||
152        maxgap < 1                 ) RETURN(-1) ;  /* bad things */
153 
154    kind = DSET_BRICK_TYPE(dset,ival) ;
155    if( kind != MRI_short && kind != MRI_byte && kind != MRI_float ) RETURN(-1) ;  /* bad */
156 
157    nrow = THD_get_dset_rowcount( dset , dcode ) ;
158    if( nrow < 1 ) RETURN(-1) ;  /* bad */
159 
160    nx = DSET_NX(dset) ;
161    ny = DSET_NY(dset) ;
162    nz = DSET_NZ(dset) ;
163 
164    xtop = ytop = ztop = 1 ;
165    switch( dcode ){
166       case 1: case -1: ytop=ny ; ztop=nz ; break ;
167       case 2: case -2: xtop=nx ; ztop=nz ; break ;
168       case 3: case -3: xtop=nx ; ytop=ny ; break ;
169    }
170 
171    switch( kind ){
172 
173       case MRI_short:{
174          short * row ;
175          for( zz=0 ; zz < ztop ; zz++ )
176           for( yy=0 ; yy < ytop ; yy++ )
177             for( xx=0 ; xx < xtop ; xx++ ){
178                row = THD_get_dset_row( dset,ival , dcode,xx,yy,zz ) ;
179                nff = THD_rowfillin_short( nrow , row , maxgap ) ;
180                if( nff > 0 ){
181                   THD_put_dset_row( dset,ival , dcode,xx,yy,zz , row ) ;
182                   nfftot += nff ;
183                }
184                free(row) ;
185             }
186       }
187       break ;
188 
189       case MRI_byte:{
190          byte * row ;
191          for( zz=0 ; zz < ztop ; zz++ )
192           for( yy=0 ; yy < ytop ; yy++ )
193             for( xx=0 ; xx < xtop ; xx++ ){
194                row = THD_get_dset_row( dset,ival , dcode,xx,yy,zz ) ;
195                nff = THD_rowfillin_byte( nrow , row , maxgap ) ;
196                if( nff > 0 ){
197                   THD_put_dset_row( dset,ival , dcode,xx,yy,zz , row ) ;
198                   nfftot += nff ;
199                }
200                free(row) ;
201             }
202       }
203       break ;
204 
205       case MRI_float:{
206          float * row ;
207          for( zz=0 ; zz < ztop ; zz++ )
208           for( yy=0 ; yy < ytop ; yy++ )
209             for( xx=0 ; xx < xtop ; xx++ ){
210                row = THD_get_dset_row( dset,ival , dcode,xx,yy,zz ) ;
211                nff = THD_rowfillin_float( nrow , row , maxgap ) ;
212                if( nff > 0 ){
213                   THD_put_dset_row( dset,ival , dcode,xx,yy,zz , row ) ;
214                   nfftot += nff ;
215                }
216                free(row) ;
217             }
218       }
219       break ;
220 
221    }
222 
223    RETURN(nfftot) ;
224 }
225