#include "mrilib.h" /*------------------------------------------------------------------- Create a new dataset from an old one, with zero padding around the edges. For example, add_I = number of zero planes to add at inferior edge (if < 0, number of data planes to cut off inferior edge) 09 Feb 2001 - added "flag" input, which is the OR (|) of ZPAD_EMPTY = produce only an padded "empty copy" of the input ZPAD_PURGE = purge input dataset bricks after they are copied ZPAD_MM = increments are mm instead of slice counts (at least 'add_?' mm will be added/subtracted) ZPAD_IJK = increments are relative to dset axes I0--I1 J0--J1 K0--K1 not I--S A--P L---R 14 May 2002: if inputs crops are all zero, return something anyway ---------------------------------------------------------------------*/ THD_3dim_dataset * 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 ) { THD_3dim_dataset *outset ; int nxold,nyold,nzold , nxnew,nynew,nznew , nxyold,nxynew , nxbot=0,nxtop=0 , nybot=0,nytop=0 , nzbot=0,nztop=0 ; int ii,jj,kk , iv , iibot,iitop , jjbot,jjtop , kkbot,kktop; int empty_flag = (flag & ZPAD_EMPTY) ; /* 09 Feb 2001 */ int purge_flag = (flag & ZPAD_PURGE) ; /* 09 Feb 2001 */ int mm_flag = (flag & ZPAD_MM ) ; /* 13 Feb 2001 */ int ijk_flag = (flag & ZPAD_IJK ); /* ZSS: 23 Dec The year of the war on Christmas */ THD_ivec3 iv_nxyz ; THD_fvec3 fv_xyzorg ; MRI_IMAGE *oldim ; void *vnew ; ENTRY("THD_zeropad") ; /*-- check inputs --*/ if( !ISVALID_DSET(inset) ) RETURN( NULL ) ; if( !THD_filename_ok(prefix) ) prefix = "zeropad" ; if( add_I==0 && add_S==0 && add_P==0 && add_A==0 && add_L==0 && add_R==0 ){ INFO_message("THD_zeropad: all pad values are zero - just copying dataset") ; if( !empty_flag ){ outset = EDIT_full_copy( inset , prefix ) ; /* 14 May 2002 */ } else { outset = EDIT_empty_copy( inset ) ; /* 15 Sep 2020 - oops */ EDIT_dset_items( outset , ADN_prefix , prefix , ADN_none ) ; EDIT_dset_items( outset , ADN_prefix , prefix , ADN_none ) ; } RETURN( outset ); } /*-- map add_? values into dataset xyz coordinate directions --*/ nxold = DSET_NX(inset) ; nyold = DSET_NY(inset) ; nzold = DSET_NZ(inset) ; if (ijk_flag) { /* ZSS Dec 23 05 */ nxbot = add_I; nxtop = add_S; nybot = add_A; nytop = add_P; nzbot = add_L; nztop = add_R; } else { /* comput n?top and n?bot, the number of planes to add at the top and bottom of the ? direction, for ? = x, y, or z */ switch( inset->daxes->xxorient ){ default: fprintf(stderr,"*** THD_zeropad: Unknown orientation codes!\n") ; RETURN( NULL ); case ORI_R2L_TYPE: nxtop = add_L ; nxbot = add_R ; break ; case ORI_L2R_TYPE: nxtop = add_R ; nxbot = add_L ; break ; case ORI_P2A_TYPE: nxtop = add_A ; nxbot = add_P ; break ; case ORI_A2P_TYPE: nxtop = add_P ; nxbot = add_A ; break ; case ORI_I2S_TYPE: nxtop = add_S ; nxbot = add_I ; break ; case ORI_S2I_TYPE: nxtop = add_I ; nxbot = add_S ; break ; } switch( inset->daxes->yyorient ){ default: fprintf(stderr,"*** THD_zeropad: Unknown orientation codes!\n") ; RETURN( NULL ); case ORI_R2L_TYPE: nytop = add_L ; nybot = add_R ; break ; case ORI_L2R_TYPE: nytop = add_R ; nybot = add_L ; break ; case ORI_P2A_TYPE: nytop = add_A ; nybot = add_P ; break ; case ORI_A2P_TYPE: nytop = add_P ; nybot = add_A ; break ; case ORI_I2S_TYPE: nytop = add_S ; nybot = add_I ; break ; case ORI_S2I_TYPE: nytop = add_I ; nybot = add_S ; break ; } switch( inset->daxes->zzorient ){ default: fprintf(stderr,"*** THD_zeropad: Unknown orientation codes!\n") ; RETURN( NULL ); case ORI_R2L_TYPE: nztop = add_L ; nzbot = add_R ; break ; case ORI_L2R_TYPE: nztop = add_R ; nzbot = add_L ; break ; case ORI_P2A_TYPE: nztop = add_A ; nzbot = add_P ; break ; case ORI_A2P_TYPE: nztop = add_P ; nzbot = add_A ; break ; case ORI_I2S_TYPE: nztop = add_S ; nzbot = add_I ; break ; case ORI_S2I_TYPE: nztop = add_I ; nzbot = add_S ; break ; } /* 13 Feb 2001: round to millimeters? */ #undef RMM #define RMM(n,d) \ do{ if( (n) > 0 ) (n) = (int)( (n)/fabs(d) + 0.999 ) ; \ else if( (n) < 0 ) (n) = (int)( (n)/fabs(d) - 0.999 ) ; } while(0) ; if( mm_flag ){ RMM(nxtop,inset->daxes->xxdel) ; RMM(nxbot,inset->daxes->xxdel) ; RMM(nytop,inset->daxes->yydel) ; RMM(nybot,inset->daxes->yydel) ; RMM(nztop,inset->daxes->zzdel) ; RMM(nzbot,inset->daxes->zzdel) ; } } nxnew = nxold + nxbot + nxtop ; /* dimensions of new bricks */ nynew = nyold + nybot + nytop ; nznew = nzold + nzbot + nztop ; nxyold = nxold * nyold ; /* for computing subscripts */ nxynew = nxnew * nynew ; iibot = MAX(0,-nxbot) ; iitop = MIN(nxold,nxold+nxtop) ; /* range of data */ jjbot = MAX(0,-nybot) ; jjtop = MIN(nyold,nyold+nytop) ; /* in old dataset */ kkbot = MAX(0,-nzbot) ; kktop = MIN(nzold,nzold+nztop) ; if( nxnew < 1 || iibot > iitop || /* check for reasonable sizes */ nynew < 1 || jjbot > jjtop || /* and ranges of dataset */ nznew < 1 || kkbot > kktop ){ ERROR_message("*** THD_zeropad: Can't cut dataset down too much!\n" " new: %d %d %d\n" " ii: %d..%d\n" " jj: %d..%d\n" " kk: %d..%d\n", nxnew, nynew, nznew, iibot, iitop, jjbot, jjtop, kkbot, kktop ) ; RETURN( NULL ); } #if 0 if( nxnew < 2 || iibot >= iitop || /* check for reasonable sizes */ nynew < 2 || jjbot >= jjtop || /* and ranges of dataset */ nznew < 2 || kkbot >= kktop ){ fprintf(stderr,"*** WARNING - THD_zeropad: dataset cut down to %dx%dx%d\n", nxnew,nynew,nznew) ; } #endif /*-- create the shell of the new dataset --*/ outset = EDIT_empty_copy( inset ) ; LOAD_IVEC3( iv_nxyz , nxnew,nynew,nznew ) ; LOAD_FVEC3( fv_xyzorg, inset->daxes->xxorg - nxbot * inset->daxes->xxdel, inset->daxes->yyorg - nybot * inset->daxes->yydel, inset->daxes->zzorg - nzbot * inset->daxes->zzdel ); #if 0 DUMP_IVEC3("THD_zeropad new dimensions",iv_nxyz) ; DUMP_FVEC3(" new origin ",fv_xyzorg) ; INFO_message(" grid spacings = %g %g %g", inset->daxes->xxdel,inset->daxes->yydel,inset->daxes->zzdel) ; #endif STATUS("setting new dimensions") ; EDIT_dset_items( outset , ADN_prefix , prefix , ADN_nxyz , iv_nxyz , ADN_xyzorg , fv_xyzorg , ADN_none ) ; /* Changing dimensions means old anat parent is no longer valid! */ EDIT_ZERO_ANATOMY_PARENT_ID( outset ) ; outset->anat_parent_name[0] = '\0' ; #if 0 /* if changing number of slices, can't keep slice-dependent time shifts! */ if( (nzbot!=0 || nztop!=0) && outset->taxis != NULL && outset->taxis->nsl > 0 ){ EDIT_dset_items( outset , ADN_nsl , 0 , ADN_none ) ; fprintf(stderr, "*** THD_zeropad: warning - slice-dependent time shifts have been removed!\n") ; } #else /* 31 Jan 2001: OK, lets keeps the slice-dependent time shifts (but we'll have to mangle them) -- RWCox */ if( (nzbot!=0 || nztop!=0) && outset->taxis != NULL && outset->taxis->nsl > 0 ){ int old_nsl , new_nsl , ii, nkeep,kbot,ibot ; float old_zorg_sl , *old_toff_sl , new_zorg_sl , *new_toff_sl ; /* copy current conditions to local variables */ old_nsl = outset->taxis->nsl ; old_zorg_sl = outset->taxis->zorg_sl ; old_toff_sl = (float *) malloc(sizeof(float)*old_nsl) ; memcpy( old_toff_sl , outset->taxis->toff_sl , sizeof(float)*old_nsl ) ; /* compute new values */ new_nsl = nznew ; new_zorg_sl = outset->daxes->zzorg ; /* cf. to3d.c */ new_toff_sl = (float *) malloc(sizeof(float)*new_nsl) ; for( ii=0 ; ii < new_nsl ; ii++ ) new_toff_sl[ii] = 0.0 ; /* extras are 0 */ nkeep = old_nsl ; /* how many to keep from the old list */ if( nzbot < 0 ) nkeep += nzbot ; /* lost this many at the bottom */ if( nztop < 0 ) nkeep += nztop ; /* lost this many at the top */ if( nzbot < 0 ){ kbot = -nzbot ; /* which old one to start with */ ibot = 0 ; /* and where it goes in new list */ } else { kbot = 0 ; ibot = nzbot ; } memcpy( new_toff_sl+ibot , old_toff_sl+kbot , sizeof(float)*nkeep ) ; /* set new values in dataset */ STATUS("setting new time-offsets") ; EDIT_dset_items( outset , ADN_nsl , new_nsl , ADN_toff_sl , new_toff_sl , ADN_zorg_sl , new_zorg_sl , ADN_none ) ; free(new_toff_sl) ; free(old_toff_sl) ; } #endif if( empty_flag ) RETURN(outset) ; /* 09 Feb 2001 */ /*-- now read the old dataset in, and make bricks for the new dataset --*/ STATUS("reading dataset in") ; DSET_load(inset) ; if( !DSET_LOADED(inset) ){ fprintf(stderr,"*** THD_zeropad: Can't load input dataset BRIK!\n"); DSET_delete(outset) ; RETURN( NULL ); } STATUS("padding") ; for( iv=0 ; iv < DSET_NVALS(inset) ; iv++ ){ /* create a brick of zeros */ oldim = DSET_BRICK(inset,iv) ; /* image structure of old brick */ vnew = (void*)calloc( nxnew*nynew*nznew , oldim->pixel_size ) ; /* new brick */ if( vnew == NULL ){ fprintf(stderr, "*** THD_zeropad: Can't malloc space for new sub-brick %d\n", iv) ; DSET_delete(outset) ; RETURN( NULL ); } /* macros for computing 1D subscripts from 3D indices */ #undef SNEW /* in case was defined in some stupid .h file */ #undef SOLD #define SNEW(i,j,k) ((i+nxbot)+(j+nybot)*nxnew+(k+nzbot)*nxynew) #define SOLD(i,j,k) (i+j*nxold+k*nxyold) switch( oldim->kind ){ /* copy rows of old into new */ default: break ; case MRI_byte:{ byte * bnew = (byte *) vnew, * bold = mri_data_pointer(oldim) ; for( kk=kkbot ; kk < kktop ; kk++ ) for( jj=jjbot ; jj < jjtop ; jj++ ) for( ii=iibot ; ii < iitop ; ii++ ) bnew[SNEW(ii,jj,kk)] = bold[SOLD(ii,jj,kk)] ; } break ; case MRI_short:{ short * bnew = (short *) vnew, * bold = mri_data_pointer(oldim) ; for( kk=kkbot ; kk < kktop ; kk++ ) for( jj=jjbot ; jj < jjtop ; jj++ ) for( ii=iibot ; ii < iitop ; ii++ ) bnew[SNEW(ii,jj,kk)] = bold[SOLD(ii,jj,kk)] ; } break ; case MRI_float:{ float * bnew = (float *) vnew, * bold = mri_data_pointer(oldim) ; for( kk=kkbot ; kk < kktop ; kk++ ) for( jj=jjbot ; jj < jjtop ; jj++ ) for( ii=iibot ; ii < iitop ; ii++ ) bnew[SNEW(ii,jj,kk)] = bold[SOLD(ii,jj,kk)] ; } break ; case MRI_complex:{ complex * bnew = (complex *) vnew, * bold = mri_data_pointer(oldim) ; for( kk=kkbot ; kk < kktop ; kk++ ) for( jj=jjbot ; jj < jjtop ; jj++ ) for( ii=iibot ; ii < iitop ; ii++ ) bnew[SNEW(ii,jj,kk)] = bold[SOLD(ii,jj,kk)] ; } break ; case MRI_int:{ int * bnew = (int *) vnew, * bold = mri_data_pointer(oldim) ; for( kk=kkbot ; kk < kktop ; kk++ ) for( jj=jjbot ; jj < jjtop ; jj++ ) for( ii=iibot ; ii < iitop ; ii++ ) bnew[SNEW(ii,jj,kk)] = bold[SOLD(ii,jj,kk)] ; } break ; case MRI_rgb:{ rgbyte * bnew = (rgbyte *) vnew, * bold = mri_data_pointer(oldim) ; for( kk=kkbot ; kk < kktop ; kk++ ) for( jj=jjbot ; jj < jjtop ; jj++ ) for( ii=iibot ; ii < iitop ; ii++ ) bnew[SNEW(ii,jj,kk)] = bold[SOLD(ii,jj,kk)] ; } break ; } /* end of switch on sub-brick type */ if( purge_flag) DSET_unload_one(inset,iv) ; /* 09 Feb 2001 */ EDIT_substitute_brick( outset , iv , oldim->kind , vnew ) ; } /* end of loop on sub-brick index */ #if 0 if( purge_flag ) DSET_unload(inset) ; /* 09 Feb 2001 */ #endif RETURN( outset ); } /*-------------------------------------------------------------------*/ /* Create a 1-volume dataset built from (a) a master dataset 'mset' -- for the geometry (b) an image 'imin' containing the data (c) padding amounts that reckon with the dimensional differences between imin and mset -- that is, for positive 'pad' values, imin should be bigger than mset. This function is for use in 3dAllineate. [25 Jan 2021 - RWC] *//*-----------------------------------------------------------------*/ THD_3dim_dataset * 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 ) { THD_3dim_dataset *dset ; MRI_IMAGE *dim ; int nxx,nyy,nzz ; ENTRY("THD_volume_to_dataset") ; if( mset == NULL || imin == NULL ) RETURN(NULL) ; if( !THD_filename_ok(prefix) ) prefix = "volumized" ; /* size of dataset after it gets padded */ nxx = DSET_NX(mset) + pad_xm + pad_xp ; nyy = DSET_NY(mset) + pad_ym + pad_yp ; nzz = DSET_NZ(mset) + pad_zm + pad_zp ; /* should match input image */ if( nxx != imin->nx || nyy != imin->ny || nzz != imin->nz ){ WARNING_message("THD_volume_to_dataset mismatch:\n" " Dataset = %3d %3d %3d\n" " x-pad = %3d %3d\n" " y-pad = %3d %3d\n" " z-pad = %3d %3d\n" " Image = %3d %3d %3d" , DSET_NX(mset),DSET_NY(mset),DSET_NZ(mset) , pad_xm , pad_xp , pad_ym , pad_yp , pad_zm , pad_zp , imin->nx , imin->ny , imin->nz ) ; RETURN(NULL) ; } /* create an empty (no data) dataset properly padded from mset */ dset = THD_zeropad( mset , pad_xm , pad_xp , pad_ym , pad_yp , pad_zm , pad_zp , prefix , ZPAD_EMPTY | ZPAD_IJK ) ; EDIT_dset_items( dset , ADN_nvals , 1 , ADN_ntt , 0 , ADN_datum_all , imin->kind , ADN_none ) ; /* copy the data */ dim = mri_copy(imin) ; /* shove it into the new dataset */ EDIT_substitute_brick( dset , 0 , dim->kind , mri_data_pointer(dim) ) ; /* erase the shell of the image copy */ mri_clear_and_free( dim ) ; /* get the hell out of Dodge */ RETURN(dset) ; } /*-------------------------------------------------------------------*/ /* Create a dataset from an MRI_IMARR. All sub-images must be the same dimensions. All sub-images will be converted (if necessary) to the same datum type as the [0] image. NULL is returned if inputs are stooopid. [RWC - 05 Mar 2021] *//*-----------------------------------------------------------------*/ THD_3dim_dataset * THD_imarr_to_dataset( MRI_IMARR *imar, char *prefix ) { int nx,ny,nz , nvals , vv ; MRI_IMAGE *qim ; MRI_TYPE qtyp ; THD_3dim_dataset *qset ; THD_ivec3 nxyz ; THD_fvec3 dxyz ; ENTRY("THD_imarr_to_dataset") ; if( imar == NULL || IMARR_COUNT(imar) == 0 ) RETURN(NULL) ; qim = IMARR_SUBIM(imar,0); if( qim == NULL ) RETURN(NULL) ; /* dimensions */ nx = qim->nx ; ny = qim->ny ; nz = qim->nz ; if( nx < 2 && ny < 2 && nz < 2 ) RETURN(NULL) ; LOAD_IVEC3(nxyz,nx,ny,nz) ; LOAD_FVEC3(dxyz,1.0f,1.0f,1.0f) ; qtyp = qim->kind ; nvals = IMARR_COUNT(imar) ; /* dimension check */ for( vv=1 ; vv < nvals ; vv++ ){ qim = IMARR_SUBIM(imar,vv) ; if( qim == NULL ) RETURN(NULL) ; if( qim->nx != nx || qim->ny != ny || qim->nz != nz ) RETURN(NULL) ; } /* create empty shell of dataset and set dimensions */ qset = EDIT_empty_copy(NULL) ; if( !THD_filename_ok(prefix) ) prefix = "fromIMARR" ; EDIT_dset_items( qset , ADN_datum_all , (int)qtyp , ADN_nvals , nvals , ADN_nxyz , nxyz , ADN_xyzdel , dxyz , ADN_prefix , prefix , ADN_none ) ; /* make a copy of each image, and stuff its data into the dataset, then clear out copy's data (now inside dataset) and free image struct */ for( vv=0 ; vv < nvals ; vv++ ){ qim = mri_to_mri( (int)qtyp , IMARR_SUBIM(imar,vv) ) ; EDIT_substitute_brick( qset , vv , (int)qtyp , mri_data_pointer(qim) ) ; mri_clear_and_free( qim ) ; } RETURN(qset) ; } /*-------------------------------------------------------------------*/ THD_3dim_dataset * THD_image_to_dataset( MRI_IMAGE *imin, char *prefix ) { MRI_IMARR *imar ; THD_3dim_dataset *qset ; ENTRY("THD_image_to_dataset") ; if( imin == NULL ) RETURN(NULL) ; INIT_IMARR(imar) ; ADDTO_IMARR(imar,imin) ; qset = THD_imarr_to_dataset( imar , prefix ) ; FREE_IMARR(imar) ; RETURN(qset) ; }