1 #include "mrilib.h"
2
3 /*------------------------------------------------------------------- Create a new dataset from an old one, with zero padding around
4 the edges. For example,
5 add_I = number of zero planes to add at inferior edge
6 (if < 0, number of data planes to cut off inferior edge)
7
8 09 Feb 2001 - added "flag" input, which is the OR (|) of
9 ZPAD_EMPTY = produce only an padded "empty copy" of the input
10 ZPAD_PURGE = purge input dataset bricks after they are copied
11 ZPAD_MM = increments are mm instead of slice counts
12 (at least 'add_?' mm will be added/subtracted)
13 ZPAD_IJK = increments are relative to dset axes I0--I1 J0--J1
14 K0--K1 not I--S A--P L---R
15 14 May 2002: if inputs crops are all zero, return something anyway
16 ---------------------------------------------------------------------*/
17
THD_zeropad(THD_3dim_dataset * inset,int add_I,int add_S,int add_A,int add_P,int add_L,int add_R,char * prefix,int flag)18 THD_3dim_dataset * THD_zeropad( THD_3dim_dataset *inset ,
19 int add_I , int add_S , int add_A ,
20 int add_P , int add_L , int add_R ,
21 char *prefix , int flag )
22 {
23 THD_3dim_dataset *outset ;
24 int nxold,nyold,nzold , nxnew,nynew,nznew , nxyold,nxynew ,
25 nxbot=0,nxtop=0 , nybot=0,nytop=0 , nzbot=0,nztop=0 ;
26 int ii,jj,kk , iv , iibot,iitop , jjbot,jjtop , kkbot,kktop;
27
28 int empty_flag = (flag & ZPAD_EMPTY) ; /* 09 Feb 2001 */
29 int purge_flag = (flag & ZPAD_PURGE) ; /* 09 Feb 2001 */
30 int mm_flag = (flag & ZPAD_MM ) ; /* 13 Feb 2001 */
31 int ijk_flag = (flag & ZPAD_IJK ); /* ZSS: 23 Dec The year of the war on Christmas */
32
33 THD_ivec3 iv_nxyz ;
34 THD_fvec3 fv_xyzorg ;
35
36 MRI_IMAGE *oldim ;
37 void *vnew ;
38
39 ENTRY("THD_zeropad") ;
40
41 /*-- check inputs --*/
42
43 if( !ISVALID_DSET(inset) ) RETURN( NULL ) ;
44
45 if( !THD_filename_ok(prefix) ) prefix = "zeropad" ;
46
47 if( add_I==0 && add_S==0 && add_P==0 &&
48 add_A==0 && add_L==0 && add_R==0 ){
49
50 INFO_message("THD_zeropad: all pad values are zero - just copying dataset") ;
51 if( !empty_flag ){
52 outset = EDIT_full_copy( inset , prefix ) ; /* 14 May 2002 */
53 } else {
54 outset = EDIT_empty_copy( inset ) ; /* 15 Sep 2020 - oops */
55 EDIT_dset_items( outset , ADN_prefix , prefix , ADN_none ) ;
56 EDIT_dset_items( outset , ADN_prefix , prefix , ADN_none ) ;
57 }
58 RETURN( outset );
59 }
60
61 /*-- map add_? values into dataset xyz coordinate directions --*/
62
63 nxold = DSET_NX(inset) ;
64 nyold = DSET_NY(inset) ;
65 nzold = DSET_NZ(inset) ;
66
67 if (ijk_flag) { /* ZSS Dec 23 05 */
68 nxbot = add_I; nxtop = add_S;
69 nybot = add_A; nytop = add_P;
70 nzbot = add_L; nztop = add_R;
71 } else {
72 /* comput n?top and n?bot, the number of planes to add at
73 the top and bottom of the ? direction, for ? = x, y, or z */
74
75 switch( inset->daxes->xxorient ){
76 default:
77 fprintf(stderr,"*** THD_zeropad: Unknown orientation codes!\n") ;
78 RETURN( NULL );
79
80 case ORI_R2L_TYPE: nxtop = add_L ; nxbot = add_R ; break ;
81 case ORI_L2R_TYPE: nxtop = add_R ; nxbot = add_L ; break ;
82 case ORI_P2A_TYPE: nxtop = add_A ; nxbot = add_P ; break ;
83 case ORI_A2P_TYPE: nxtop = add_P ; nxbot = add_A ; break ;
84 case ORI_I2S_TYPE: nxtop = add_S ; nxbot = add_I ; break ;
85 case ORI_S2I_TYPE: nxtop = add_I ; nxbot = add_S ; break ;
86 }
87
88 switch( inset->daxes->yyorient ){
89 default:
90 fprintf(stderr,"*** THD_zeropad: Unknown orientation codes!\n") ;
91 RETURN( NULL );
92
93 case ORI_R2L_TYPE: nytop = add_L ; nybot = add_R ; break ;
94 case ORI_L2R_TYPE: nytop = add_R ; nybot = add_L ; break ;
95 case ORI_P2A_TYPE: nytop = add_A ; nybot = add_P ; break ;
96 case ORI_A2P_TYPE: nytop = add_P ; nybot = add_A ; break ;
97 case ORI_I2S_TYPE: nytop = add_S ; nybot = add_I ; break ;
98 case ORI_S2I_TYPE: nytop = add_I ; nybot = add_S ; break ;
99 }
100
101 switch( inset->daxes->zzorient ){
102 default:
103 fprintf(stderr,"*** THD_zeropad: Unknown orientation codes!\n") ;
104 RETURN( NULL );
105
106 case ORI_R2L_TYPE: nztop = add_L ; nzbot = add_R ; break ;
107 case ORI_L2R_TYPE: nztop = add_R ; nzbot = add_L ; break ;
108 case ORI_P2A_TYPE: nztop = add_A ; nzbot = add_P ; break ;
109 case ORI_A2P_TYPE: nztop = add_P ; nzbot = add_A ; break ;
110 case ORI_I2S_TYPE: nztop = add_S ; nzbot = add_I ; break ;
111 case ORI_S2I_TYPE: nztop = add_I ; nzbot = add_S ; break ;
112 }
113
114 /* 13 Feb 2001: round to millimeters? */
115
116 #undef RMM
117 #define RMM(n,d) \
118 do{ if( (n) > 0 ) (n) = (int)( (n)/fabs(d) + 0.999 ) ; \
119 else if( (n) < 0 ) (n) = (int)( (n)/fabs(d) - 0.999 ) ; } while(0) ;
120
121 if( mm_flag ){
122 RMM(nxtop,inset->daxes->xxdel) ; RMM(nxbot,inset->daxes->xxdel) ;
123 RMM(nytop,inset->daxes->yydel) ; RMM(nybot,inset->daxes->yydel) ;
124 RMM(nztop,inset->daxes->zzdel) ; RMM(nzbot,inset->daxes->zzdel) ;
125 }
126 }
127
128 nxnew = nxold + nxbot + nxtop ; /* dimensions of new bricks */
129 nynew = nyold + nybot + nytop ;
130 nznew = nzold + nzbot + nztop ;
131
132 nxyold = nxold * nyold ; /* for computing subscripts */
133 nxynew = nxnew * nynew ;
134
135 iibot = MAX(0,-nxbot) ; iitop = MIN(nxold,nxold+nxtop) ; /* range of data */
136 jjbot = MAX(0,-nybot) ; jjtop = MIN(nyold,nyold+nytop) ; /* in old dataset */
137 kkbot = MAX(0,-nzbot) ; kktop = MIN(nzold,nzold+nztop) ;
138
139 if( nxnew < 1 || iibot > iitop || /* check for reasonable sizes */
140 nynew < 1 || jjbot > jjtop || /* and ranges of dataset */
141 nznew < 1 || kkbot > kktop ){
142
143 ERROR_message("*** THD_zeropad: Can't cut dataset down too much!\n"
144 " new: %d %d %d\n"
145 " ii: %d..%d\n"
146 " jj: %d..%d\n"
147 " kk: %d..%d\n",
148 nxnew, nynew, nznew,
149 iibot, iitop, jjbot, jjtop, kkbot, kktop ) ;
150 RETURN( NULL );
151 }
152
153 #if 0
154 if( nxnew < 2 || iibot >= iitop || /* check for reasonable sizes */
155 nynew < 2 || jjbot >= jjtop || /* and ranges of dataset */
156 nznew < 2 || kkbot >= kktop ){
157
158 fprintf(stderr,"*** WARNING - THD_zeropad: dataset cut down to %dx%dx%d\n",
159 nxnew,nynew,nznew) ;
160 }
161 #endif
162
163 /*-- create the shell of the new dataset --*/
164
165 outset = EDIT_empty_copy( inset ) ;
166
167 LOAD_IVEC3( iv_nxyz , nxnew,nynew,nznew ) ;
168
169 LOAD_FVEC3( fv_xyzorg, inset->daxes->xxorg - nxbot * inset->daxes->xxdel,
170 inset->daxes->yyorg - nybot * inset->daxes->yydel,
171 inset->daxes->zzorg - nzbot * inset->daxes->zzdel );
172
173 #if 0
174 DUMP_IVEC3("THD_zeropad new dimensions",iv_nxyz) ;
175 DUMP_FVEC3(" new origin ",fv_xyzorg) ;
176 INFO_message(" grid spacings = %g %g %g",
177 inset->daxes->xxdel,inset->daxes->yydel,inset->daxes->zzdel) ;
178 #endif
179
180 STATUS("setting new dimensions") ;
181
182 EDIT_dset_items( outset ,
183 ADN_prefix , prefix ,
184 ADN_nxyz , iv_nxyz ,
185 ADN_xyzorg , fv_xyzorg ,
186 ADN_none ) ;
187
188 /* Changing dimensions means old anat parent is no longer valid! */
189
190 EDIT_ZERO_ANATOMY_PARENT_ID( outset ) ;
191 outset->anat_parent_name[0] = '\0' ;
192
193 #if 0
194 /* if changing number of slices, can't keep slice-dependent time shifts! */
195
196 if( (nzbot!=0 || nztop!=0) && outset->taxis != NULL && outset->taxis->nsl > 0 ){
197 EDIT_dset_items( outset , ADN_nsl , 0 , ADN_none ) ;
198 fprintf(stderr,
199 "*** THD_zeropad: warning - slice-dependent time shifts have been removed!\n") ;
200 }
201 #else
202 /* 31 Jan 2001: OK, lets keeps the slice-dependent time shifts
203 (but we'll have to mangle them) -- RWCox */
204
205 if( (nzbot!=0 || nztop!=0) && outset->taxis != NULL && outset->taxis->nsl > 0 ){
206 int old_nsl , new_nsl , ii, nkeep,kbot,ibot ;
207 float old_zorg_sl , *old_toff_sl , new_zorg_sl , *new_toff_sl ;
208
209 /* copy current conditions to local variables */
210
211 old_nsl = outset->taxis->nsl ;
212 old_zorg_sl = outset->taxis->zorg_sl ;
213 old_toff_sl = (float *) malloc(sizeof(float)*old_nsl) ;
214 memcpy( old_toff_sl , outset->taxis->toff_sl , sizeof(float)*old_nsl ) ;
215
216 /* compute new values */
217
218 new_nsl = nznew ;
219 new_zorg_sl = outset->daxes->zzorg ; /* cf. to3d.c */
220 new_toff_sl = (float *) malloc(sizeof(float)*new_nsl) ;
221 for( ii=0 ; ii < new_nsl ; ii++ ) new_toff_sl[ii] = 0.0 ; /* extras are 0 */
222
223 nkeep = old_nsl ; /* how many to keep from the old list */
224 if( nzbot < 0 ) nkeep += nzbot ; /* lost this many at the bottom */
225 if( nztop < 0 ) nkeep += nztop ; /* lost this many at the top */
226
227 if( nzbot < 0 ){
228 kbot = -nzbot ; /* which old one to start with */
229 ibot = 0 ; /* and where it goes in new list */
230 } else {
231 kbot = 0 ;
232 ibot = nzbot ;
233 }
234
235 memcpy( new_toff_sl+ibot , old_toff_sl+kbot , sizeof(float)*nkeep ) ;
236
237 /* set new values in dataset */
238
239 STATUS("setting new time-offsets") ;
240
241 EDIT_dset_items( outset ,
242 ADN_nsl , new_nsl ,
243 ADN_toff_sl , new_toff_sl ,
244 ADN_zorg_sl , new_zorg_sl ,
245 ADN_none ) ;
246
247 free(new_toff_sl) ; free(old_toff_sl) ;
248 }
249 #endif
250
251 if( empty_flag ) RETURN(outset) ; /* 09 Feb 2001 */
252
253 /*-- now read the old dataset in, and make bricks for the new dataset --*/
254
255 STATUS("reading dataset in") ;
256
257 DSET_load(inset) ;
258 if( !DSET_LOADED(inset) ){
259 fprintf(stderr,"*** THD_zeropad: Can't load input dataset BRIK!\n");
260 DSET_delete(outset) ;
261 RETURN( NULL );
262 }
263
264 STATUS("padding") ;
265
266 for( iv=0 ; iv < DSET_NVALS(inset) ; iv++ ){
267
268 /* create a brick of zeros */
269
270 oldim = DSET_BRICK(inset,iv) ; /* image structure of old brick */
271
272 vnew = (void*)calloc( nxnew*nynew*nznew , oldim->pixel_size ) ; /* new brick */
273 if( vnew == NULL ){
274 fprintf(stderr,
275 "*** THD_zeropad: Can't malloc space for new sub-brick %d\n",
276 iv) ;
277 DSET_delete(outset) ; RETURN( NULL );
278 }
279
280 /* macros for computing 1D subscripts from 3D indices */
281
282 #undef SNEW /* in case was defined in some stupid .h file */
283 #undef SOLD
284 #define SNEW(i,j,k) ((i+nxbot)+(j+nybot)*nxnew+(k+nzbot)*nxynew)
285 #define SOLD(i,j,k) (i+j*nxold+k*nxyold)
286
287 switch( oldim->kind ){ /* copy rows of old into new */
288
289 default: break ;
290
291 case MRI_byte:{
292 byte * bnew = (byte *) vnew, * bold = mri_data_pointer(oldim) ;
293 for( kk=kkbot ; kk < kktop ; kk++ )
294 for( jj=jjbot ; jj < jjtop ; jj++ )
295 for( ii=iibot ; ii < iitop ; ii++ )
296 bnew[SNEW(ii,jj,kk)] = bold[SOLD(ii,jj,kk)] ;
297 }
298 break ;
299
300 case MRI_short:{
301 short * bnew = (short *) vnew, * bold = mri_data_pointer(oldim) ;
302 for( kk=kkbot ; kk < kktop ; kk++ )
303 for( jj=jjbot ; jj < jjtop ; jj++ )
304 for( ii=iibot ; ii < iitop ; ii++ )
305 bnew[SNEW(ii,jj,kk)] = bold[SOLD(ii,jj,kk)] ;
306 }
307 break ;
308
309 case MRI_float:{
310 float * bnew = (float *) vnew, * bold = mri_data_pointer(oldim) ;
311 for( kk=kkbot ; kk < kktop ; kk++ )
312 for( jj=jjbot ; jj < jjtop ; jj++ )
313 for( ii=iibot ; ii < iitop ; ii++ )
314 bnew[SNEW(ii,jj,kk)] = bold[SOLD(ii,jj,kk)] ;
315 }
316 break ;
317
318 case MRI_complex:{
319 complex * bnew = (complex *) vnew, * bold = mri_data_pointer(oldim) ;
320 for( kk=kkbot ; kk < kktop ; kk++ )
321 for( jj=jjbot ; jj < jjtop ; jj++ )
322 for( ii=iibot ; ii < iitop ; ii++ )
323 bnew[SNEW(ii,jj,kk)] = bold[SOLD(ii,jj,kk)] ;
324 }
325 break ;
326
327 case MRI_int:{
328 int * bnew = (int *) vnew, * bold = mri_data_pointer(oldim) ;
329 for( kk=kkbot ; kk < kktop ; kk++ )
330 for( jj=jjbot ; jj < jjtop ; jj++ )
331 for( ii=iibot ; ii < iitop ; ii++ )
332 bnew[SNEW(ii,jj,kk)] = bold[SOLD(ii,jj,kk)] ;
333 }
334 break ;
335
336 case MRI_rgb:{
337 rgbyte * bnew = (rgbyte *) vnew, * bold = mri_data_pointer(oldim) ;
338 for( kk=kkbot ; kk < kktop ; kk++ )
339 for( jj=jjbot ; jj < jjtop ; jj++ )
340 for( ii=iibot ; ii < iitop ; ii++ )
341 bnew[SNEW(ii,jj,kk)] = bold[SOLD(ii,jj,kk)] ;
342 }
343 break ;
344
345 } /* end of switch on sub-brick type */
346
347 if( purge_flag) DSET_unload_one(inset,iv) ; /* 09 Feb 2001 */
348
349 EDIT_substitute_brick( outset , iv , oldim->kind , vnew ) ;
350
351 } /* end of loop on sub-brick index */
352
353 #if 0
354 if( purge_flag ) DSET_unload(inset) ; /* 09 Feb 2001 */
355 #endif
356
357 RETURN( outset );
358 }
359
360 /*-------------------------------------------------------------------*/
361 /* Create a 1-volume dataset built from
362 (a) a master dataset 'mset' -- for the geometry
363 (b) an image 'imin' containing the data
364 (c) padding amounts that reckon with the dimensional
365 differences between imin and mset -- that is, for positive
366 'pad' values, imin should be bigger than mset.
367 This function is for use in 3dAllineate. [25 Jan 2021 - RWC]
368 *//*-----------------------------------------------------------------*/
369
THD_volume_to_dataset(THD_3dim_dataset * mset,MRI_IMAGE * imin,char * prefix,int pad_xm,int pad_xp,int pad_ym,int pad_yp,int pad_zm,int pad_zp)370 THD_3dim_dataset * THD_volume_to_dataset( THD_3dim_dataset *mset ,
371 MRI_IMAGE *imin ,
372 char *prefix ,
373 int pad_xm , int pad_xp ,
374 int pad_ym , int pad_yp ,
375 int pad_zm , int pad_zp )
376 {
377 THD_3dim_dataset *dset ;
378 MRI_IMAGE *dim ;
379 int nxx,nyy,nzz ;
380
381 ENTRY("THD_volume_to_dataset") ;
382
383 if( mset == NULL || imin == NULL ) RETURN(NULL) ;
384
385 if( !THD_filename_ok(prefix) ) prefix = "volumized" ;
386
387 /* size of dataset after it gets padded */
388
389 nxx = DSET_NX(mset) + pad_xm + pad_xp ;
390 nyy = DSET_NY(mset) + pad_ym + pad_yp ;
391 nzz = DSET_NZ(mset) + pad_zm + pad_zp ;
392
393 /* should match input image */
394
395 if( nxx != imin->nx || nyy != imin->ny || nzz != imin->nz ){
396 WARNING_message("THD_volume_to_dataset mismatch:\n"
397 " Dataset = %3d %3d %3d\n"
398 " x-pad = %3d %3d\n"
399 " y-pad = %3d %3d\n"
400 " z-pad = %3d %3d\n"
401 " Image = %3d %3d %3d" ,
402 DSET_NX(mset),DSET_NY(mset),DSET_NZ(mset) ,
403 pad_xm , pad_xp ,
404 pad_ym , pad_yp ,
405 pad_zm , pad_zp ,
406 imin->nx , imin->ny , imin->nz ) ;
407 RETURN(NULL) ;
408 }
409
410 /* create an empty (no data) dataset properly padded from mset */
411
412 dset = THD_zeropad( mset ,
413 pad_xm , pad_xp ,
414 pad_ym , pad_yp ,
415 pad_zm , pad_zp ,
416 prefix , ZPAD_EMPTY | ZPAD_IJK ) ;
417
418 EDIT_dset_items( dset ,
419 ADN_nvals , 1 ,
420 ADN_ntt , 0 ,
421 ADN_datum_all , imin->kind ,
422 ADN_none ) ;
423
424 /* copy the data */
425
426 dim = mri_copy(imin) ;
427
428 /* shove it into the new dataset */
429
430 EDIT_substitute_brick( dset , 0 , dim->kind , mri_data_pointer(dim) ) ;
431
432 /* erase the shell of the image copy */
433
434 mri_clear_and_free( dim ) ;
435
436 /* get the hell out of Dodge */
437
438 RETURN(dset) ;
439 }
440
441 /*-------------------------------------------------------------------*/
442 /* Create a dataset from an MRI_IMARR.
443 All sub-images must be the same dimensions.
444 All sub-images will be converted (if necessary) to the same
445 datum type as the [0] image.
446 NULL is returned if inputs are stooopid.
447 [RWC - 05 Mar 2021]
448 *//*-----------------------------------------------------------------*/
449
THD_imarr_to_dataset(MRI_IMARR * imar,char * prefix)450 THD_3dim_dataset * THD_imarr_to_dataset( MRI_IMARR *imar, char *prefix )
451 {
452 int nx,ny,nz , nvals , vv ;
453 MRI_IMAGE *qim ;
454 MRI_TYPE qtyp ;
455 THD_3dim_dataset *qset ;
456 THD_ivec3 nxyz ;
457 THD_fvec3 dxyz ;
458
459 ENTRY("THD_imarr_to_dataset") ;
460
461 if( imar == NULL || IMARR_COUNT(imar) == 0 ) RETURN(NULL) ;
462
463 qim = IMARR_SUBIM(imar,0); if( qim == NULL ) RETURN(NULL) ;
464
465 /* dimensions */
466
467 nx = qim->nx ; ny = qim->ny ; nz = qim->nz ;
468 if( nx < 2 && ny < 2 && nz < 2 ) RETURN(NULL) ;
469 LOAD_IVEC3(nxyz,nx,ny,nz) ;
470 LOAD_FVEC3(dxyz,1.0f,1.0f,1.0f) ;
471
472 qtyp = qim->kind ;
473 nvals = IMARR_COUNT(imar) ;
474
475 /* dimension check */
476
477 for( vv=1 ; vv < nvals ; vv++ ){
478 qim = IMARR_SUBIM(imar,vv) ; if( qim == NULL ) RETURN(NULL) ;
479 if( qim->nx != nx || qim->ny != ny || qim->nz != nz ) RETURN(NULL) ;
480 }
481
482 /* create empty shell of dataset and set dimensions */
483
484 qset = EDIT_empty_copy(NULL) ;
485
486 if( !THD_filename_ok(prefix) ) prefix = "fromIMARR" ;
487
488 EDIT_dset_items( qset ,
489 ADN_datum_all , (int)qtyp ,
490 ADN_nvals , nvals ,
491 ADN_nxyz , nxyz ,
492 ADN_xyzdel , dxyz ,
493 ADN_prefix , prefix ,
494 ADN_none ) ;
495
496 /* make a copy of each image,
497 and stuff its data into the dataset,
498 then clear out copy's data (now inside dataset) and free image struct */
499
500 for( vv=0 ; vv < nvals ; vv++ ){
501 qim = mri_to_mri( (int)qtyp , IMARR_SUBIM(imar,vv) ) ;
502 EDIT_substitute_brick( qset , vv , (int)qtyp , mri_data_pointer(qim) ) ;
503 mri_clear_and_free( qim ) ;
504 }
505
506 RETURN(qset) ;
507 }
508
509 /*-------------------------------------------------------------------*/
510
THD_image_to_dataset(MRI_IMAGE * imin,char * prefix)511 THD_3dim_dataset * THD_image_to_dataset( MRI_IMAGE *imin, char *prefix )
512 {
513 MRI_IMARR *imar ;
514 THD_3dim_dataset *qset ;
515
516 ENTRY("THD_image_to_dataset") ;
517
518 if( imin == NULL ) RETURN(NULL) ;
519
520 INIT_IMARR(imar) ; ADDTO_IMARR(imar,imin) ;
521
522 qset = THD_imarr_to_dataset( imar , prefix ) ;
523
524 FREE_IMARR(imar) ;
525
526 RETURN(qset) ;
527 }
528