1 #include "mrilib.h"
2
3 typedef struct {
4 int nmask[3] ;
5 byte * mask[3] ;
6 } Tmask ;
7
free_Tmask(Tmask * tm)8 void free_Tmask( Tmask * tm )
9 {
10 if( tm != NULL ){
11 free(tm->mask[0]) ; free(tm->mask[1]) ; free(tm->mask[2]) ; free(tm) ;
12 }
13 return ;
14 }
15
16 #define IXY 2 /* fixdir-1 for each plane */
17 #define IYZ 0
18 #define IZX 1
19
create_Tmask(int nx,int ny,int nz,byte * vol)20 Tmask * create_Tmask( int nx, int ny, int nz, byte * vol )
21 {
22 Tmask * tm ;
23 int ii,jj,kk,vv , nxy,nyz,nzx ;
24 byte * bz , *xym,*yzm,*zxm , *bxy,*byz,*bzx ;
25
26 tm = (Tmask *) malloc(sizeof(Tmask)) ;
27 tm->nmask[IXY] = nxy = nx*ny ;
28 tm->nmask[IYZ] = nyz = ny*nz ;
29 tm->nmask[IZX] = nzx = nz*nx ;
30
31 tm->mask[IXY] = xym = (byte *) calloc(1,sizeof(byte)*nxy) ;
32 tm->mask[IYZ] = yzm = (byte *) calloc(1,sizeof(byte)*nyz) ;
33 tm->mask[IZX] = zxm = (byte *) calloc(1,sizeof(byte)*nzx) ;
34
35 for( byz=yzm,kk=0 ; kk < nz ; kk++,byz+=ny ){
36 bz = vol + kk*nxy ;
37 for( bxy=xym,jj=0 ; jj < ny ; jj++,bz+=nx,bxy+=nx ){
38 for( bzx=zxm,ii=0 ; ii < nx ; ii++,bzx+=nz ){
39 if( bz[ii] ){ bxy[ii] = byz[jj] = bzx[kk] = 1 ; }
40 }
41 }
42 }
43
44 return tm ;
45 }
46
47 /*===========================================================================
48 Functions to extract a plane of shifted bytes from a 3D volume.
49 nx, ny, nz = dimensions of vol
50 vol = input 3D volume of bytes
51
52 fixdir = fixed direction (1=x, 2=y, 3=z)
53 fixijk = fixed index
54 da, db = shift in planar coordinaes (non-fixed directions)
55 ma, mb = dimensions of im
56 im = output 2D image
57
58 Goal is im[a,b] = vol[ P(a-da,b-db,c=fixijk) ] for a=0..ma-1, b=0..mb-1,
59 where P(a,b,c) is the permutation of (a,b,c) that goes with fixdir:
60 P(x,y,z) = (y,z,x) for fixdir == 1
61 P(x,y,z) = (z,x,y) for fixdir == 2
62 P(x,y,z) = (x,y,z) for fixdir == 3
63 For values outside the range of vol[], im[] is set to 0.
64
65 The five interpolation routines that follow are:
66 _nn = nearest neigbhor "interpolation"
67 _lifl = linear interpolation, with floating point arithmetic
68 _liby = linear interpolation, with byte arithmetic
69 _ts = two-step interpolation
70 _fs = four-step interpolation
71 =============================================================================*/
72
73 /* macros for offsets in vol[] to corners of the interpolation square */
74
75 #undef LL
76 #undef LR
77 #undef UL
78 #undef UR
79
80 #define LL 0 /* lower left */
81 #define LR astep /* lower right */
82 #define UL bstep /* upper left */
83 #define UR (astep+bstep) /* upper right */
84
85 #define ASSIGN_DIRECTIONS \
86 do{ switch( fixdir ){ \
87 default: \
88 case 1: /* x-direction: (a,b,c) = (y,z,x) */ \
89 astep = nx ; bstep = nxy ; cstep = 1 ; \
90 na = ny ; nb = nz ; nc = nx ; \
91 break ; \
92 \
93 case 2: /* y-direction: (a,b,c) = (z,x,y) */ \
94 astep = nxy ; bstep = 1 ; cstep = nx ; \
95 na = nz ; nb = nx ; nc = ny ; \
96 break ; \
97 \
98 case 3: /* z-direction: (a,b,c) = (x,y,z) */ \
99 astep = 1 ; bstep = nx ; cstep = nxy ; \
100 na = nx ; nb = ny ; nc = nz ; \
101 break ; \
102 } } while(0)
103
104 /*-----------------------------------------------------------------------*/
105
extract_assign_directions(int nx,int ny,int nz,int fixdir,int * Astep,int * Bstep,int * Cstep,int * Na,int * Nb,int * Nc)106 void extract_assign_directions( int nx, int ny, int nz, int fixdir ,
107 int *Astep, int *Bstep, int *Cstep ,
108 int *Na , int *Nb , int *Nc )
109 {
110 int astep,bstep,cstep , na,nb,nc , nxy=nx*ny ;
111
112 ASSIGN_DIRECTIONS ;
113
114 *Astep = astep ; *Bstep = bstep ; *Cstep = cstep ;
115 *Na = na ; *Nb = nb ; *Nc = nc ; return ;
116 }
117
118 /*-----------------------------------------------------------------------
119 NN "interpolation"
120 -------------------------------------------------------------------------*/
121
extract_byte_nn(int nx,int ny,int nz,byte * vol,Tmask * tm,int fixdir,int fixijk,float da,float db,int ma,int mb,byte * im)122 void extract_byte_nn( int nx , int ny , int nz , byte * vol ,
123 Tmask * tm ,
124 int fixdir , int fixijk , float da , float db ,
125 int ma , int mb , byte * im )
126 {
127 int adel,bdel , abot,atop , bb,bbot,btop , nxy=nx*ny ;
128 register int aa , ijkoff , aoff,boff ;
129 int astep,bstep,cstep , na,nb,nc ;
130 byte * mask ;
131
132 memset( im , 0 , ma*mb ) ; /* initialize output to zero */
133
134 if( fixijk < 0 ) return ;
135
136 ASSIGN_DIRECTIONS ;
137
138 if( fixijk >= nc ) return ;
139
140 da += 0.5 ; adel = (int) da ; if( da < 0.0 ) adel-- ; /* floor(da+0.5) */
141 db += 0.5 ; bdel = (int) db ; if( db < 0.0 ) bdel-- ; /* floor(db+0.5) */
142
143 abot = 0 ; if( abot < adel ) abot = adel ; /* range in im[] */
144 atop = na+adel ; if( atop > ma ) atop = ma ;
145
146 bbot = 0 ; if( bbot < bdel ) bbot = bdel ;
147 btop = nb+bdel ; if( btop > mb ) btop = mb ;
148
149 ijkoff = fixijk*cstep + (abot-adel)*astep + (bbot-bdel)*bstep ;
150 boff = bbot * ma ;
151
152 mask = (tm == NULL) ? NULL
153 : tm->mask[fixdir%3] + (fixijk*nb - bdel) ;
154
155 for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
156 if( mask == NULL || mask[bb] )
157 for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
158 im[aa+boff] = vol[aoff+ijkoff] ;
159
160 return ;
161 }
162
163 /*---------------------------------------------------------------------------
164 Linear interpolation with floating point arithmetic
165 -----------------------------------------------------------------------------*/
166
extract_byte_lifl(int nx,int ny,int nz,byte * vol,Tmask * tm,int fixdir,int fixijk,float da,float db,int ma,int mb,byte * im)167 void extract_byte_lifl( int nx , int ny , int nz , byte * vol ,
168 Tmask * tm ,
169 int fixdir , int fixijk , float da , float db ,
170 int ma , int mb , byte * im )
171 {
172 int adel,bdel , abot,atop , bb,bbot,btop , nxy=nx*ny ;
173 register int aa , ijkoff , aoff,boff ;
174 int astep,bstep,cstep , na,nb,nc ;
175 float fa , fb ;
176 float f_a_b , f_ap_b , f_a_bp , f_ap_bp ;
177 byte * mask ;
178
179 memset( im , 0 , ma*mb ) ; /* initialize output to zero */
180
181 if( fixijk < 0 ) return ;
182
183 ASSIGN_DIRECTIONS ;
184
185 if( fixijk >= nc ) return ;
186
187 adel = (int) da ; if( da < 0.0 ) adel-- ; /* floor(da) */
188 bdel = (int) db ; if( db < 0.0 ) bdel-- ; /* floor(db) */
189
190 fa = da - adel ; /* fractional part of dj */
191 fb = db - bdel ; /* fractional part of dk */
192
193 adel++ ; bdel++ ;
194
195 f_a_b = fa * fb ;
196 f_ap_b = (1.0-fa)* fb ;
197 f_a_bp = fa *(1.0-fb) ;
198 f_ap_bp = (1.0-fa)*(1.0-fb) ;
199
200 abot = 0 ; if( abot < adel ) abot = adel ; /* range in im[] */
201 atop = na+adel-1 ; if( atop > ma ) atop = ma ;
202
203 bbot = 0 ; if( bbot < bdel ) bbot = bdel ;
204 btop = nb+bdel-1 ; if( btop > mb ) btop = mb ;
205
206 ijkoff = fixijk*cstep + (abot-adel)*astep + (bbot-bdel)*bstep ;
207 boff = bbot * ma ;
208
209 mask = (tm == NULL) ? NULL
210 : tm->mask[fixdir%3] + (fixijk*nb - bdel) ;
211
212 for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
213 if( mask == NULL || mask[bb] || mask[bb+1] )
214 for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
215 im[aa+boff] = (byte)( f_a_b * vol[aoff+ijkoff]
216 + f_ap_b * vol[aoff+(ijkoff+LR)]
217 + f_a_bp * vol[aoff+(ijkoff+UL)]
218 + f_ap_bp * vol[aoff+(ijkoff+UR)] ) ;
219 return ;
220 }
221
222 /*---------------------------------------------------------------------------
223 Linear interpolation with fixed point arithmetic
224 -----------------------------------------------------------------------------*/
225
extract_byte_liby(int nx,int ny,int nz,byte * vol,Tmask * tm,int fixdir,int fixijk,float da,float db,int ma,int mb,byte * im)226 void extract_byte_liby( int nx , int ny , int nz , byte * vol ,
227 Tmask * tm ,
228 int fixdir , int fixijk , float da , float db ,
229 int ma , int mb , byte * im )
230 {
231 int adel,bdel , abot,atop , bb,bbot,btop , nxy=nx*ny ;
232 register int aa , ijkoff , aoff,boff ;
233 int astep,bstep,cstep , na,nb,nc ;
234 float fa , fb ;
235 float f_a_b , f_ap_b , f_a_bp , f_ap_bp ;
236 byte b_a_b , b_ap_b , b_a_bp , b_ap_bp ;
237 byte * mask ;
238
239 memset( im , 0 , ma*mb ) ; /* initialize output to zero */
240
241 if( fixijk < 0 ) return ;
242
243 ASSIGN_DIRECTIONS ;
244
245 if( fixijk >= nc ) return ;
246
247 adel = (int) da ; if( da < 0.0 ) adel-- ; /* floor(da) */
248 bdel = (int) db ; if( db < 0.0 ) bdel-- ; /* floor(db) */
249
250 fa = da - adel ; /* fractional part of dj */
251 fb = db - bdel ; /* fractional part of dk */
252
253 adel++ ; bdel++ ;
254
255 f_a_b = fa * fb ;
256 f_ap_b = (1.0-fa)* fb ;
257 f_a_bp = fa *(1.0-fb) ;
258 f_ap_bp = (1.0-fa)*(1.0-fb) ;
259
260 bb = (int)(256*f_a_b + 0.499) ; if( bb == 256 ) bb-- ; b_a_b = (byte) bb ;
261 bb = (int)(256*f_ap_b + 0.499) ; if( bb == 256 ) bb-- ; b_ap_b = (byte) bb ;
262 bb = (int)(256*f_a_bp + 0.499) ; if( bb == 256 ) bb-- ; b_a_bp = (byte) bb ;
263 bb = (int)(256*f_ap_bp+ 0.499) ; if( bb == 256 ) bb-- ; b_ap_bp= (byte) bb ;
264
265 abot = 0 ; if( abot < adel ) abot = adel ; /* range in im[] */
266 atop = na+adel-1 ; if( atop > ma ) atop = ma ;
267
268 bbot = 0 ; if( bbot < bdel ) bbot = bdel ;
269 btop = nb+bdel-1 ; if( btop > mb ) btop = mb ;
270
271 ijkoff = fixijk*cstep + (abot-adel)*astep + (bbot-bdel)*bstep ;
272 boff = bbot * ma ;
273
274 mask = (tm == NULL) ? NULL
275 : tm->mask[fixdir%3] + (fixijk*nb - bdel) ;
276
277 for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
278 if( mask == NULL || mask[bb] || mask[bb+1] )
279 for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
280 im[aa+boff] = (byte)(( b_a_b * vol[aoff+ijkoff]
281 + b_ap_b * vol[aoff+(ijkoff+LR)]
282 + b_a_bp * vol[aoff+(ijkoff+UL)]
283 + b_ap_bp * vol[aoff+(ijkoff+UR)] ) >> 8 ) ;
284 return ;
285 }
286
287 /*---------------------------------------------------------------------------
288 Two-step interpolation
289 -----------------------------------------------------------------------------*/
290
291 #if 0
292 # define TSBOT 0.3
293 # define TSTOP 0.7
294 #else
295 # define TSBOT 0.25
296 # define TSTOP 0.75
297 #endif
298
extract_byte_ts(int nx,int ny,int nz,byte * vol,Tmask * tm,int fixdir,int fixijk,float da,float db,int ma,int mb,byte * im)299 void extract_byte_ts( int nx , int ny , int nz , byte * vol ,
300 Tmask * tm ,
301 int fixdir , int fixijk , float da , float db ,
302 int ma , int mb , byte * im )
303 {
304 int adel,bdel , abot,atop , bb,bbot,btop , nxy=nx*ny ;
305 register int aa , ijkoff , aoff,boff ;
306 int astep,bstep,cstep , na,nb,nc , nts,dts1,dts2 ;
307 float fa , fb ;
308 byte * mask ;
309
310 memset( im , 0 , ma*mb ) ; /* initialize output to zero */
311
312 if( fixijk < 0 ) return ;
313
314 ASSIGN_DIRECTIONS ;
315
316 if( fixijk >= nc ) return ;
317
318 adel = (int) da ; if( da < 0.0 ) adel-- ; /* floor(da) */
319 bdel = (int) db ; if( db < 0.0 ) bdel-- ; /* floor(db) */
320
321 fa = da - adel ; /* fractional part of dj */
322 fb = db - bdel ; /* fractional part of dk */
323
324 fa = 1.0-fa ; fb = 1.0-fb ;
325
326 if( fa < TSBOT ){ /*- Left 30% -*/
327 if( fb < TSBOT ){ /*- Lower 30% -*/
328 nts = 1 ; dts1 = LL ; /* [0,0] */
329 } else if( fb > TSTOP ){ /*- Upper 30% -*/
330 nts = 1 ; dts1 = UL ; /* [0,1] */
331 } else { /*- Middle 40% -*/
332 nts = 2 ; dts1 = LL ; dts2 = UL ; /* mid of [0,0] and [0,1] */
333 }
334 } else if( fa > TSTOP ){ /*- Right 30% -*/
335 if( fb < TSBOT ){ /*- Lower 30% -*/
336 nts = 1 ; dts1 = LR ; /* [1,0] */
337 } else if( fb > TSTOP ){ /*- Upper 30% -*/
338 nts = 1 ; dts1 = UR ; /* [1,1] */
339 } else { /*- Middle 40% -*/
340 nts = 2 ; dts1 = LR ; dts2 = UR ; /* mid of [1,0] and [1,1] */
341 }
342 } else { /*- Middle 40% -*/
343 if( fb < TSBOT ){ /*- Lower 30% -*/
344 nts = 2 ; dts1 = LL ; dts2 = LR ; /* mid of [0,0] and [1,0] */
345 } else if( fb > TSTOP ){ /*- Upper 30% -*/
346 nts = 2 ; dts1 = UL ; dts2 = UR ; /* mid of [0,1] and [1,1] */
347 } else { /*- Middle 40% -*/
348 nts = 4 ; /* mid of all 4 points */
349 }
350 }
351
352 adel++ ; bdel++ ;
353
354 abot = 0 ; if( abot < adel ) abot = adel ; /* range in im[] */
355 atop = na+adel-1 ; if( atop > ma ) atop = ma ;
356
357 bbot = 0 ; if( bbot < bdel ) bbot = bdel ;
358 btop = nb+bdel-1 ; if( btop > mb ) btop = mb ;
359
360 ijkoff = fixijk*cstep + (abot-adel)*astep + (bbot-bdel)*bstep ;
361 boff = bbot * ma ;
362
363 mask = (tm == NULL) ? NULL
364 : tm->mask[fixdir%3] + (fixijk*nb - bdel) ;
365
366 switch( nts ){
367
368 case 1:
369 ijkoff += dts1 ;
370 for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
371 if( mask == NULL || mask[bb] || mask[bb+1] )
372 for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
373 im[aa+boff] = vol[aoff+ijkoff] ;
374 break ;
375
376 case 2:
377 ijkoff += dts1 ; dts2 -= dts1 ;
378 for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
379 if( mask == NULL || mask[bb] || mask[bb+1] )
380 for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
381 im[aa+boff] = (vol[aoff+ijkoff] + vol[aoff+(ijkoff+dts2)]) >> 1;
382 break ;
383
384 case 4:
385 for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
386 if( mask == NULL || mask[bb] || mask[bb+1] )
387 for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
388 im[aa+boff] = ( vol[aoff+ijkoff] +vol[aoff+(ijkoff+LR)]
389 +vol[aoff+(ijkoff+UL)]+vol[aoff+(ijkoff+UR)]) >> 2;
390 break ;
391 }
392
393 return ;
394 }
395
396 /*---------------------------------------------------------------------------
397 Four-step interpolation
398 -----------------------------------------------------------------------------*/
399
400 #if 0
401 # define FSA 0.175
402 # define FSB 0.400
403 # define FSC 0.600
404 # define FSD 0.825
405 #else
406 # define FSA 0.125
407 # define FSB 0.375
408 # define FSC 0.625
409 # define FSD 0.875
410 #endif
411
extract_byte_fs(int nx,int ny,int nz,byte * vol,Tmask * tm,int fixdir,int fixijk,float da,float db,int ma,int mb,byte * im)412 void extract_byte_fs( int nx , int ny , int nz , byte * vol ,
413 Tmask * tm ,
414 int fixdir , int fixijk , float da , float db ,
415 int ma , int mb , byte * im )
416 {
417 int adel,bdel , abot,atop , bb,bbot,btop , nxy=nx*ny ;
418 register int aa , ijkoff , aoff,boff ;
419 int astep,bstep,cstep , na,nb,nc , nfs,dfs1,dfs2,dfs3,dfs4 , ap,bp ;
420 float fa , fb ;
421 byte * mask ;
422
423 memset( im , 0 , ma*mb ) ; /* initialize output to zero */
424
425 if( fixijk < 0 ) return ;
426
427 ASSIGN_DIRECTIONS ;
428
429 if( fixijk >= nc ) return ;
430
431 adel = (int) da ; if( da < 0.0 ) adel-- ; /* floor(da) */
432 bdel = (int) db ; if( db < 0.0 ) bdel-- ; /* floor(db) */
433
434 fa = da - adel ; /* fractional part of dj */
435 fb = db - bdel ; /* fractional part of dk */
436
437 fa = 1.0-fa ; fb = 1.0-fb ; /* weights for right/upper sides */
438
439 if( fa < FSA ) ap = 0 ; /* left-right position */
440 else if( fa < FSB ) ap = 1 ;
441 else if( fa < FSC ) ap = 2 ;
442 else if( fa < FSD ) ap = 3 ;
443 else ap = 4 ;
444
445 if( fb < FSA ) bp = 0 ; /* down-up position */
446 else if( fb < FSB ) bp = 1 ;
447 else if( fb < FSC ) bp = 2 ;
448 else if( fb < FSD ) bp = 3 ;
449 else bp = 4 ;
450
451 /*----- 5x5 grid of possible interpolation cases (nfs): -----------------
452
453 bp = 4| 1 3 2 3 1 04 14 24 34 44 <- grid of
454 3| 3 4 5 4 3 03 13 23 33 43 <- 10*ap + bp
455 2| 2 5 6 5 2 02 12 22 32 42 <- values
456 1| 3 4 5 4 3 01 11 21 31 41
457 0| 1 3 2 3 1 00 10 20 30 40
458 -----------
459 ap = 0 1 2 3 4
460
461 ----- The indices and nfs cases are assigned in the switch below. -----*/
462
463
464 dfs2=dfs3=dfs4=-1 ;
465 switch( 10*ap + bp ){
466
467 default: return ; /* should never be executed */
468
469 case 00: nfs = 1 ; dfs1 = LL ; break ; /* 1 point */
470 case 04: nfs = 1 ; dfs1 = UL ; break ;
471 case 40: nfs = 1 ; dfs1 = LR ; break ;
472 case 44: nfs = 1 ; dfs1 = UR ; break ;
473
474 case 20: nfs = 2 ; dfs1 = LL ; dfs2 = LR ; break ; /* 2 points: */
475 case 02: nfs = 2 ; dfs1 = LL ; dfs2 = UL ; break ; /* 1/2 = dfs1 */
476 case 24: nfs = 2 ; dfs1 = UL ; dfs2 = UR ; break ; /* 1/2 = dfs2 */
477 case 42: nfs = 2 ; dfs1 = LR ; dfs2 = UR ; break ;
478
479 case 10: nfs = 3 ; dfs1 = LL ; dfs2 = LR ; break ; /* 2 points: */
480 case 30: nfs = 3 ; dfs1 = LR ; dfs2 = LL ; break ; /* 3/4 = dfs1 */
481 case 01: nfs = 3 ; dfs1 = LL ; dfs2 = UL ; break ; /* 1/4 = dfs2 */
482 case 03: nfs = 3 ; dfs1 = UL ; dfs2 = LL ; break ;
483 case 14: nfs = 3 ; dfs1 = UL ; dfs2 = UR ; break ;
484 case 34: nfs = 3 ; dfs1 = UR ; dfs2 = UL ; break ;
485 case 41: nfs = 3 ; dfs1 = LR ; dfs2 = UR ; break ;
486 case 43: nfs = 3 ; dfs1 = UR ; dfs2 = LR ; break ;
487
488 case 11: nfs = 4 ; dfs1 = LL ; dfs2 = LR ; /* 4 points: */
489 dfs3 = UL ; dfs4 = UR ; break ; /* 9/16 = dfs1 */
490 case 13: nfs = 4 ; dfs1 = UL ; dfs2 = UR ; /* 3/16 = dfs2 */
491 dfs3 = LL ; dfs4 = LR ; break ; /* 3/16 = dfs3 */
492 case 31: nfs = 4 ; dfs1 = LR ; dfs2 = LL ; /* 1/16 = dfs4 */
493 dfs3 = UR ; dfs4 = UL ; break ;
494 case 33: nfs = 4 ; dfs1 = UR ; dfs2 = UL ;
495 dfs3 = LR ; dfs4 = LL ; break ;
496
497 case 12: nfs = 5 ; dfs1 = LL ; dfs2 = UL ; /* 4 points: */
498 dfs3 = LR ; dfs4 = UR ; break ; /* 3/8 = dfs1 */
499 case 21: nfs = 5 ; dfs1 = LL ; dfs2 = LR ; /* 3/8 = dfs2 */
500 dfs3 = UL ; dfs4 = UR ; break ; /* 1/8 = dfs3 */
501 case 23: nfs = 5 ; dfs1 = UL ; dfs2 = UR ; /* 1/8 = dfs4 */
502 dfs3 = LL ; dfs4 = LR ; break ;
503 case 32: nfs = 5 ; dfs1 = LR ; dfs2 = UR ;
504 dfs3 = LL ; dfs4 = UL ; break ;
505
506 case 22: nfs = 6 ; dfs1 = LL ; dfs2 = LR ; /* 4 points: */
507 dfs3 = UL ; dfs4 = UR ; break ; /* 1/4 = all */
508 }
509
510 adel++ ; bdel++ ;
511
512 abot = 0 ; if( abot < adel ) abot = adel ; /* range in im[] */
513 atop = na+adel-1 ; if( atop > ma ) atop = ma ;
514
515 bbot = 0 ; if( bbot < bdel ) bbot = bdel ;
516 btop = nb+bdel-1 ; if( btop > mb ) btop = mb ;
517
518 ijkoff = fixijk*cstep + (abot-adel)*astep + (bbot-bdel)*bstep ;
519 boff = bbot * ma ;
520
521 #if 0
522 printf("fixijk=%3d nfs=%d dfs1=%d dfs2=%d dfs3=%d dfs4=%d\n",
523 fixijk,nfs,dfs1,dfs2,dfs3,dfs4);
524 #endif
525
526 mask = (tm == NULL) ? NULL
527 : tm->mask[fixdir%3] + (fixijk*nb - bdel) ;
528
529 switch( nfs ){
530
531 case 1: /* 1 point (NN copy) */
532 ijkoff += dfs1 ;
533 for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
534 if( mask == NULL || mask[bb] || mask[bb+1] )
535 for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
536 im[aa+boff] = vol[aoff+ijkoff] ;
537 break ;
538
539 case 2: /* 2 points (1/2+1/2) */
540 ijkoff += dfs1 ; dfs2 -= dfs1 ;
541 for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
542 if( mask == NULL || mask[bb] || mask[bb+1] )
543 for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
544 im[aa+boff] = (vol[aoff+ijkoff] + vol[aoff+(ijkoff+dfs2)]) >> 1 ;
545 break ;
546
547 case 3: /* 2 points (3/4+1/4) */
548 ijkoff += dfs1 ; dfs2 -= dfs1 ;
549 for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
550 if( mask == NULL || mask[bb] || mask[bb+1] )
551 for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
552 im[aa+boff] = ( (vol[aoff+ijkoff] << 1) + vol[aoff+ijkoff]
553 + vol[aoff+(ijkoff+dfs2)] ) >> 2 ;
554 break ;
555
556 case 4: /* 4 points (9/16+3/16+3/16+1/16) */
557 ijkoff += dfs1 ; dfs2 -= dfs1 ; dfs3 -= dfs1 ; dfs4 -= dfs1 ;
558 for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
559 if( mask == NULL || mask[bb] || mask[bb+1] )
560 for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
561 im[aa+boff] = ( (vol[aoff+ijkoff] << 3)
562 + vol[aoff+ijkoff]
563 +((vol[aoff+(ijkoff+dfs2)] + vol[aoff+(ijkoff+dfs3)]) << 1)
564 + (vol[aoff+(ijkoff+dfs2)] + vol[aoff+(ijkoff+dfs3)])
565 + vol[aoff+(ijkoff+dfs4)] ) >> 4 ;
566 break ;
567
568 case 5: /* 4 points (3/8+3/8+1/8+1/8) */
569 ijkoff += dfs1 ; dfs2 -= dfs1 ; dfs3 -= dfs1 ; dfs4 -= dfs1 ;
570 for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
571 if( mask == NULL || mask[bb] || mask[bb+1] )
572 for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
573 im[aa+boff] = ( ((vol[aoff+ijkoff] + vol[aoff+(ijkoff+dfs2)]) << 1)
574 + (vol[aoff+ijkoff] + vol[aoff+(ijkoff+dfs2)])
575 + vol[aoff+(ijkoff+dfs3)] + vol[aoff+(ijkoff+dfs4)] ) >> 3 ;
576 break;
577
578 case 6: /* 4 points (1/4+1/4+1/4+1/4) */
579 ijkoff += dfs1 ; dfs2 -= dfs1 ; dfs3 -= dfs1 ; dfs4 -= dfs1 ;
580 for( bb=bbot ; bb < btop ; bb++,boff+=ma,ijkoff+=bstep )
581 if( mask == NULL || mask[bb] || mask[bb+1] )
582 for( aa=abot,aoff=0 ; aa < atop ; aa++,aoff+=astep )
583 im[aa+boff] = ( vol[aoff+ijkoff] + vol[aoff+(ijkoff+dfs2)]
584 + vol[aoff+(ijkoff+dfs3)] + vol[aoff+(ijkoff+dfs4)] ) >> 2 ;
585 break;
586 }
587
588 return ;
589 }
590
591 /*---------------------------------------------------------------------------
592 Test the speeds of the above routines:
593 nrep = number of repetitions to execute
594 ct = float [5] array (must be allocated by caller)
595 ct[0] = CPU time for _nn
596 ct[1] = CPU time for _lifl
597 ct[2] = CPU time for _liby
598 ct[3] = CPU time for _ts
599 ct[4] = CPU time for _fs
600 -----------------------------------------------------------------------------*/
601
extract_byte_speedtest(int nrep,int fixdir,float * ct)602 void extract_byte_speedtest( int nrep , int fixdir , float * ct )
603 {
604 double cputim ;
605 int pp , nx=161,ny=191,nz=141,nxy=nx*ny ,
606 kk , ma,mb,mab , apad,bpad ;
607 float aa=0.347 , bb=-0.521 , da,db ;
608 byte * vin , * vout ;
609 int astep,bstep,cstep , na,nb,nc ;
610
611 ASSIGN_DIRECTIONS ;
612
613 /* setup bricks */
614
615 da = fabs( 0.5*aa*(nc-1.0) ) ; db = fabs( 0.5*bb*(nc-1.0) ) ;
616 apad = (int)(2.0+da) ; bpad = (int)(2.0+db) ;
617 ma = na + 2*apad ; mb = nb + 2*bpad ; mab = ma*mb ;
618
619 vin = (byte *) malloc( sizeof(byte) * (na*nb*nc) ) ;
620 if( vin == NULL ) return ;
621
622 vout = (byte *) malloc( sizeof(byte) * (ma*mb*nc) ) ;
623 if( vout == NULL ){ free(vin) ; return ; }
624
625 vin[0] = 1 ;
626 for( kk=1 ; kk < na*nb*nc ; kk++ ) vin[kk] = (byte)((3*vin[kk-1]+7) % 256) ;
627
628 #undef BTEST
629 #define BTEST(func) do{ cputim = COX_cpu_time() ; \
630 for( pp=0 ; pp < nrep ; pp++ ){ \
631 for( kk=0 ; kk < nc ; kk++ ){ \
632 da = aa*(kk - 0.5*(nc-1.0)) + apad ; \
633 db = bb*(kk - 0.5*(nc-1.0)) + bpad ; \
634 func( nx,ny,nz , vin , \
635 NULL , \
636 fixdir , kk , da , db , \
637 ma , mb , vout + kk*mab ) ; \
638 } \
639 } \
640 cputim = COX_cpu_time() - cputim ; } while(0)
641
642 BTEST(extract_byte_nn) ; ct[0] = cputim ;
643 BTEST(extract_byte_lifl) ; ct[1] = cputim ;
644 BTEST(extract_byte_liby) ; ct[2] = cputim ;
645 BTEST(extract_byte_ts) ; ct[3] = cputim ;
646 BTEST(extract_byte_fs) ; ct[4] = cputim ;
647
648 #undef BTEST
649
650 free(vin) ; free(vout) ; return ;
651 }
652
653 /*-----------------------------------------------------------------------
654 Simple get/put of a fixed plane (no shifting, zero padding).
655 -------------------------------------------------------------------------*/
656
getplane_byte(int nx,int ny,int nz,byte * vol,int fixdir,int fixijk,byte * im)657 void getplane_byte( int nx , int ny , int nz , byte * vol ,
658 int fixdir , int fixijk , byte * im )
659 {
660 int bb , nxy=nx*ny ;
661 register int aa , ijkoff , aoff,boff ;
662 int astep,bstep,cstep , na,nb,nc ;
663
664 if( fixijk < 0 ) return ;
665
666 ASSIGN_DIRECTIONS ;
667
668 if( fixijk >= nc ) return ;
669
670 ijkoff = fixijk*cstep ;
671
672 for( bb=0,boff=0 ; bb < nb ; bb++,boff+=na,ijkoff+=bstep )
673 for( aa=0,aoff=0 ; aa < na ; aa++,aoff+=astep )
674 im[aa+boff] = vol[aoff+ijkoff] ;
675
676 return ;
677 }
678
putplane_byte(int nx,int ny,int nz,byte * vol,int fixdir,int fixijk,byte * im)679 void putplane_byte( int nx , int ny , int nz , byte * vol ,
680 int fixdir , int fixijk , byte * im )
681 {
682 int bb , nxy=nx*ny ;
683 register int aa , ijkoff , aoff,boff ;
684 int astep,bstep,cstep , na,nb,nc ;
685
686 if( fixijk < 0 ) return ;
687
688 ASSIGN_DIRECTIONS ;
689
690 if( fixijk >= nc ) return ;
691
692 ijkoff = fixijk*cstep ;
693
694 for( bb=0,boff=0 ; bb < nb ; bb++,boff+=na,ijkoff+=bstep )
695 for( aa=0,aoff=0 ; aa < na ; aa++,aoff+=astep )
696 vol[aoff+ijkoff] = im[aa+boff] ;
697
698 return ;
699 }
700
701 /******************************************************************************
702 ******************************************************************************
703 ******************************************************************************/
704
705 typedef void gfun( int , int , int , byte * , Tmask * ,
706 int , int , float , float , int , int , byte * ) ;
707
main(int argc,char * argv[])708 int main( int argc , char * argv[] )
709 {
710 THD_3dim_dataset * in_dset ;
711 Tmask * tmask ;
712 int nx,ny,nz,nxy , kk,ii , ma,mb,mab ,
713 apad,bpad , pp,ploop=1,fixdir;
714 float aa , bb , da,db ;
715 THD_ivec3 iv ;
716 byte * vin , * vout , * vmax ;
717 MRI_IMAGE * imout , * immax ;
718 double cputim ;
719 gfun * func = extract_byte_nn ;
720 char * cfun = "nn" ;
721 int astep,bstep,cstep , na,nb,nc , use_tmask ;
722
723 if( argc < 3 ){
724 printf("Usage 1: extor fixdir A B bytedset [loops [suffix]]\n") ;
725 printf("Usage 2: extor fixdir loops\n") ;
726 exit(0) ;
727 }
728
729 fixdir = strtol(argv[1],NULL,10) ;
730 if( fixdir < 0 ){ use_tmask = 1 ; fixdir = -fixdir ; }
731 if( fixdir<1 || fixdir>3 ){fprintf(stderr,"fixdir=%d?\n",fixdir);exit(1);}
732
733 if( argc == 3 ){
734 float ct[5] ;
735 ploop = strtol(argv[2],NULL,10) ;
736 if( ploop < 1 ){ fprintf(stderr,"loop=%d?\n",ploop);exit(1);}
737 extract_byte_speedtest( ploop , fixdir , ct ) ;
738 printf("Speed test with fixdir=%d\n"
739 "_nn = %g (%g/rep)\n"
740 "_lifl = %g (%g/rep)\n"
741 "_liby = %g (%g/rep)\n"
742 "_ts = %g (%g/rep)\n"
743 "_fs = %g (%g/rep)\n" ,
744 fixdir ,
745 ct[0],ct[0]/ploop, ct[1],ct[1]/ploop,
746 ct[2],ct[2]/ploop, ct[3],ct[3]/ploop, ct[4],ct[4]/ploop ) ;
747 exit(1) ;
748 }
749
750 aa = strtod(argv[2],NULL) ;
751 bb = strtod(argv[3],NULL) ;
752 if( aa == 0.0 && bb == 0.0 ){fprintf(stderr,"A=B=0?\n");exit(1);}
753
754 if( argc > 5 ){
755 ploop = strtol(argv[5],NULL,10) ;
756 if( ploop < 1 ){ fprintf(stderr,"loop=%d?\n",ploop);exit(1); }
757 }
758
759 if( argc > 6 ){
760 cfun = argv[6] ;
761 if( strstr(argv[6],"nn") != NULL )
762 func = extract_byte_nn ;
763 else if( strstr(argv[6],"lifl") != NULL )
764 func = extract_byte_lifl ;
765 else if( strstr(argv[6],"liby") != NULL )
766 func = extract_byte_liby ;
767 else if( strstr(argv[6],"ts") != NULL )
768 func = extract_byte_ts ;
769 else if( strstr(argv[6],"fs") != NULL )
770 func = extract_byte_fs ;
771 else {
772 fprintf(stderr,"Unknown func suffix\n");exit(1);
773 }
774 }
775
776 in_dset = THD_open_dataset( argv[4] ) ;
777 if( in_dset == NULL ){fprintf(stderr,"can't open dataset?\n");exit(1);}
778 if( DSET_NVALS(in_dset) > 1 ){fprintf(stderr,"nvals > 1?\n");exit(1);}
779 if( DSET_BRICK_TYPE(in_dset,0) != MRI_byte ){fprintf(stderr,"not byte?\n");exit(1);}
780
781 nx = DSET_NX(in_dset) ;
782 ny = DSET_NY(in_dset) ;
783 nz = DSET_NZ(in_dset) ; nxy = nx*ny ;
784
785 ASSIGN_DIRECTIONS ;
786
787 da = fabs( 0.5*aa*(nc-1.0) ) ; db = fabs( 0.5*bb*(nc-1.0) ) ;
788 if( da < 1.0 && db < 1.0 ){fprintf(stderr,"da=%g db=%g ?\n",da,db);exit(1);}
789
790 apad = (int)(2.0+da) ; bpad = (int)(2.0+db) ;
791 ma = na + 2*apad ; mb = nb + 2*bpad ; mab = ma*mb ;
792
793 DSET_load(in_dset) ;
794 vin = DSET_BRICK_ARRAY(in_dset,0) ;
795
796 imout = mri_new( ma,mb , MRI_byte ) ; vout = MRI_BYTE_PTR(imout) ;
797 immax = mri_new( ma,mb , MRI_byte ) ; vmax = MRI_BYTE_PTR(immax) ;
798
799 tmask = (use_tmask) ? create_Tmask(nx,ny,nz,vin) : NULL ;
800
801 cputim = COX_cpu_time() ;
802
803 for( pp=0 ; pp < ploop ; pp++ ){
804 memset( vmax , 0 , mab ) ;
805 for( kk=0 ; kk < nc ; kk++ ){
806 da = aa*(kk - 0.5*(nc-1.0)) + apad ;
807 db = bb*(kk - 0.5*(nc-1.0)) + bpad ;
808
809 func( nx,ny,nz , vin,tmask , fixdir,kk , da,db , ma,mb , vout ) ;
810
811 for( ii=0 ; ii < mab ; ii++ )
812 if( vout[ii] > vmax[ii] ) vmax[ii] = vout[ii] ;
813 }
814 }
815
816 cputim = (COX_cpu_time() - cputim)/ploop ;
817 fprintf(stderr,"CPU time per loop = %g [%s]\n",cputim,cfun) ;
818
819 { char fname[128] = "exim_" ;
820 strcat(fname,cfun) ; strcat(fname,".pgm") ;
821 mri_write( fname , immax ) ;
822 }
823
824 exit(0) ;
825 }
826