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 #include <string.h>
00009 
00010 void AB_interp( int , double ) ;
00011 
00012 #define LINEAR 7
00013 #define BLOCKY 8
00014 
00015 void Syntax(void)
00016 {
00017    printf(
00018      "ABUT:  put noncontiguous FMRI slices together [for to3d]\n"
00019      "\n"
00020      "method: put zero valued slices in the gaps, then\n"
00021      "        replicate images to simulate thinner slices\n"
00022      "\n"
00023      "Usage:\n"
00024      "   abut [-dzin thickness] [-dzout thickness] [-root name]\n"
00025      "        [-linear | -blocky] [-verbose] [-skip n+gap] ... images ...\n"
00026      "\n"
00027      "   -dzin   the thickness value in mm;  if not given,\n"
00028      "             taken to be 1.0 (in which case, the output\n"
00029      "             thickness and gap sizes are simply relative\n"
00030      "             to the slice thickness, rather than absolute)\n"
00031      "\n"
00032      "   -dzout  the output slice thickness, usually smaller than\n"
00033      "             the input thickness;  if not given, the program\n"
00034      "             will compute a value (the smaller the ratio\n"
00035      "             dzout/dzin is, the more slices will be output)\n"
00036      "\n"
00037      "   -root   'name' is the root (or prefix) for the output\n"
00038      "             slice filename;  for example, '-root fred.'\n"
00039      "             will result in files fred.0001, fred.0002, ...\n"
00040      "\n"
00041      "   -linear if present, this flag indicates that subdivided slices\n"
00042      "             will be linearly interpolated rather than simply\n"
00043      "             replicated -- this will make the results smoother\n"
00044      "             in the through-slice direction (if dzout < dzin)\n"
00045      "\n"
00046      "   -blocky similar to -linear, but uses AFNI's 'blocky' interpolation\n"
00047      "             when possible to put out intermediate slices.\n"
00048      "             Both interpolation options only apply when dzout < dzin\n"
00049      "             and when an output slice has a non-gappy neighbor.\n"
00050      "\n"
00051      "   -skip   'n+gap' indicates that a gap is to be inserted\n"
00052      "             between input slices #n and #n+1, where n=1,2,...;\n"
00053      "             for example, -skip 6+5.5 means put a gap of 5.5 mm\n"
00054      "             between slices 6 and 7.\n"
00055      "\n"
00056      "   More than one -skip option is allowed.  They must all occur\n"
00057      "   before the list of input image filenames.\n"
00058    ) ;
00059 
00060    exit(0) ;
00061 }
00062 
00063 /*-------------------------------------------------------------------*/
00064 
00065 #define SLICES_MAX 128
00066 
00067 static MRI_IMAGE * imin[SLICES_MAX] , * zero_im ;
00068 static double gap[SLICES_MAX] ;
00069 static double zbot[SLICES_MAX] , ztop[SLICES_MAX] ;
00070 
00071 static double dzin = 1.0 , dzout = 0.0 , gap_max , gap_sum ;
00072 static char prefix[256] = "abut." ;
00073 static int interp = 0 ;
00074 static int verbose = 0 ;
00075 static int iout , nfiles , nx , ny ;
00076 static char fname[256] ;
00077 
00078 int main( int argc , char * argv[] )
00079 {
00080    int isl , nopt , ii ;
00081    double fff , zout ;
00082    short * zar ;
00083 
00084    /*----- initialize -----*/
00085 
00086    if( argc < 3 || strncmp(argv[1],"-help",2) == 0 ) Syntax() ;
00087 
00088    machdep() ;
00089 
00090    for( isl=0 ; isl < SLICES_MAX ; isl++ ){
00091       imin[isl] = NULL ;
00092       gap[isl]  = 0.0 ;
00093    }
00094 
00095    /*----- read command line options -----*/
00096 
00097    nopt = 1 ;
00098 
00099    while( nopt < argc && argv[nopt][0] == '-' ){
00100 
00101       /*--- -verbose ---*/
00102 
00103       if( strncmp(argv[nopt],"-verbose",4) == 0 ){
00104          verbose = 1 ;
00105          nopt++ ; continue ;
00106       }
00107 
00108       /*--- -linear ---*/
00109 
00110       if( strncmp(argv[nopt],"-linear",4) == 0 ){
00111          interp = LINEAR ;
00112          nopt++ ; continue ;
00113       }
00114 
00115       /*--- -blocky ---*/
00116 
00117       if( strncmp(argv[nopt],"-blocky",4) == 0 ){
00118          interp = BLOCKY ;
00119          nopt++ ; continue ;
00120       }
00121 
00122       /*--- -dzin thickness ---*/
00123 
00124       if( strncmp(argv[nopt],"-dzin",4) == 0 ){
00125          if( ++nopt >= argc ){
00126             fprintf(stderr,"\n*** no argument for -dzin?\n\n") ;
00127             exit(1) ;
00128          }
00129          fff = strtod( argv[nopt] , NULL ) ;
00130          if( fff <= 0.0 ){
00131             fprintf(stderr,"\n*** illegal argument for -dzin: %f\n\n",fff) ;
00132             exit(1) ;
00133          }
00134          dzin = fff ;
00135          nopt++ ; continue ;
00136       }
00137 
00138       /*--- -dzout thickness ---*/
00139 
00140       if( strncmp(argv[nopt],"-dzout",4) == 0 ){
00141          if( ++nopt >= argc ){
00142             fprintf(stderr,"\n*** no argument for -dzout?\n\n") ;
00143             exit(1) ;
00144          }
00145          fff = strtod( argv[nopt] , NULL ) ;
00146          if( fff <= 0.0 ){
00147             fprintf(stderr,"\n*** illegal argument for -dzout: %f\n\n",fff) ;
00148             exit(1) ;
00149          }
00150          dzout = fff ;
00151          nopt++ ; continue ;
00152       }
00153 
00154       /*--- -root name ---*/
00155 
00156       if( strncmp(argv[nopt],"-root",4) == 0 ){
00157          if( ++nopt >= argc ){
00158             fprintf(stderr,"\n*** no argument for -root?\n\n") ;
00159             exit(1) ;
00160          }
00161          strcpy( prefix , argv[nopt] ) ;
00162          ii = strlen(prefix) ;
00163          if( prefix[ii-1] != '.' ){  /* add a period? */
00164             prefix[ii]   = '.' ;
00165             prefix[ii+1] = '\0' ;
00166          }
00167          nopt++ ; continue ;
00168       }
00169 
00170       /*--- -skip n+gap ---*/
00171 
00172       if( strncmp(argv[nopt],"-skip",4) == 0 ){
00173          int nok , nnn = -1 ;
00174          double ggg = -1.0 ;
00175          if( ++nopt >= argc ){
00176             fprintf(stderr,"\n*** no argument for -skip?\n\n") ;
00177             exit(1) ;
00178          }
00179 
00180          nok = sscanf( argv[nopt] , "%d+%lf" , &nnn , &ggg ) ;
00181 
00182          if( nok != 2 || nnn <= 0 || nnn > SLICES_MAX || ggg < 0 ){
00183             fprintf(stderr,
00184                     "\n*** illegal argument for -skip: %s\n\n",argv[nopt]) ;
00185             exit(1) ;
00186          }
00187 
00188          gap[nnn-1] = ggg ;    /* external slice indices are 1 more */
00189          nopt++ ; continue ;   /* than the internal index system!   */
00190       }
00191    }
00192 
00193    /*----- read input images -----*/
00194 
00195    if( nopt >= argc ){
00196       fprintf(stderr,"\n*** no input image files!!!\n\n") ;
00197       exit(1) ;
00198    }
00199 
00200    nfiles = argc - nopt ;
00201    if( nfiles > SLICES_MAX ){
00202       fprintf(stderr,
00203               "\n*** You input %d files, but max allowed is %d!!!\n\n" ,
00204               nfiles , SLICES_MAX ) ;
00205       exit(1) ;
00206    } else if( nfiles < 1 ){
00207       fprintf(stderr,"\n*** no input image files!!!\n\n") ;
00208       exit(1) ;
00209    } else if( nfiles == 1 ){
00210       fprintf(stderr,"\n*** only one input image file!!!\n\n") ;
00211       exit(1) ;
00212    }
00213 
00214    imin[0] = mri_read_just_one( argv[nopt] ) ;
00215    if( imin[0] == NULL ){
00216       fprintf(stderr,"\n*** Cannot read 1st image from file %s\n",argv[nopt]);
00217       exit(1) ;
00218    }
00219    if( imin[0]->kind != MRI_short ){
00220       fprintf(stderr,"\n*** Cannot deal with non-short images!\n") ;
00221       exit(1) ;
00222    }
00223    nx = imin[0]->nx ; ny = imin[0]->ny ;
00224 
00225    for( isl=1 ; isl < nfiles ; isl++ ){
00226       imin[isl] = mri_read_just_one( argv[nopt+isl] ) ;
00227       if( imin[isl] == NULL ){
00228          fprintf(stderr,"\n*** Cannot read 1st image from file %s\n",argv[nopt+isl]) ;
00229          exit(1) ;
00230       }
00231 
00232       if( imin[isl]->nx != nx || imin[isl]->ny != ny ){
00233          fprintf(stderr,
00234                  "\n*** first file is %d x %d but file %s is %d x %d!\n\n",
00235                  nx,ny , argv[nopt+isl] , imin[isl]->nx , imin[isl]->ny ) ;
00236          exit(1) ;
00237       }
00238 
00239       if( imin[isl]->kind != MRI_short ){
00240          fprintf(stderr,
00241                  "\n*** image file %s is not an image of shorts!\n",
00242                  argv[nopt+isl] ) ;
00243          exit(1) ;
00244       }
00245    }  /* all images input when this loop exits */
00246 
00247    zero_im = mri_new( nx , ny , MRI_short ) ;
00248    zar     = mri_data_pointer( zero_im ) ;
00249    for( ii=0 ; ii < (nx*ny) ; ii++ ) zar[ii] = 0 ;
00250 
00251    /*----- adjust gaps to be relative to dzin -----*/
00252 
00253    gap_max = 0 ;
00254    for( isl=0 ; isl < nfiles-1 ; isl++ ){
00255       gap[isl] = gap[isl] / dzin ;
00256       gap_max  = (gap_max < gap[isl]) ? (gap[isl]) : (gap_max) ;
00257    }
00258 
00259    /*----- compute dzout (relative to dzin) if not already given -----*/
00260 
00261    if( dzout > 0.0 ){
00262       dzout = dzout / dzin ;
00263    } else if( gap_max == 0.0 ){
00264       dzout = 1.0 ;
00265    } else {
00266 #define NBASE 5
00267       int ibase , ibest ;
00268       double fit , worst_fit[NBASE+1];
00269 
00270       for( ibase=1 ; ibase <= NBASE ; ibase++ ){   /* find worst error */
00271          worst_fit[ibase] = 0.0 ;                  /* for dzout = 1/ibase */
00272 
00273          for( isl=0 ; isl < nfiles-1 ; isl++ ){
00274             if( gap[isl] > 0.0 ){
00275                fit = gap[isl] * ibase ;
00276                fit = fabs( fit - floor(fit+0.5) ) ;
00277                if( fit > worst_fit[ibase] ) worst_fit[ibase] = fit ;
00278             }
00279          }
00280       }
00281 
00282       for( ibase=1 ; ibase <= NBASE ; ibase++ )
00283          if( worst_fit[ibase] < 0.05 ) break ;   /* first one < 5% */
00284 
00285       ibest = ibase ;
00286       if( ibest > NBASE ){                       /* otherwise take best */
00287          fit   = worst_fit[1] ;
00288          ibest = 1 ;
00289          for( ibase=2 ; ibase <= NBASE ; ibase++ )
00290             if( worst_fit[ibase] < fit ){
00291                ibest = ibase ;
00292                fit   = worst_fit[ibase] ;
00293             }
00294       }
00295 
00296       dzout = ((double) 1.0) / (double) ibest ;
00297 
00298       printf("-- computed dzout = %f\n" , dzout * dzin ) ;
00299    }
00300 
00301    /*----- adjust gaps to be integer multiples of dzout -----*/
00302 
00303    for( isl=0 ; isl < nfiles-1 ; isl++ ){
00304 
00305       fff      = gap[isl] ;
00306       iout     = (int)( fff / dzout + 0.5 ) ;
00307       gap[isl] = iout * dzout ;
00308 
00309       if( fabs(gap[isl]-fff) > 0.001 )
00310          printf("-- adjusted gap after slice %d to %f\n",
00311                 isl+1 , gap[isl]*dzin ) ;
00312    }
00313 
00314    /*----- loop through output points and output stuff -----*/
00315 
00316    isl  = 0 ;
00317    zout = 0 ;
00318 
00319    gap_sum = 0.0 ;
00320    for( isl=0 ; isl < nfiles ; isl++ ){
00321       zbot[isl] = zout ;
00322       ztop[isl] = zout + 1.0 ;
00323       zout     += 1.0 + gap[isl] ;  /* start of next slice */
00324    }
00325 
00326    zout = 0.001 ;
00327    iout = 0 ;
00328    isl  = 0 ;
00329 
00330    do {
00331 
00332      iout++ ;
00333      sprintf(fname,"%s%04d",prefix,iout) ;
00334 
00335      if( zout >= zbot[isl] && zout < ztop[isl] ){  /* inside this slice */
00336         AB_interp( isl , zout ) ;
00337         continue ;
00338      }
00339 
00340      if( zout >= ztop[isl] ){  /* move to next slice */
00341         isl++ ;
00342      }
00343 
00344      if( zout < zbot[isl] ){  /* before slice ==> in a gap */
00345         mri_write( fname , zero_im ) ;
00346         if(verbose) printf("  -- new slice %d is all zeros\n",iout) ;
00347         continue ;
00348      }
00349 
00350      AB_interp( isl , zout ) ;  /* into the next slice now */
00351 
00352    } while ( (zout += dzout) < ztop[nfiles-1] ) ;
00353 
00354    /*----- DONE -----*/
00355 
00356    printf("-- wrote %d output slices\n",iout) ;
00357    exit(0) ;
00358 }
00359 
00360 void AB_interp( int isl , double zout )
00361 {
00362    int nb , ii ;
00363    double wt , fnow,fnb;
00364    short * snow , * snb , * sint ;
00365    MRI_IMAGE * imint ;
00366 
00367    if( interp != LINEAR && interp != BLOCKY ){    /* this is easy! */
00368       mri_write( fname , imin[isl] ) ;
00369       if(verbose)
00370         printf("  -- new slice %d is a replica of input %d\n",iout,isl+1) ;
00371       return ;
00372    }
00373 
00374    /* must deal with possibility of interpolation */
00375 
00376    /* possible to interpolate below? */
00377 
00378    if( zout-zbot[isl] <= 0.5 &&               /* closer to below */
00379        isl > 0 &&                             /* there is a below */
00380        fabs(ztop[isl-1]-zbot[isl]) < 0.001 ){ /* below is adjacent */
00381 
00382       nb = isl-1 ;                 /* neighbor index */
00383       wt = zout-zbot[isl] + 0.5 ;  /* distance from below mid to zout */
00384 
00385    /* possible to interpolate above? */
00386 
00387    } else if( ztop[isl]-zout <= 0.5 &&               /* closer to above */
00388               isl < nfiles-1 &&                      /* there is an above */
00389               fabs(zbot[isl+1]-ztop[isl]) < 0.001 ){ /* above is adjacent */
00390 
00391       nb = isl+1 ;                  /* neighbor index */
00392       wt = ztop[isl]-zout + 0.5 ;   /* distance from above mid to zout */
00393    } else {
00394       mri_write( fname , imin[isl] ) ;  /* no interp possible */
00395       if(verbose)
00396         printf("  -- new slice %d is a replica of input %d\n",iout,isl+1) ;
00397       return ;
00398    }
00399 
00400    snow  = mri_data_pointer( imin[isl] ) ;
00401    snb   = mri_data_pointer( imin[nb]  ) ;
00402    imint = mri_new( nx , ny , MRI_short ) ;
00403    sint  = mri_data_pointer( imint ) ;
00404 
00405    if( interp == LINEAR ) fnow = wt ;  /* linear */
00406    else                   fnow = 1.0 - 8.0 * pow(1.0-wt,4.0) ;
00407    fnb = 1.0 - fnow ;
00408 
00409    for( ii=0 ; ii < (nx*ny) ; ii++ )
00410       sint[ii] = fnow * snow[ii] + fnb * snb[ii] ;
00411 
00412    if(verbose)
00413      printf(
00414        "  -- new slice %d interpolated from input %d(%6.3f) and %d(%6.3f)\n",
00415        iout , isl+1,fnow,nb+1,fnb ) ;
00416 
00417    mri_write(fname,imint) ;
00418    mri_free(imint) ;
00419    return ;
00420 }