00001 #include "mrilib.h"
00002 
00003 /*-------------------------------------------------------------------
00004   Create a new dataset from an old one, with zero padding around
00005   the edges.  For example,
00006     add_I = number of zero planes to add at inferior edge
00007             (if < 0, number of data planes to cut off inferior edge)
00008 
00009   09 Feb 2001 - added "flag" input, which is the OR of
00010     ZPAD_EMPTY = produce only an padded "empty copy" of the input
00011     ZPAD_PURGE = purge input dataset bricks after they are copied
00012     ZPAD_MM    = increments are mm instead of slice counts
00013                  (at least 'add_?' mm will be added/subtracted)
00014 
00015   14 May 2002: if inputs crops are all zero, return something anyway
00016 ---------------------------------------------------------------------*/
00017 
00018 THD_3dim_dataset * THD_zeropad( THD_3dim_dataset * inset ,
00019                                 int add_I , int add_S , int add_A ,
00020                                 int add_P , int add_L , int add_R ,
00021                                 char * prefix , int flag )
00022 {
00023    THD_3dim_dataset *outset ;
00024    int nxold,nyold,nzold , nxnew,nynew,nznew , nxyold,nxynew ,
00025        nxbot=0,nxtop=0 , nybot=0,nytop=0 , nzbot=0,nztop=0    ;
00026    int ii,jj,kk , iv , iibot,iitop , jjbot,jjtop , kkbot,kktop ;
00027 
00028    int empty_flag = (flag & ZPAD_EMPTY) ;  /* 09 Feb 2001 */
00029    int purge_flag = (flag & ZPAD_PURGE) ;  /* 09 Feb 2001 */
00030    int mm_flag    = (flag & ZPAD_MM   ) ;  /* 13 Feb 2001 */
00031 
00032    THD_ivec3 iv_nxyz ;
00033    THD_fvec3 fv_xyzorg ;
00034 
00035    MRI_IMAGE * oldim ;
00036    void * vnew ;
00037 
00038 ENTRY("THD_zeropad") ;
00039 
00040    /*-- check inputs --*/
00041 
00042    if( !ISVALID_DSET(inset) ) RETURN( NULL ) ;
00043 
00044    if( add_I==0 && add_S==0 && add_P==0 &&
00045        add_A==0 && add_L==0 && add_R==0    ){
00046 
00047       fprintf(stderr,"++ THD_zeropad: all pad values are zero!\n") ;
00048 
00049       outset = EDIT_full_copy( inset , prefix ) ;  /* 14 May 2002 */
00050       RETURN( outset );
00051    }
00052 
00053    if( !THD_filename_ok(prefix) ) prefix = "zeropad" ;
00054 
00055    /*-- map add_? values into dataset xyz coordinate directions --*/
00056 
00057    nxold = DSET_NX(inset) ;
00058    nyold = DSET_NY(inset) ;
00059    nzold = DSET_NZ(inset) ;
00060 
00061    /* comput n?top and n?bot, the number of planes to add at
00062       the top and bottom of the ? direction, for ? = x, y, or z */
00063 
00064    switch( inset->daxes->xxorient ){
00065       default:
00066         fprintf(stderr,"*** THD_zeropad: Unknown orientation codes!\n") ;
00067         RETURN( NULL );
00068 
00069       case ORI_R2L_TYPE: nxtop = add_L ; nxbot = add_R ; break ;
00070       case ORI_L2R_TYPE: nxtop = add_R ; nxbot = add_L ; break ;
00071       case ORI_P2A_TYPE: nxtop = add_A ; nxbot = add_P ; break ;
00072       case ORI_A2P_TYPE: nxtop = add_P ; nxbot = add_A ; break ;
00073       case ORI_I2S_TYPE: nxtop = add_S ; nxbot = add_I ; break ;
00074       case ORI_S2I_TYPE: nxtop = add_I ; nxbot = add_S ; break ;
00075    }
00076 
00077    switch( inset->daxes->yyorient ){
00078       default:
00079         fprintf(stderr,"*** THD_zeropad: Unknown orientation codes!\n") ;
00080         RETURN( NULL );
00081 
00082       case ORI_R2L_TYPE: nytop = add_L ; nybot = add_R ; break ;
00083       case ORI_L2R_TYPE: nytop = add_R ; nybot = add_L ; break ;
00084       case ORI_P2A_TYPE: nytop = add_A ; nybot = add_P ; break ;
00085       case ORI_A2P_TYPE: nytop = add_P ; nybot = add_A ; break ;
00086       case ORI_I2S_TYPE: nytop = add_S ; nybot = add_I ; break ;
00087       case ORI_S2I_TYPE: nytop = add_I ; nybot = add_S ; break ;
00088    }
00089 
00090    switch( inset->daxes->zzorient ){
00091       default:
00092         fprintf(stderr,"*** THD_zeropad: Unknown orientation codes!\n") ;
00093         RETURN( NULL );
00094 
00095       case ORI_R2L_TYPE: nztop = add_L ; nzbot = add_R ; break ;
00096       case ORI_L2R_TYPE: nztop = add_R ; nzbot = add_L ; break ;
00097       case ORI_P2A_TYPE: nztop = add_A ; nzbot = add_P ; break ;
00098       case ORI_A2P_TYPE: nztop = add_P ; nzbot = add_A ; break ;
00099       case ORI_I2S_TYPE: nztop = add_S ; nzbot = add_I ; break ;
00100       case ORI_S2I_TYPE: nztop = add_I ; nzbot = add_S ; break ;
00101    }
00102 
00103    /* 13 Feb 2001: round to millimeters? */
00104 
00105 #undef  RMM
00106 #define RMM(n,d)                                                           \
00107   do{      if( (n) > 0 ) (n) = (int)( (n)/fabs(d) + 0.999 ) ;              \
00108       else if( (n) < 0 ) (n) = (int)( (n)/fabs(d) - 0.999 ) ; } while(0) ;
00109 
00110    if( mm_flag ){
00111       RMM(nxtop,inset->daxes->xxdel) ; RMM(nxbot,inset->daxes->xxdel) ;
00112       RMM(nytop,inset->daxes->yydel) ; RMM(nybot,inset->daxes->yydel) ;
00113       RMM(nztop,inset->daxes->zzdel) ; RMM(nzbot,inset->daxes->zzdel) ;
00114    }
00115 
00116    nxnew = nxold + nxbot + nxtop ;  /* dimensions of new bricks */
00117    nynew = nyold + nybot + nytop ;
00118    nznew = nzold + nzbot + nztop ;
00119 
00120    nxyold = nxold * nyold ;         /* for computing subscripts */
00121    nxynew = nxnew * nynew ;
00122 
00123    iibot = MAX(0,-nxbot) ; iitop = MIN(nxold,nxold+nxtop) ;  /* range of data */
00124    jjbot = MAX(0,-nybot) ; jjtop = MIN(nyold,nyold+nytop) ;  /* in old dataset */
00125    kkbot = MAX(0,-nzbot) ; kktop = MIN(nzold,nzold+nztop) ;
00126 
00127    if( nxnew < 1 || iibot > iitop ||   /* check for reasonable sizes */
00128        nynew < 1 || jjbot > jjtop ||   /* and ranges of dataset     */
00129        nznew < 1 || kkbot > kktop   ){
00130 
00131       fprintf(stderr,"*** THD_zeropad: Can't cut dataset down too much!\n") ;
00132       RETURN( NULL );
00133    }
00134 
00135 #if 0
00136    if( nxnew < 2 || iibot >= iitop ||   /* check for reasonable sizes */
00137        nynew < 2 || jjbot >= jjtop ||   /* and ranges of dataset     */
00138        nznew < 2 || kkbot >= kktop   ){
00139 
00140       fprintf(stderr,"*** WARNING - THD_zeropad: dataset cut down to %dx%dx%d\n",
00141                       nxnew,nynew,nznew) ;
00142    }
00143 #endif
00144 
00145    /*-- create the shell of the new dataset --*/
00146 
00147    outset = EDIT_empty_copy( inset ) ;
00148 
00149    LOAD_IVEC3( iv_nxyz , nxnew,nynew,nznew ) ;
00150 
00151    LOAD_FVEC3( fv_xyzorg, inset->daxes->xxorg - nxbot * inset->daxes->xxdel,
00152                           inset->daxes->yyorg - nybot * inset->daxes->yydel,
00153                           inset->daxes->zzorg - nzbot * inset->daxes->zzdel );
00154 
00155 STATUS("setting new dimensions") ;
00156 
00157    EDIT_dset_items( outset ,
00158                        ADN_prefix , prefix    ,
00159                        ADN_nxyz   , iv_nxyz   ,
00160                        ADN_xyzorg , fv_xyzorg ,
00161                     ADN_none ) ;
00162 
00163    /* Changing dimensions means old anat parent is no longer valid! */
00164 
00165    EDIT_ZERO_ANATOMY_PARENT_ID( outset ) ;
00166    outset->anat_parent_name[0] = '\0' ;
00167 
00168 #if 0
00169    /* if changing number of slices, can't keep slice-dependent time shifts! */
00170 
00171    if( (nzbot!=0 || nztop!=0) && outset->taxis != NULL && outset->taxis->nsl > 0 ){
00172       EDIT_dset_items( outset , ADN_nsl , 0 , ADN_none ) ;
00173       fprintf(stderr,
00174               "*** THD_zeropad: warning - slice-dependent time shifts have been removed!\n") ;
00175    }
00176 #else
00177    /* 31 Jan 2001: OK, lets keeps the slice-dependent time shifts
00178                    (but we'll have to mangle them) -- RWCox      */
00179 
00180    if( (nzbot!=0 || nztop!=0) && outset->taxis != NULL && outset->taxis->nsl > 0 ){
00181       int old_nsl , new_nsl , ii, nkeep,kbot,ibot ;
00182       float old_zorg_sl , *old_toff_sl , new_zorg_sl , *new_toff_sl ;
00183 
00184       /* copy current conditions to local variables */
00185 
00186       old_nsl     = outset->taxis->nsl ;
00187       old_zorg_sl = outset->taxis->zorg_sl ;
00188       old_toff_sl = (float *) malloc(sizeof(float)*old_nsl) ;
00189       memcpy( old_toff_sl , outset->taxis->toff_sl , sizeof(float)*old_nsl ) ;
00190 
00191       /* compute new values */
00192 
00193       new_nsl     = nznew ;
00194       new_zorg_sl = outset->daxes->zzorg ;                      /* cf. to3d.c */
00195       new_toff_sl = (float *) malloc(sizeof(float)*new_nsl) ;
00196       for( ii=0 ; ii < new_nsl ; ii++ ) new_toff_sl[ii] = 0.0 ; /* extras are 0 */
00197 
00198       nkeep = old_nsl ;                 /* how many to keep from the old list */
00199       if( nzbot < 0 ) nkeep += nzbot ;  /* lost this many at the bottom */
00200       if( nztop < 0 ) nkeep += nztop ;  /* lost this many at the top */
00201 
00202       if( nzbot < 0 ){
00203          kbot = -nzbot ;   /* which old one to start with */
00204          ibot = 0 ;        /* and where it goes in new list */
00205       } else {
00206          kbot = 0 ;
00207          ibot = nzbot ;
00208       }
00209 
00210       memcpy( new_toff_sl+ibot , old_toff_sl+kbot , sizeof(float)*nkeep ) ;
00211 
00212       /* set new values in dataset */
00213 
00214 STATUS("setting new time-offsets") ;
00215 
00216       EDIT_dset_items( outset ,
00217                          ADN_nsl     , new_nsl     ,
00218                          ADN_toff_sl , new_toff_sl ,
00219                          ADN_zorg_sl , new_zorg_sl ,
00220                        ADN_none ) ;
00221 
00222       free(new_toff_sl) ; free(old_toff_sl) ;
00223    }
00224 #endif
00225 
00226    if( empty_flag ) RETURN(outset) ;  /* 09 Feb 2001 */
00227 
00228    /*-- now read the old dataset in, and make bricks for the new dataset --*/
00229 
00230 STATUS("reading dataset in") ;
00231 
00232    DSET_load(inset) ;
00233    if( !DSET_LOADED(inset) ){
00234       fprintf(stderr,"*** THD_zeropad: Can't load input dataset BRIK!\n");
00235       DSET_delete(outset) ;
00236       RETURN( NULL );
00237    }
00238 
00239 STATUS("padding") ;
00240 
00241    for( iv=0 ; iv < DSET_NVALS(inset) ; iv++ ){
00242 
00243       /* create a brick of zeros */
00244 
00245       oldim = DSET_BRICK(inset,iv) ;  /* image structure of old brick */
00246 
00247       vnew  = (void*)calloc( nxnew*nynew*nznew , 
00248                              oldim->pixel_size ) ; /* new brick */
00249       if( vnew == NULL ){
00250          fprintf(stderr,
00251                  "*** THD_zeropad: Can't malloc space for new sub-brick %d\n",
00252                  iv) ;
00253          DSET_delete(outset) ; RETURN( NULL );
00254       }
00255 
00256       /* macros for computing 1D subscripts from 3D indices */
00257 
00258 #undef  SNEW  /* in case was defined in some stupid .h file */
00259 #undef  SOLD
00260 #define SNEW(i,j,k) ((i+nxbot)+(j+nybot)*nxnew+(k+nzbot)*nxynew)
00261 #define SOLD(i,j,k) (i+j*nxold+k*nxyold)
00262 
00263       switch( oldim->kind ){  /* copy rows of old into new */
00264 
00265          case MRI_byte:{
00266             byte * bnew = (byte *) vnew, * bold = mri_data_pointer(oldim) ;
00267             for( kk=kkbot ; kk < kktop ; kk++ )
00268                for( jj=jjbot ; jj < jjtop ; jj++ )
00269                   for( ii=iibot ; ii < iitop ; ii++ )
00270                      bnew[SNEW(ii,jj,kk)] = bold[SOLD(ii,jj,kk)] ;
00271          }
00272          break ;
00273 
00274          case MRI_short:{
00275             short * bnew = (short *) vnew, * bold = mri_data_pointer(oldim) ;
00276             for( kk=kkbot ; kk < kktop ; kk++ )
00277                for( jj=jjbot ; jj < jjtop ; jj++ )
00278                   for( ii=iibot ; ii < iitop ; ii++ )
00279                      bnew[SNEW(ii,jj,kk)] = bold[SOLD(ii,jj,kk)] ;
00280          }
00281          break ;
00282 
00283          case MRI_float:{
00284             float * bnew = (float *) vnew, * bold = mri_data_pointer(oldim) ;
00285             for( kk=kkbot ; kk < kktop ; kk++ )
00286                for( jj=jjbot ; jj < jjtop ; jj++ )
00287                   for( ii=iibot ; ii < iitop ; ii++ )
00288                      bnew[SNEW(ii,jj,kk)] = bold[SOLD(ii,jj,kk)] ;
00289          }
00290          break ;
00291 
00292          case MRI_complex:{
00293             complex * bnew = (complex *) vnew, * bold = mri_data_pointer(oldim) ;
00294             for( kk=kkbot ; kk < kktop ; kk++ )
00295                for( jj=jjbot ; jj < jjtop ; jj++ )
00296                   for( ii=iibot ; ii < iitop ; ii++ )
00297                      bnew[SNEW(ii,jj,kk)] = bold[SOLD(ii,jj,kk)] ;
00298          }
00299          break ;
00300 
00301          case MRI_int:{
00302             int * bnew = (int *) vnew, * bold = mri_data_pointer(oldim) ;
00303             for( kk=kkbot ; kk < kktop ; kk++ )
00304                for( jj=jjbot ; jj < jjtop ; jj++ )
00305                   for( ii=iibot ; ii < iitop ; ii++ )
00306                      bnew[SNEW(ii,jj,kk)] = bold[SOLD(ii,jj,kk)] ;
00307          }
00308          break ;
00309 
00310          case MRI_rgb:{
00311             rgbyte * bnew = (rgbyte *) vnew, * bold = mri_data_pointer(oldim) ;
00312             for( kk=kkbot ; kk < kktop ; kk++ )
00313                for( jj=jjbot ; jj < jjtop ; jj++ )
00314                   for( ii=iibot ; ii < iitop ; ii++ )
00315                      bnew[SNEW(ii,jj,kk)] = bold[SOLD(ii,jj,kk)] ;
00316          }
00317          break ;
00318 
00319       } /* end of switch on sub-brick type */
00320 
00321       if( purge_flag) DSET_unload_one(inset,iv) ; /* 09 Feb 2001 */
00322 
00323       EDIT_substitute_brick( outset , iv , oldim->kind , vnew ) ;
00324 
00325    } /* end of loop on sub-brick index */
00326 
00327 #if 0
00328    if( purge_flag ) DSET_unload(inset) ; /* 09 Feb 2001 */
00329 #endif
00330 
00331    RETURN( outset );
00332 }