00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #include "mrilib.h"
00008 
00009 /*** Not 7D safe ***/
00010 
00011 #define MAX_NUP 32
00012 
00013 static void (*usammer)(int,int,float *,float *) = upsample_7 ;
00014 
00015 static int usammer_mode = 7 ;
00016 
00017 static void upsample_1by2( int, byte *, byte * ) ;  /* at end of file */
00018 static void upsample_1by3( int, byte *, byte * ) ;
00019 static void upsample_1by4( int, byte *, byte * ) ;
00020 
00021 static MRI_IMAGE * mri_dup2D_rgb4( MRI_IMAGE *) ;   /* 14 Mar 2002 */
00022 static MRI_IMAGE * mri_dup2D_rgb3( MRI_IMAGE *) ;   /* 14 Mar 2002 */
00023 static MRI_IMAGE * mri_dup2D_rgb2( MRI_IMAGE *) ;   /* 22 Mar 2002 */
00024 
00025 static MRI_IMAGE * mri_dup2D_rgb_NN( MRI_IMAGE *, int);/* 22 Feb 2004 [rickr] */
00026 
00027 /*-----------------------------------------------------------------------------
00028   Set resample mode for mri_dup2D (1 or 7)
00029 -------------------------------------------------------------------------------*/
00030 
00031 void mri_dup2D_mode( int mm )
00032 {
00033    switch( mm ){
00034       case 1:  usammer = upsample_1 ; break ;
00035      default:  usammer = upsample_7 ; break ;
00036    }
00037    usammer_mode = mm ;
00038 }
00039 
00040 /*------------------------------------------------------------------------------
00041    Blow up a 2D image nup times, using 7th order polynomial for interpolation.
00042 --------------------------------------------------------------------------------*/
00043 
00044 MRI_IMAGE * mri_dup2D( int nup , MRI_IMAGE *imin )
00045 {
00046    MRI_IMAGE *flim , *newim ;
00047    float     *flar , *newar , *cold , *cnew ;
00048    int nx,ny , nxup,nyup , ii,jj,kk, NNmode = 0 ;
00049 
00050 ENTRY("mri_dup2D") ;
00051    /*-- sanity checks --*/
00052 
00053    if( nup < 1 || nup > MAX_NUP || imin == NULL ) RETURN( NULL );
00054 
00055    if( nup == 1 ){ newim = mri_to_mri( imin->kind, imin ); RETURN(newim); }
00056 
00057    /* does the user want neighbor interpolation?     22 Feb 2004 [rickr] */
00058    if ( AFNI_yesenv("AFNI_IMAGE_ZOOM_NN") ) {
00059       mri_dup2D_mode(0);
00060       NNmode = 1;
00061    }
00062 
00063    /*-- complex-valued images: do each part separately --*/
00064 
00065    if( imin->kind == MRI_complex ){
00066       MRI_IMARR *impair ; MRI_IMAGE *rim, *iim, *tim ;
00067 
00068       impair = mri_complex_to_pair( imin ) ;
00069       if( impair == NULL ){
00070          fprintf(stderr,"*** mri_complex_to_pair fails in mri_dup2D!\n"); EXIT(1);
00071       }
00072       rim = IMAGE_IN_IMARR(impair,0) ;
00073       iim = IMAGE_IN_IMARR(impair,1) ;  FREE_IMARR(impair) ;
00074       tim = mri_dup2D( nup, rim ); mri_free( rim ); rim = tim ;
00075       tim = mri_dup2D( nup, iim ); mri_free( iim ); iim = tim ;
00076       newim = mri_pair_to_complex( rim , iim ) ;
00077       mri_free( rim ) ; mri_free( iim ) ;
00078       MRI_COPY_AUX(newim,imin) ;
00079       RETURN(newim) ;
00080    }
00081 
00082    /*-- 14 Mar 2002: RGB image up by 2..4, all colors at once --*/
00083 
00084    if( imin->kind == MRI_rgb ){
00085      MRI_IMAGE *qqim=NULL ;
00086      if ( NNmode )
00087          qqim = mri_dup2D_rgb_NN(imin, nup);  /* 22 Feb 2004 [rickr] */
00088      else {
00089        switch(nup){
00090          case 4: qqim = mri_dup2D_rgb4(imin); break; /* special purpose fast codes */
00091          case 3: qqim = mri_dup2D_rgb3(imin); break; /* using fixed pt arithmetic  */
00092          case 2: qqim = mri_dup2D_rgb2(imin); break;
00093 
00094         /*-- other factors: do each color separately as a byte image --*/
00095         default:{
00096            MRI_IMARR *imtriple ; MRI_IMAGE *rim, *gim, *bim, *tim ;
00097 
00098            imtriple = mri_rgb_to_3byte( imin ) ;
00099            if( imtriple == NULL ){
00100              fprintf(stderr,"*** mri_rgb_to_3byte fails in mri_dup2D!\n"); RETURN(NULL);
00101            }
00102            rim = IMAGE_IN_IMARR(imtriple,0) ;
00103            gim = IMAGE_IN_IMARR(imtriple,1) ;
00104            bim = IMAGE_IN_IMARR(imtriple,2) ; FREE_IMARR(imtriple) ;
00105            tim = mri_dup2D( nup, rim ); mri_free(rim); rim = tim;
00106            tim = mri_dup2D( nup, gim ); mri_free(gim); gim = tim;
00107            tim = mri_dup2D( nup, bim ); mri_free(bim); bim = tim;
00108            newim = mri_3to_rgb( rim, gim, bim ) ;
00109            mri_free(rim) ; mri_free(gim) ; mri_free(bim) ;
00110            MRI_COPY_AUX(newim,imin) ;
00111            qqim = newim ;
00112          }
00113          break ;
00114        }
00115      }
00116 
00117      RETURN(qqim) ;
00118    }
00119 
00120    /*-- Special case: byte-valued image upsampled by 2/3/4 [13 Mar 2002] --*/
00121 
00122    if( imin->kind == MRI_byte && nup <= 4 ){
00123      void (*usbyte)(int,byte *,byte *) = NULL ;
00124      byte *bar=MRI_BYTE_PTR(imin) , *bnew , *cold, *cnew ;
00125      nx = imin->nx; ny = imin->ny; nxup = nx*nup; nyup = ny*nup ;
00126      newim = mri_new( nxup,nyup , MRI_byte ); bnew = MRI_BYTE_PTR(newim);
00127      switch( nup ){
00128        case 2: usbyte = upsample_1by2 ; break ;  /* special fast codes */
00129        case 3: usbyte = upsample_1by3 ; break ;
00130        case 4: usbyte = upsample_1by4 ; break ;
00131      }
00132      for( jj=0 ; jj < ny ; jj++ )                      /* upsample rows */
00133        usbyte( nx , bar+jj*nx , bnew+jj*nxup ) ;
00134      cold = (byte *) malloc( sizeof(byte) * ny   ) ;
00135      cnew = (byte *) malloc( sizeof(byte) * nyup ) ;
00136      for( ii=0 ; ii < nxup ; ii++ ){                   /* upsample cols */
00137        for( jj=0 ; jj < ny ; jj++ ) cold[jj] = bnew[ii+jj*nxup] ;
00138        usbyte( ny , cold , cnew ) ;
00139        for( jj=0 ; jj < nyup ; jj++ ) bnew[ii+jj*nxup] = cnew[jj] ;
00140      }
00141      free(cold); free(cnew); MRI_COPY_AUX(newim,imin); RETURN(newim);
00142    }
00143 
00144    /*-- otherwise, make sure we operate on a float image --*/
00145 
00146    if( imin->kind == MRI_float ) flim = imin ;
00147    else                          flim = mri_to_float( imin ) ;
00148 
00149    flar = MRI_FLOAT_PTR(flim) ;
00150 
00151    nx = flim->nx ; ny = flim->ny ; nxup = nx*nup ; nyup = ny*nup ;
00152    newim = mri_new( nxup , nyup , MRI_float ) ;
00153    newar = MRI_FLOAT_PTR(newim) ;
00154 
00155    /*-- upsample rows --*/
00156 
00157    for( jj=0 ; jj < ny ; jj++ )
00158       usammer( nup , nx , flar + jj*nx , newar + jj*nxup ) ;
00159 
00160    if( flim != imin ) mri_free(flim) ;
00161 
00162    /*-- upsample columns --*/
00163 
00164    cold = (float *) malloc( sizeof(float) * ny ) ;
00165    cnew = (float *) malloc( sizeof(float) * nyup ) ;
00166    if( cold == NULL || cnew == NULL ){
00167       fprintf(stderr,"*** mri_dup2D malloc failure!\n"); EXIT(1);
00168    }
00169 
00170    for( ii=0 ; ii < nxup ; ii++ ){
00171       for( jj=0 ; jj < ny ; jj++ ) cold[jj] = newar[ii + jj*nxup] ;
00172       usammer( nup , ny , cold , cnew ) ;
00173       for( jj=0 ; jj < nyup ; jj++ ) newar[ii+jj*nxup] = cnew[jj] ;
00174    }
00175 
00176    free(cold) ; free(cnew) ;
00177 
00178    /*-- type convert output, if necessary --*/
00179 
00180    switch( imin->kind ){
00181 
00182       case MRI_byte:{
00183          byte * bar ; MRI_IMAGE * bim ; float fmin , fmax ;
00184 
00185          bim = mri_new( nxup,nyup , MRI_byte ) ; bar = MRI_BYTE_PTR(bim) ;
00186          fmin = mri_min(imin) ; fmax = mri_max(imin) ;
00187          for( ii=0 ; ii < newim->nvox ; ii++ )
00188             bar[ii] =  (newar[ii] < fmin) ? fmin
00189                      : (newar[ii] > fmax) ? fmax : newar[ii] ;
00190          mri_free(newim) ; newim = bim ;
00191       }
00192       break ;
00193 
00194       case MRI_short:{
00195          short * sar ; MRI_IMAGE * sim ; float fmin , fmax ;
00196 
00197          sim = mri_new( nxup,nyup , MRI_short ) ; sar = MRI_SHORT_PTR(sim) ;
00198          fmin = mri_min(imin) ; fmax = mri_max(imin) ;
00199          for( ii=0 ; ii < newim->nvox ; ii++ )
00200             sar[ii] =  (newar[ii] < fmin) ? fmin
00201                      : (newar[ii] > fmax) ? fmax : newar[ii] ;
00202          mri_free(newim) ; newim = sim ;
00203       }
00204       break ;
00205 
00206       case MRI_float:{
00207          float fmin , fmax ;
00208 
00209          fmin = mri_min(imin) ; fmax = mri_max(imin) ;
00210          for( ii=0 ; ii < newim->nvox ; ii++ )
00211                  if( newar[ii] < fmin ) newar[ii] = fmin ;
00212             else if( newar[ii] > fmax ) newar[ii] = fmax ;
00213       }
00214    }
00215 
00216    /*-- finito --*/
00217 
00218    MRI_COPY_AUX(newim,imin) ;
00219    RETURN( newim );
00220 }
00221 
00222 /*======================================================================*/
00223 
00224 /*-- seventh order Lagrange polynomials --*/
00225 
00226 #define S_M3(x) (x*(x*x-1.0)*(x*x-4.0)*(x-3.0)*(4.0-x)*0.0001984126984)
00227 #define S_M2(x) (x*(x*x-1.0)*(x-2.0)*(x*x-9.0)*(x-4.0)*0.001388888889)
00228 #define S_M1(x) (x*(x-1.0)*(x*x-4.0)*(x*x-9.0)*(4.0-x)*0.004166666667)
00229 #define S_00(x) ((x*x-1.0)*(x*x-4.0)*(x*x-9.0)*(x-4.0)*0.006944444444)
00230 #define S_P1(x) (x*(x+1.0)*(x*x-4.0)*(x*x-9.0)*(4.0-x)*0.006944444444)
00231 #define S_P2(x) (x*(x*x-1.0)*(x+2.0)*(x*x-9.0)*(x-4.0)*0.004166666667)
00232 #define S_P3(x) (x*(x*x-1.0)*(x*x-4.0)*(x+3.0)*(4.0-x)*0.001388888889)
00233 #define S_P4(x) (x*(x*x-1.0)*(x*x-4.0)*(x*x-9.0)*0.0001984126984)
00234 
00235 #ifdef ZFILL
00236 #  define FINS(i) ( ((i)<0 || (i)>=nar) ? 0.0 : far[(i)] )
00237 #else
00238 #  define FINS(i) ( ((i)<0) ? far[0] : ((i)>=nar) ? far[nar-1] : far[(i)] )
00239 #endif
00240 
00241 #define INT7(k,i) (  fm3[k] * far[i-3] + fm2[k] * far[i-2] \
00242                    + fm1[k] * far[i-1] + f00[k] * far[i  ] \
00243                    + fp1[k] * far[i+1] + fp2[k] * far[i+2] \
00244                    + fp3[k] * far[i+3] + fp4[k] * far[i+4]  )
00245 
00246 #define FINT7(k,i) (  fm3[k] * FINS(i-3) + fm2[k] * FINS(i-2) \
00247                     + fm1[k] * FINS(i-1) + f00[k] * FINS(i  ) \
00248                     + fp1[k] * FINS(i+1) + fp2[k] * FINS(i+2) \
00249                     + fp3[k] * FINS(i+3) + fp4[k] * FINS(i+4)  )
00250 
00251 /*----------------------------------------------------------------------------
00252   Up sample an array far[0..nar-1] nup times to produce fout[0..nar*nup-1].
00253   Uses 7th order polynomial interpolation.
00254 ------------------------------------------------------------------------------*/
00255 
00256 void upsample_7( int nup , int nar , float * far , float * fout )
00257 {
00258    int kk,ii , ibot,itop ;
00259    static int nupold = -1 ;
00260    static float fm3[MAX_NUP], fm2[MAX_NUP], fm1[MAX_NUP], f00[MAX_NUP],
00261                 fp1[MAX_NUP], fp2[MAX_NUP], fp3[MAX_NUP], fp4[MAX_NUP] ;
00262 
00263    /*-- sanity checks --*/
00264 
00265    if( nup < 1 || nup > MAX_NUP || nar < 2 || far == NULL || fout == NULL ) return ;
00266 
00267    if( nup == 1 ){ memcpy( fout, far, sizeof(float)*nar ); return; }
00268 
00269    /*-- initialize interpolation coefficient, if nup has changed --*/
00270 
00271    if( nup != nupold ){
00272      float val ;
00273      for( kk=0 ; kk < nup ; kk++ ){
00274        val = ((float)kk) / ((float)nup) ;
00275        fm3[kk] = S_M3(val); fm2[kk] = S_M2(val); fm1[kk] = S_M1(val);
00276        f00[kk] = S_00(val); fp1[kk] = S_P1(val); fp2[kk] = S_P2(val);
00277        fp3[kk] = S_P3(val); fp4[kk] = S_P4(val);
00278      }
00279      nupold = nup ;
00280    }
00281 
00282    /*-- interpolate the intermediate places --*/
00283 
00284    ibot = 3 ; itop = nar-5 ;
00285 
00286    switch( nup ){
00287       default:
00288         for( ii=ibot ; ii <= itop ; ii++ )
00289           for( kk=0 ; kk < nup ; kk++ ) fout[kk+ii*nup] = INT7(kk,ii) ;
00290       break ;
00291 
00292       case 2:
00293         for( ii=ibot ; ii <= itop ; ii++ ){
00294           fout[ii*nup]   = INT7(0,ii) ; fout[ii*nup+1] = INT7(1,ii) ;
00295         }
00296       break ;
00297 
00298       case 3:
00299         for( ii=ibot ; ii <= itop ; ii++ ){
00300           fout[ii*nup]   = INT7(0,ii) ; fout[ii*nup+1] = INT7(1,ii) ;
00301           fout[ii*nup+2] = INT7(2,ii) ;
00302         }
00303       break ;
00304 
00305       case 4:
00306         for( ii=ibot ; ii <= itop ; ii++ ){
00307           fout[ii*nup]   = INT7(0,ii) ; fout[ii*nup+1] = INT7(1,ii) ;
00308           fout[ii*nup+2] = INT7(2,ii) ; fout[ii*nup+3] = INT7(3,ii) ;
00309         }
00310       break ;
00311    }
00312 
00313    /*-- interpolate the outside edges --*/
00314 
00315    for( ii=0 ; ii < ibot ; ii++ )
00316      for( kk=0 ; kk < nup ; kk++ ) fout[kk+ii*nup] = FINT7(kk,ii) ;
00317 
00318    for( ii=itop+1 ; ii < nar ; ii++ )
00319      for( kk=0 ; kk < nup ; kk++ ) fout[kk+ii*nup] =  FINT7(kk,ii) ;
00320 
00321    return ;
00322 }
00323 
00324 /*======================================================================*/
00325 
00326 #define INT1(k,i)  (f00[k]*far[i]  + fp1[k]*far[i+1] )
00327 
00328 #define FINT1(k,i) (f00[k]*FINS(i) + fp1[k]*FINS(i+1))
00329 
00330 /*----------------------------------------------------------------------------
00331   Up sample an array far[0..nar-1] nup times to produce fout[0..nar*nup-1].
00332   Uses linear polynomial interpolation.
00333 ------------------------------------------------------------------------------*/
00334 
00335 void upsample_1( int nup , int nar , float * far , float * fout )
00336 {
00337    int kk,ii , ibot,itop ;
00338    static int nupold=-1 ;
00339    static float f00[MAX_NUP], fp1[MAX_NUP] ;
00340 
00341    /*-- sanity checks --*/
00342 
00343    if( nup < 1 || nup > MAX_NUP || nar < 2 || far == NULL || fout == NULL ) return ;
00344 
00345    if( nup == 1 ){ memcpy( fout, far, sizeof(float)*nar ); return; }
00346 
00347    /*-- initialize interpolation coefficient, if nup has changed --*/
00348 
00349    if( nup != nupold ){
00350      float val ;
00351      for( kk=0 ; kk < nup ; kk++ ){
00352        val = ((float)kk) / ((float)nup) ;
00353        f00[kk] = 1.0 - val ; fp1[kk] = val ;
00354      }
00355      nupold = nup ;
00356    }
00357 
00358    /*-- interpolate the intermediate places --*/
00359 
00360    ibot = 0 ; itop = nar-2 ;
00361 
00362    switch( nup ){
00363       default:
00364         for( ii=ibot ; ii <= itop ; ii++ )
00365           for( kk=0 ; kk < nup ; kk++ ) fout[kk+ii*nup] = INT1(kk,ii) ;
00366       break ;
00367 
00368       case 2:
00369         for( ii=ibot ; ii <= itop ; ii++ ){
00370           fout[ii*nup]   = INT1(0,ii) ; fout[ii*nup+1] = INT1(1,ii) ;
00371         }
00372       break ;
00373 
00374       case 3:
00375         for( ii=ibot ; ii <= itop ; ii++ ){
00376           fout[ii*nup]   = INT1(0,ii) ; fout[ii*nup+1] = INT1(1,ii) ;
00377           fout[ii*nup+2] = INT1(2,ii) ;
00378         }
00379       break ;
00380 
00381       case 4:
00382         for( ii=ibot ; ii <= itop ; ii++ ){
00383           fout[ii*nup]   = INT1(0,ii) ; fout[ii*nup+1] = INT1(1,ii) ;
00384           fout[ii*nup+2] = INT1(2,ii) ; fout[ii*nup+3] = INT1(3,ii) ;
00385         }
00386       break ;
00387    }
00388 
00389    /*-- interpolate the outside edges --*/
00390 
00391 #if 0                             /* nugatory */
00392    for( ii=0 ; ii < ibot ; ii++ )
00393      for( kk=0 ; kk < nup ; kk++ ) fout[kk+ii*nup] = FINT1(kk,ii) ;
00394 #endif
00395 
00396    for( ii=itop+1 ; ii < nar ; ii++ )
00397      for( kk=0 ; kk < nup ; kk++ ) fout[kk+ii*nup] =  FINT1(kk,ii) ;
00398 
00399    return ;
00400 }
00401 
00402 /*------------------------------------------------------------------*/
00403 /*! Upsample an array of bytes by exactly 2.
00404 --------------------------------------------------------------------*/
00405 
00406 static void upsample_1by2( int nar, byte *bar , byte *bout )
00407 {
00408    int ii ;
00409    if( nar < 1 || bar == NULL || bout == NULL ) return ;
00410 
00411    for( ii=0 ; ii < nar-1 ; ii++ ){
00412       bout[2*ii]   = bar[ii] ;
00413       bout[2*ii+1] = (bar[ii]+bar[ii+1]) >> 1 ;
00414    }
00415    bout[2*nar-1] = bout[2*nar-2] = bar[nar-1] ;
00416 }
00417 
00418 /*------------------------------------------------------------------*/
00419 /*! Upsample an array of bytes by exactly 3.
00420 --------------------------------------------------------------------*/
00421 
00422 static void upsample_1by3( int nar, byte *bar , byte *bout )
00423 {
00424    int ii ;
00425    if( nar < 1 || bar == NULL || bout == NULL ) return ;
00426 
00427    /* Note that 85/256 is about 1/3 and 171/256 is about 2/3;  */
00428    /* by using this trick, we avoid division and so are faster */
00429 
00430    for( ii=0 ; ii < nar-1 ; ii++ ){
00431       bout[3*ii]   = bar[ii] ;
00432       bout[3*ii+1] = (171*bar[ii]+ 85*bar[ii+1]) >> 8 ;
00433       bout[3*ii+2] = ( 85*bar[ii]+171*bar[ii+1]) >> 8 ;
00434    }
00435    bout[3*nar-1] = bout[3*nar-2] = bout[3*nar-3] = bar[nar-1] ;
00436 }
00437 
00438 /*------------------------------------------------------------------*/
00439 /*! Upsample an array of bytes by exactly 4.
00440 --------------------------------------------------------------------*/
00441 
00442 static void upsample_1by4( int nar, byte *bar , byte *bout )
00443 {
00444    int ii ;
00445    if( nar < 1 || bar == NULL || bout == NULL ) return ;
00446 
00447    for( ii=0 ; ii < nar-1 ; ii++ ){
00448       bout[4*ii]   = bar[ii] ;
00449       bout[4*ii+1] = (3*bar[ii]+  bar[ii+1]) >> 2 ;
00450       bout[4*ii+2] = (  bar[ii]+  bar[ii+1]) >> 1 ;
00451       bout[4*ii+3] = (  bar[ii]+3*bar[ii+1]) >> 2 ;
00452    }
00453    bout[4*nar-1] = bout[4*nar-2] = bout[4*nar-3] = bout[4*nar-4] = bar[nar-1];
00454 }
00455 
00456 /*************************************************************************/
00457 /*************************************************************************/
00458 
00459 /* added option for nearest neighbor (if AFNI_IMAGE_ZOOM_NN is Y/y)
00460                                                      22 Feb 2004 [rickr] */
00461 
00462 static MRI_IMAGE * mri_dup2D_rgb_NN( MRI_IMAGE *inim, int nup )
00463 {
00464    rgbyte *bin , *bout , *bin1 ;
00465    MRI_IMAGE *outim ;
00466    int ii,jj,kk,ll , nx,ny , nxup,nyup ;
00467 
00468 ENTRY("mri_dup2D_rgb_NN") ;
00469    if( inim == NULL || inim->kind != MRI_rgb ) RETURN(NULL);
00470 
00471    bin = (rgbyte *) MRI_RGB_PTR(inim); if( bin == NULL ) RETURN(NULL);
00472 
00473    /* make output image **/
00474 
00475    nx = inim->nx ; ny = inim->ny ; nxup = nup*nx ; nyup = nup*ny ;
00476    outim = mri_new( nxup , nyup , MRI_rgb ) ;
00477    bout  = (rgbyte *) MRI_RGB_PTR(outim) ;
00478 
00479    for( jj=0 ; jj < ny ; jj++ ){   /* loop over input rows */
00480      
00481      for ( kk= 0; kk < nup; kk++ ) { /* do rows nup times */
00482 
00483        bin1 = bin;
00484 
00485        for( ii=0 ; ii < nx ; ii++ ){
00486 
00487          for ( ll= 0; ll < nup; ll++ ) {
00488              *bout++ = *bin1;
00489          }
00490 
00491          bin1++;
00492        }
00493 
00494      }
00495      bin += nx ;
00496    }
00497 
00498    MRI_COPY_AUX(outim,inim) ;
00499    RETURN(outim) ;
00500 }
00501 
00502 /*************************************************************************/
00503 /*************************************************************************/
00504 
00505 static MRI_IMAGE * mri_dup2D_rgb4( MRI_IMAGE *inim )
00506 {
00507    rgbyte *bin , *bout , *bin1,*bin2 , *bout1,*bout2,*bout3,*bout4 ;
00508    MRI_IMAGE *outim ;
00509    int ii,jj , nx,ny , nxup,nyup ;
00510 
00511 ENTRY("mri_dup2D_rgb4") ;
00512    if( inim == NULL || inim->kind != MRI_rgb ) RETURN(NULL);
00513 
00514    bin = (rgbyte *) MRI_RGB_PTR(inim); if( bin == NULL ) RETURN(NULL);
00515 
00516    /* make output image **/
00517 
00518    nx = inim->nx ; ny = inim->ny ; nxup = 4*nx ; nyup = 4*ny ;
00519    outim = mri_new( nxup , nyup , MRI_rgb ) ;
00520    bout  = (rgbyte *) MRI_RGB_PTR(outim) ;
00521 
00522    /** macros for the 16 different interpolations
00523        between the four corners:   ul    ur >   00 10 20 30
00524                                             >=> 01 11 21 31
00525                                             >=> 02 12 22 32
00526                                    ll    lr >   03 13 23 33 **/
00527 
00528 #define BOUT(ul,ur,ll,lr,a,b,c,d) ((a*ul+b*ur+c*ll+d*lr) >> 6)
00529 
00530 #define BOUT_00(ul,ur,ll,lr) BOUT(ul,ur,ll,lr,49, 7, 7, 1)
00531 #define BOUT_10(ul,ur,ll,lr) BOUT(ul,ur,ll,lr,35,21, 5, 3)
00532 #define BOUT_20(ul,ur,ll,lr) BOUT(ul,ur,ll,lr,21,35, 3, 5)
00533 #define BOUT_30(ul,ur,ll,lr) BOUT(ul,ur,ll,lr, 7,49, 1, 7)
00534 
00535 #define BOUT_01(ul,ur,ll,lr) BOUT(ul,ur,ll,lr,35, 5,21, 3)
00536 #define BOUT_11(ul,ur,ll,lr) BOUT(ul,ur,ll,lr,25,15,15, 9)
00537 #define BOUT_21(ul,ur,ll,lr) BOUT(ul,ur,ll,lr,15,25, 9,15)
00538 #define BOUT_31(ul,ur,ll,lr) BOUT(ul,ur,ll,lr, 5,35, 3,21)
00539 
00540 #define BOUT_02(ul,ur,ll,lr) BOUT(ul,ur,ll,lr,21, 3,35, 5)
00541 #define BOUT_12(ul,ur,ll,lr) BOUT(ul,ur,ll,lr,15, 9,25,15)
00542 #define BOUT_22(ul,ur,ll,lr) BOUT(ul,ur,ll,lr, 9,15,15,25)
00543 #define BOUT_32(ul,ur,ll,lr) BOUT(ul,ur,ll,lr, 3,21, 5,35)
00544 
00545 #define BOUT_03(ul,ur,ll,lr) BOUT(ul,ur,ll,lr, 7, 1,49, 7)
00546 #define BOUT_13(ul,ur,ll,lr) BOUT(ul,ur,ll,lr, 5, 3,35,21)
00547 #define BOUT_23(ul,ur,ll,lr) BOUT(ul,ur,ll,lr, 3, 5,21,35)
00548 #define BOUT_33(ul,ur,ll,lr) BOUT(ul,ur,ll,lr, 1, 7, 7,49)
00549 
00550 
00551   /** do 16 interpolations between  ul ur
00552                                     ll lr  for index #k, color #c **/
00553 
00554 #define FOUR_ROWS(k,c,ul,ur,ll,lr)                    \
00555    { bout1[4*k+2].c = BOUT_00(ul.c,ur.c,ll.c,lr.c) ;  \
00556      bout1[4*k+3].c = BOUT_10(ul.c,ur.c,ll.c,lr.c) ;  \
00557      bout1[4*k+4].c = BOUT_20(ul.c,ur.c,ll.c,lr.c) ;  \
00558      bout1[4*k+5].c = BOUT_30(ul.c,ur.c,ll.c,lr.c) ;  \
00559      bout2[4*k+2].c = BOUT_01(ul.c,ur.c,ll.c,lr.c) ;  \
00560      bout2[4*k+3].c = BOUT_11(ul.c,ur.c,ll.c,lr.c) ;  \
00561      bout2[4*k+4].c = BOUT_21(ul.c,ur.c,ll.c,lr.c) ;  \
00562      bout2[4*k+5].c = BOUT_31(ul.c,ur.c,ll.c,lr.c) ;  \
00563      bout3[4*k+2].c = BOUT_02(ul.c,ur.c,ll.c,lr.c) ;  \
00564      bout3[4*k+3].c = BOUT_12(ul.c,ur.c,ll.c,lr.c) ;  \
00565      bout3[4*k+4].c = BOUT_22(ul.c,ur.c,ll.c,lr.c) ;  \
00566      bout3[4*k+5].c = BOUT_32(ul.c,ur.c,ll.c,lr.c) ;  \
00567      bout4[4*k+2].c = BOUT_03(ul.c,ur.c,ll.c,lr.c) ;  \
00568      bout4[4*k+3].c = BOUT_13(ul.c,ur.c,ll.c,lr.c) ;  \
00569      bout4[4*k+4].c = BOUT_23(ul.c,ur.c,ll.c,lr.c) ;  \
00570      bout4[4*k+5].c = BOUT_33(ul.c,ur.c,ll.c,lr.c) ; }
00571 
00572   /** do the above for all 3 colors */
00573 
00574 #define FOUR_RGB(k,ul,ur,ll,lr)  { FOUR_ROWS(k,r,ul,ur,ll,lr) ; \
00575                                    FOUR_ROWS(k,g,ul,ur,ll,lr) ; \
00576                                    FOUR_ROWS(k,b,ul,ur,ll,lr) ;  }
00577 
00578    bin1  = bin    ;                  /* 2 input rows */
00579    bin2  = bin+nx ;
00580 
00581    bout1 = bout+2*nxup ;             /* 4 output rows */
00582    bout2 = bout1+nxup  ;
00583    bout3 = bout2+nxup  ;
00584    bout4 = bout3+nxup  ;
00585 
00586    for( jj=0 ; jj < ny-1 ; jj++ ){   /* loop over input rows */
00587 
00588      for( ii=0 ; ii < nx-1 ; ii++ ){
00589         FOUR_RGB(ii,bin1[ii],bin1[ii+1],bin2[ii],bin2[ii+1]) ;
00590      }
00591 
00592      /* at this point, output rows have elements [2..4*nx-3] filled;
00593         now copy [2] into [0..1] and [4*nx-3] into [4*nx-2..4*nx-1] */
00594 
00595      bout1[0] = bout1[1] = bout1[2] ;
00596      bout2[0] = bout2[1] = bout2[2] ;
00597      bout3[0] = bout3[1] = bout3[2] ;
00598      bout4[0] = bout4[1] = bout4[2] ;
00599 
00600      bout1[nxup-2] = bout1[nxup-1] = bout1[nxup-2] ;
00601      bout2[nxup-2] = bout2[nxup-1] = bout2[nxup-2] ;
00602      bout3[nxup-2] = bout3[nxup-1] = bout3[nxup-2] ;
00603      bout4[nxup-2] = bout4[nxup-1] = bout4[nxup-2] ;
00604 
00605      /* advance input and output rows */
00606 
00607      bin1 = bin2; bin2 += nx ;
00608      bout1 = bout4+nxup; bout2 = bout1+nxup;
00609      bout3 = bout2+nxup; bout4 = bout3+nxup;
00610    }
00611 
00612    /* copy row 2 into rows 0 and row 1 */
00613 
00614    bout1 = bout ; bout2 = bout1+nxup ; bout3 = bout2+nxup ;
00615    for( ii=0 ; ii < nxup ; ii++ )
00616       bout1[ii] = bout2[ii] = bout3[ii] ;
00617 
00618    /* copy rown nyup-3 into rows nyup-2 and nyup-1 */
00619 
00620    bout1 = bout + (nyup-3)*nxup ; bout2 = bout1+nxup ; bout3 = bout2+nxup ;
00621    for( ii=0 ; ii < nxup ; ii++ )
00622       bout2[ii] = bout3[ii] = bout1[ii] ;
00623 
00624    MRI_COPY_AUX(outim,inim) ;
00625    RETURN(outim) ;
00626 }
00627 
00628 /*************************************************************************/
00629 
00630 static MRI_IMAGE * mri_dup2D_rgb3( MRI_IMAGE *inim )
00631 {
00632    rgbyte *bin , *bout , *bin1,*bin2 , *bout1,*bout2,*bout3 ;
00633    MRI_IMAGE *outim ;
00634    int ii,jj , nx,ny , nxup,nyup ;
00635 
00636 ENTRY("mri_dup2D_rgb3") ;
00637    if( inim == NULL || inim->kind != MRI_rgb ) RETURN(NULL) ;
00638 
00639    bin = (rgbyte *) MRI_RGB_PTR(inim); if( bin == NULL ) RETURN(NULL);
00640 
00641    /* make output image **/
00642 
00643    nx = inim->nx ; ny = inim->ny ; nxup = 3*nx ; nyup = 3*ny ;
00644    outim = mri_new( nxup , nyup , MRI_rgb ) ;
00645    bout  = (rgbyte *) MRI_RGB_PTR(outim) ;
00646 
00647    /** macros for the 9 different interpolations
00648        between the four corners:   ul  ur >   00 10 20
00649                                           >=> 01 11 21
00650                                    ll  lr >   02 12 22 **/
00651 
00652 #define COUT_00(ul,ur,ll,lr) (     ul                              )
00653 #define COUT_10(ul,ur,ll,lr) ((171*ul + 85*ur                ) >> 8)
00654 #define COUT_20(ul,ur,ll,lr) (( 85*ul +171*ur                ) >> 8)
00655 
00656 #define COUT_01(ul,ur,ll,lr) ((171*ul +         85*ll        ) >> 8)
00657 #define COUT_11(ul,ur,ll,lr) ((114*ul + 57*ur + 57*ll + 28*lr) >> 8)
00658 #define COUT_21(ul,ur,ll,lr) (( 57*ul +114*ur + 28*ll + 57*lr) >> 8)
00659 
00660 #define COUT_02(ul,ur,ll,lr) (( 85*ul +        171*ll        ) >> 8)
00661 #define COUT_12(ul,ur,ll,lr) (( 57*ul + 28*ur +114*ll + 57*lr) >> 8)
00662 #define COUT_22(ul,ur,ll,lr) (( 28*ul + 57*ur + 57*ll +114*lr) >> 8)
00663 
00664   /** do 9 interpolations between  ul ur
00665                                    ll lr  for index #k, color #c **/
00666 
00667 #define THREE_ROWS(k,c,ul,ur,ll,lr)                  \
00668    { bout1[3*k+1].c = COUT_00(ul.c,ur.c,ll.c,lr.c) ; \
00669      bout1[3*k+2].c = COUT_10(ul.c,ur.c,ll.c,lr.c) ; \
00670      bout1[3*k+3].c = COUT_20(ul.c,ur.c,ll.c,lr.c) ; \
00671      bout2[3*k+1].c = COUT_01(ul.c,ur.c,ll.c,lr.c) ; \
00672      bout2[3*k+2].c = COUT_11(ul.c,ur.c,ll.c,lr.c) ; \
00673      bout2[3*k+3].c = COUT_21(ul.c,ur.c,ll.c,lr.c) ; \
00674      bout3[3*k+1].c = COUT_02(ul.c,ur.c,ll.c,lr.c) ; \
00675      bout3[3*k+2].c = COUT_12(ul.c,ur.c,ll.c,lr.c) ; \
00676      bout3[3*k+3].c = COUT_22(ul.c,ur.c,ll.c,lr.c) ;  }
00677 
00678   /** do the above for all 3 colors **/
00679 
00680 #define THREE_RGB(k,ul,ur,ll,lr)  { THREE_ROWS(k,r,ul,ur,ll,lr) ; \
00681                                     THREE_ROWS(k,g,ul,ur,ll,lr) ; \
00682                                     THREE_ROWS(k,b,ul,ur,ll,lr) ;  }
00683 
00684    bin1  = bin       ; bin2  = bin+nx    ;  /* 2 input rows */
00685    bout1 = bout +nxup; bout2 = bout1+nxup;  /* 3 output rows */
00686    bout3 = bout2+nxup;
00687 
00688    for( jj=0 ; jj < ny-1 ; jj++ ){   /* loop over input rows */
00689 
00690      for( ii=0 ; ii < nx-1 ; ii++ ){
00691         THREE_RGB(ii,bin1[ii],bin1[ii+1],bin2[ii],bin2[ii+1]) ;
00692      }
00693 
00694      /* here, ii=nx-1; can't use ii+1 */
00695      /* do the 3*ii+1 output only,    */
00696      /* and copy into 3*ii+2 column   */
00697 
00698      bout1[3*ii+1].r = COUT_00(bin1[ii].r,0,bin2[ii].r,0) ;
00699      bout2[3*ii+1].r = COUT_01(bin1[ii].r,0,bin2[ii].r,0) ;
00700      bout3[3*ii+1].r = COUT_02(bin1[ii].r,0,bin2[ii].r,0) ;
00701 
00702      bout1[3*ii+1].g = COUT_00(bin1[ii].g,0,bin2[ii].g,0) ;
00703      bout2[3*ii+1].g = COUT_01(bin1[ii].g,0,bin2[ii].g,0) ;
00704      bout3[3*ii+1].g = COUT_02(bin1[ii].g,0,bin2[ii].g,0) ;
00705 
00706      bout1[3*ii+1].b = COUT_00(bin1[ii].b,0,bin2[ii].b,0) ;
00707      bout2[3*ii+1].b = COUT_01(bin1[ii].b,0,bin2[ii].b,0) ;
00708      bout3[3*ii+1].b = COUT_02(bin1[ii].b,0,bin2[ii].b,0) ;
00709 
00710      bout1[3*ii+2] = bout1[3*ii+1] ;
00711      bout2[3*ii+2] = bout2[3*ii+1] ;
00712      bout3[3*ii+2] = bout3[3*ii+1] ;
00713 
00714      /* also copy [0] column output from column [1] */
00715 
00716      bout1[0] = bout1[1] ; bout2[0] = bout2[1] ; bout3[0] = bout3[1] ;
00717 
00718      /* advance input and output rows */
00719 
00720      bin1 = bin2; bin2 += nx ;
00721      bout1 = bout3+nxup; bout2 = bout1+nxup;
00722      bout3 = bout2+nxup;
00723    }
00724 
00725    /* here, jj=ny-1, so can't use jj+1 (bin2) */
00726    /* do the bout1 output row only, copy into bout2 */
00727 
00728    for( ii=0 ; ii < nx-1 ; ii++ ){
00729      bout1[3*ii+1].r = COUT_00(bin1[ii].r,bin1[ii+1].r,0,0) ;
00730      bout1[3*ii+2].r = COUT_10(bin1[ii].r,bin1[ii+1].r,0,0) ;
00731      bout1[3*ii+3].r = COUT_20(bin1[ii].r,bin1[ii+1].r,0,0) ;
00732 
00733      bout1[3*ii+1].g = COUT_00(bin1[ii].g,bin1[ii+1].g,0,0) ;
00734      bout1[3*ii+2].g = COUT_10(bin1[ii].g,bin1[ii+1].g,0,0) ;
00735      bout1[3*ii+3].g = COUT_20(bin1[ii].g,bin1[ii+1].g,0,0) ;
00736 
00737      bout1[3*ii+1].b = COUT_00(bin1[ii].b,bin1[ii+1].b,0,0) ;
00738      bout1[3*ii+2].b = COUT_10(bin1[ii].b,bin1[ii+1].b,0,0) ;
00739      bout1[3*ii+3].b = COUT_20(bin1[ii].b,bin1[ii+1].b,0,0) ;
00740 
00741      bout2[3*ii+1] = bout1[3*ii+1] ;
00742      bout2[3*ii+2] = bout1[3*ii+2] ;
00743      bout2[3*ii+3] = bout1[3*ii+3] ;
00744    }
00745 
00746    /* do bottom corners */
00747 
00748    bout1[0] = bout2[0] = bout2[1] = bout1[1] ;
00749    bout1[nxup-2] = bout1[nxup-1] = bout2[nxup-2] = bout2[nxup-1] = bout1[nxup-3] ;
00750 
00751    /* do first row of output as well */
00752 
00753    bout3 = bout+nxup ;
00754    for( ii=0 ; ii < nxup ; ii++ ) bout[ii] = bout3[ii] ;
00755 
00756    MRI_COPY_AUX(outim,inim) ;
00757    RETURN(outim) ;
00758 }
00759 
00760 /*************************************************************************/
00761 /*************************************************************************/
00762 
00763 static MRI_IMAGE * mri_dup2D_rgb2( MRI_IMAGE *inim )
00764 {
00765    rgbyte *bin , *bout , *bin1,*bin2 , *bout1,*bout2 ;
00766    MRI_IMAGE *outim ;
00767    int ii,jj , nx,ny , nxup,nyup ;
00768 
00769 ENTRY("mri_dup2D_rgb2") ;
00770    if( inim == NULL || inim->kind != MRI_rgb ) RETURN(NULL) ;
00771 
00772    bin = (rgbyte *) MRI_RGB_PTR(inim); if( bin == NULL ) RETURN(NULL);
00773 
00774    /* make output image **/
00775 
00776    nx = inim->nx ; ny = inim->ny ; nxup = 2*nx ; nyup = 2*ny ;
00777    outim = mri_new( nxup , nyup , MRI_rgb ) ;
00778    bout  = (rgbyte *) MRI_RGB_PTR(outim) ;
00779 
00780    /** macros for the 4 different interpolations
00781        between the four corners:   ul  ur >=> 00 10
00782                                    ll  lr >=> 01 11  **/
00783 
00784 #define DOUT(ul,ur,ll,lr,a,b,c,d) ((a*ul+b*ur+c*ll+d*lr) >> 4)
00785 
00786 #define DOUT_00(ul,ur,ll,lr) DOUT(ul,ur,ll,lr,9,3,3,1)
00787 #define DOUT_10(ul,ur,ll,lr) DOUT(ul,ur,ll,lr,3,9,1,3)
00788 
00789 #define DOUT_01(ul,ur,ll,lr) DOUT(ul,ur,ll,lr,3,1,9,3)
00790 #define DOUT_11(ul,ur,ll,lr) DOUT(ul,ur,ll,lr,1,3,3,9)
00791 
00792   /** do r interpolations between  ul ur
00793                                    ll lr  for index #k, color #c **/
00794 
00795 #define TWO_ROWS(k,c,ul,ur,ll,lr)                     \
00796    { bout1[2*k+1].c = DOUT_00(ul.c,ur.c,ll.c,lr.c) ;  \
00797      bout1[2*k+2].c = DOUT_10(ul.c,ur.c,ll.c,lr.c) ;  \
00798      bout2[2*k+1].c = DOUT_01(ul.c,ur.c,ll.c,lr.c) ;  \
00799      bout2[2*k+2].c = DOUT_11(ul.c,ur.c,ll.c,lr.c) ; }
00800 
00801   /** do the above for all 3 colors */
00802 
00803 #define TWO_RGB(k,ul,ur,ll,lr)  { TWO_ROWS(k,r,ul,ur,ll,lr) ;  \
00804                                   TWO_ROWS(k,g,ul,ur,ll,lr) ;  \
00805                                   TWO_ROWS(k,b,ul,ur,ll,lr) ; }
00806 
00807    bin1  = bin    ;                 /* 2 input rows */
00808    bin2  = bin+nx ;
00809 
00810    bout1 = bout +nxup ;             /* 2 output rows */
00811    bout2 = bout1+nxup  ;
00812 
00813    for( jj=0 ; jj < ny-1 ; jj++ ){   /* loop over input rows */
00814 
00815      for( ii=0 ; ii < nx-1 ; ii++ ){
00816         TWO_RGB(ii,bin1[ii],bin1[ii+1],bin2[ii],bin2[ii+1]) ;
00817      }
00818 
00819      /* at this point, output rows have elements [1..2*nx-2] filled;
00820         now copy [1] into [0] and [2*nx-2] into [2*nx-1]            */
00821 
00822      bout1[0] = bout1[1] ;
00823      bout2[0] = bout2[1] ;
00824 
00825      bout1[nxup-1] = bout1[nxup-2] ;
00826      bout2[nxup-1] = bout2[nxup-2] ;
00827 
00828      /* advance input and output rows */
00829 
00830      bin1 = bin2; bin2 += nx ;
00831      bout1 = bout2+nxup; bout2 = bout1+nxup;
00832    }
00833 
00834    /* copy row 1 into row 0 */
00835 
00836    bout1 = bout ; bout2 = bout1+nxup ;
00837    for( ii=0 ; ii < nxup ; ii++ )
00838       bout1[ii] = bout2[ii] ;
00839 
00840    /* copy rown nyup-2 into row nyup-1 */
00841 
00842    bout1 = bout + (nyup-2)*nxup ; bout2 = bout1+nxup ;
00843    for( ii=0 ; ii < nxup ; ii++ )
00844       bout2[ii] = bout1[ii] ;
00845 
00846    MRI_COPY_AUX(outim,inim) ;
00847    RETURN(outim) ;
00848 }