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 /*******************************************************
00008  * 3dFourier                                           *
00009  * T. Ross  and K. Heimerl 7/99                        *
00010  *******************************************************/
00011 
00012 #include "mrilib.h"
00013 #include <stdio.h>
00014 #include <stdlib.h>
00015 #include <string.h>
00016 #include <math.h>
00017 
00018 #if !defined(TRUE)
00019 #define TRUE (1==1)
00020 #define FALSE !TRUE
00021 #endif
00022 
00023 static void *My_Malloc(size_t);
00024 
00025 void Error_Exit(char *message)
00026 {
00027         fprintf(stderr, "\n\nError in 3dFourier:\n%s\n\nTry 3dFourier -help\n",message);
00028         exit(1);
00029 }
00030 
00031 
00032 /* include the filter and filter driver */
00033 #include "fourier_filter.c"
00034 
00035 /***********************************************************************************/
00036 
00037 void  help_message(void)
00038 {
00039         printf(
00040                 "3dFourier \n"
00041                 "(c) 1999 Medical College of Wisconsin\n"
00042                 "by T. Ross and K. Heimerl\n"
00043                 "Version 0.8 last modified 8-17-99\n\n"
00044                 "Usage: 3dFourier [options] dataset\n\n"
00045                 "The paramters and options are:\n"
00046                 "       dataset         an afni compatible 3d+time dataset to be operated upon\n"
00047                 "       -prefix name    output name for new 3d+time dataset [default = fourier]\n"
00048                 "       -lowpass f      low pass filter with a cutoff of f Hz\n"
00049                 "       -highpass f     high pass filter with a cutoff of f Hz\n"
00050                 "       -ignore n       ignore the first n images [default = 1]\n"
00051 /*              "       -autocorr       compute the autocorrelation of the data\n" */
00052                 "       -retrend        Any mean and linear trend are removed before filtering.\n"
00053                 "                       This will restore the trend after filtering.\n"
00054                 "\nNote that by combining the lowpass and highpass options, one can construct\n"
00055                 "bandpass and notch filters\n"
00056         );
00057 
00058         printf("\n" MASTER_SHORTHELP_STRING ) ;
00059 
00060 }
00061 
00062 /*--------------------------------------------------------------------------*/
00063 
00064 int main (int argc, char *argv[])
00065 {
00066         int j, narg=1, autocorr=FALSE, retrend=FALSE, ignore=1;
00067         char *prefix=NULL, *err;
00068         float low_fc = 0.0, high_fc = 0.0;
00069         THD_3dim_dataset *input=NULL;
00070 
00071         if (argc < 2 ||  strcmp(argv[1],"-help") == 0 ) {
00072                 help_message();
00073                 exit(1);
00074         }
00075 
00076         /* Loop over arguements and pull out what we need */
00077         while( narg < argc && argv[narg][0] == '-' ){
00078 
00079                 if( strncmp(argv[narg],"-prefix",5) == 0 ) {
00080                         narg++;
00081                         if (narg==argc)
00082                                 Error_Exit("-prefix must be followed by a file name");
00083                         j = strlen(argv[narg]);
00084                         prefix =(char *)My_Malloc((j+1)*sizeof(char));
00085                         strcpy(prefix, argv[narg++]);
00086                         continue;
00087                 }
00088 
00089                 if (strncmp(argv[narg], "-lowpass", 5) == 0) {
00090                         narg++;
00091                         if (narg==argc)
00092                                 Error_Exit("-lowpass must be followed by a frequency");
00093                         low_fc = (float)atof(argv[narg++]);
00094                         continue;
00095                 }
00096 
00097                 if (strncmp(argv[narg], "-highpass", 5) == 0) {
00098                         narg++;
00099                         if (narg==argc)
00100                                 Error_Exit("-highpass must be followed by a frequency");
00101                         high_fc = (float)atof(argv[narg++]);
00102                         continue;
00103                 }
00104 
00105                 if (strncmp(argv[narg], "-ignore", 5) == 0) {
00106                         narg++;
00107                         if (narg==argc)
00108                                 Error_Exit("-ignore must be followed by an integer");
00109                         ignore = (int)atol(argv[narg++]);
00110                         if ((ignore<0) || (ignore>10))
00111                                 Error_Exit("-ignore must be between 0 and 10, inclusive");
00112                         continue;
00113                 }
00114 /*
00115                 if (strncmp(argv[narg], "-autocorr", 5) == 0) {
00116                         narg++;
00117                         autocorr = TRUE;
00118                         continue;
00119                 }
00120 */
00121                 if (strncmp(argv[narg], "-retrend", 5) == 0) {
00122                         narg++;
00123                         retrend = TRUE;
00124                         continue;
00125                 }
00126 
00127                 if (strncmp(argv[narg], "-help", 5) == 0) {
00128                         help_message();
00129                 }
00130 
00131 
00132                 Error_Exit("Illegal or unrecoginized option");
00133 
00134         } /* end of while over arguements */
00135 
00136 #if 1
00137    if( low_fc > 0.0f && high_fc > 0.0f && low_fc <= high_fc )
00138      fprintf(stderr,"** WARNING: lowpass=%f is <= than highpass=%f\n"
00139                     "** -------: results from 3dFourier are suspect!\n"  ,
00140              low_fc,high_fc) ;
00141 #endif
00142 
00143         if( narg >= argc )
00144                 Error_Exit("No input datasets!?\n");
00145 
00146         if(prefix==NULL) {
00147                 prefix=(char *)My_Malloc(8*sizeof(char));
00148                 strcpy(prefix,"fourier");
00149         }
00150 
00151         /* try to open input dataset */
00152 
00153         input = THD_open_dataset( argv[narg] ) ; /* 16 Sep 1999 */
00154 
00155         if( input == NULL )
00156                 Error_Exit("Cannot open input dataset!") ;
00157 
00158         if (!DSET_GRAPHABLE(input))
00159                 Error_Exit("Input dataset is not a 3d+time dataset");
00160 
00161         if (DSET_BRICK_TYPE(input, 0)==MRI_complex)
00162                 Error_Exit("Sorry, I can't deal with complex 3d+time datasets");
00163 
00164         err=Fourier_Filter_Driver(input, low_fc, high_fc, ignore, autocorr, retrend, prefix);
00165         if (err!=NULL)
00166                 Error_Exit(err);
00167 
00168         return 0;  /* All went well exit code */
00169 }