00001 #include "mrilib.h"
00002 
00003 /** NOT 7D SAFE **/
00004 
00005 /*----------------------------------------------------------------------*/
00006 /*! Flip a 3D image.  The (outx,outy,outz) parameters specify the
00007     direction of the output axes relative to the input axes:
00008      - +1 => +x
00009      - -1 => -x
00010      - +2 => +y
00011      - -2 => -y
00012      - +3 => +z
00013      - -3 => -z
00014 
00015     Among (outx,outy,outz), exactly one must be either +1 or -1, one
00016     must be +2 or -2, and one must be +3 or -3.  Bad inputs result in
00017     a NULL return value.
00018 
00019     If the input image doesn't have a data array, then the output
00020     image won't either.
00021 ------------------------------------------------------------------------*/
00022 
00023 MRI_IMAGE * mri_flip3D( int outx, int outy, int outz, MRI_IMAGE *inim )
00024 {
00025    MRI_IMAGE *outim ;
00026    int ii,jj,kk , dsiz , nxin,nyin,nzin , nxout,nyout,nzout ;
00027    int nxyin , nxyout , ijk_out , ijk_in ;
00028    int ax,bx,cx,dx , ay,by,cy,dy , az,bz,cz,dz , aa,bb,cc,dd ;
00029    char *inar , *outar ;
00030    float delx,dely,delz ;
00031 
00032 ENTRY("mri_flip3D") ;
00033 
00034    /* check inputs for correctness */
00035 
00036    if( inim == NULL || outx == 0 || outy == 0 || outz == 0 ) RETURN( NULL );
00037    ii = abs(outx) ; jj = abs(outy) ; kk = abs(outz) ;
00038    if( ii > 3 || jj > 3 || kk > 3 )                          RETURN( NULL );
00039    if( ii == jj || ii == kk || jj == kk )                    RETURN( NULL );
00040    if( ii+jj+kk != 6 )                                       RETURN( NULL );
00041 
00042    if( outx==1 && outy==2 && outz==3 ) RETURN( mri_copy(inim) ); /* easy case */
00043 
00044    nxin = inim->nx ;
00045    nyin = inim->ny ; nxyin = nxin*nyin ;
00046    nzin = inim->nz ;
00047 
00048    /* setup so that i_out = ax + bx*i_in + cx*j_in + dx*k_in,
00049        for i_in=0..nxin-1, j_in=0..nyin-1, k_in=0..nzin-1,
00050       and then similarly for y and z axes                    */
00051 
00052    switch( outx ){
00053      case  1:  ax=0     ; bx= 1; cx= 0; dx= 0; nxout=nxin; delx=inim->dx; break;
00054      case -1:  ax=nxin-1; bx=-1; cx= 0; dx= 0; nxout=nxin; delx=inim->dx; break;
00055      case  2:  ax=0     ; bx= 0; cx= 1; dx= 0; nxout=nyin; delx=inim->dy; break;
00056      case -2:  ax=nyin-1; bx= 0; cx=-1; dx= 0; nxout=nyin; delx=inim->dy; break;
00057      case  3:  ax=0     ; bx= 0; cx= 0; dx= 1; nxout=nzin; delx=inim->dz; break;
00058      case -3:  ax=nzin-1; bx= 0; cx= 0; dx=-1; nxout=nzin; delx=inim->dz; break;
00059      default: RETURN( NULL );
00060    }
00061    switch( outy ){
00062      case  1:  ay=0     ; by= 1; cy= 0; dy= 0; nyout=nxin; dely=inim->dx; break;
00063      case -1:  ay=nxin-1; by=-1; cy= 0; dy= 0; nyout=nxin; dely=inim->dx; break;
00064      case  2:  ay=0     ; by= 0; cy= 1; dy= 0; nyout=nyin; dely=inim->dy; break;
00065      case -2:  ay=nyin-1; by= 0; cy=-1; dy= 0; nyout=nyin; dely=inim->dy; break;
00066      case  3:  ay=0     ; by= 0; cy= 0; dy= 1; nyout=nzin; dely=inim->dz; break;
00067      case -3:  ay=nzin-1; by= 0; cy= 0; dy=-1; nyout=nzin; dely=inim->dz; break;
00068      default: RETURN( NULL );
00069    }
00070    switch( outz ){
00071      case  1:  az=0     ; bz= 1; cz= 0; dz= 0; nzout=nxin; delz=inim->dx; break;
00072      case -1:  az=nxin-1; bz=-1; cz= 0; dz= 0; nzout=nxin; delz=inim->dx; break;
00073      case  2:  az=0     ; bz= 0; cz= 1; dz= 0; nzout=nyin; delz=inim->dy; break;
00074      case -2:  az=nyin-1; bz= 0; cz=-1; dz= 0; nzout=nyin; delz=inim->dy; break;
00075      case  3:  az=0     ; bz= 0; cz= 0; dz= 1; nzout=nzin; delz=inim->dz; break;
00076      case -3:  az=nzin-1; bz= 0; cz= 0; dz=-1; nzout=nzin; delz=inim->dz; break;
00077      default: RETURN( NULL );
00078    }
00079    nxyout = nxout*nyout ;
00080 
00081    /* 3D index ijk_out = i_out + nxout*j_out + nxyout*k_out
00082 
00083                        = (ax + bx*i_in + cx*j_in + dx*k_in)
00084                         +(ay + by*i_in + cy*j_in + dy*k_in)*nxout
00085                         +(az + bz*i_in + cz*j_in + dz*k_in)*nxyout
00086 
00087                        = (ax+ay*nxout+az*nxyout)
00088                         +(bx+by*nxout+bz*nxyout)*i_in
00089                         +(cx+cy*nxout+cz*nxyout)*j_in
00090                         +(dx+dy*nxout+dz*nxyout)*k_in
00091 
00092                        = aa + bb*i_in + cc*j_in + dd*k_in  */
00093 
00094    inar  = mri_data_pointer( inim ) ;
00095    outim = mri_new_vol( nxout,nyout,nzout , inim->kind ) ;
00096 
00097    outim->dx = delx ; outim->dy = dely ; outim->dz = delz ;
00098 
00099    outar = mri_data_pointer( outim ) ;
00100 
00101    if( inar == NULL ){       /* empty input ==> empty output */
00102      free(outar) ; mri_fix_data_pointer(NULL,outim ); RETURN(outim);
00103    }
00104 
00105    dsiz  = outim->pixel_size ;          /* size of each voxel in bytes */
00106 
00107    aa = (ax+ay*nxout+az*nxyout)*dsiz ;  /* same as aa, etc. in above */
00108    bb = (bx+by*nxout+bz*nxyout)*dsiz ;  /* comment, but scaled to be */
00109    cc = (cx+cy*nxout+cz*nxyout)*dsiz ;  /* address offset in bytes   */
00110    dd = (dx+dy*nxout+dz*nxyout)*dsiz ;
00111 
00112    for( ijk_in=kk=0 ; kk < nzin ; kk++ ){
00113      for( jj=0 ; jj < nyin ; jj++ ){
00114        for( ii=0 ; ii < nxin ; ii++,ijk_in+=dsiz ){
00115          ijk_out = aa + bb*ii + cc*jj + dd*kk ;
00116          memcpy( outar+ijk_out , inar+ijk_in , dsiz ) ;
00117        }
00118      }
00119    }
00120 
00121    RETURN(outim) ;
00122 }