00001 #include <odindata/data.h>
00002 #include <odindata/fileio.h>
00003 #include <odindata/utils.h>
00004 #include <odindata/statistics.h>
00005 #include <odindata/correlation.h>
00006 #include <odindata/filter.h>
00007
00062 void usage(const Protocol& prot) {
00063 FileReadOpts ropts;
00064 FileWriteOpts wopts;
00065 FilterFactory filter;
00066 STD_cout << "micalc: Performs basic mathematics with data sets" << STD_endl;
00067 STD_cout << " File formats are automatically identified by their file extension." << STD_endl;
00068 STD_cout << "Usage:" << STD_endl;
00069 STD_cout << " Binary operation with data sets and/or scalar numbers:" << STD_endl;
00070 STD_cout << " micalc [-if1 <input-file1> | -in1 <input-number1>] -op <operation(+,-,*,/,lcorr,kcorr)> [-if2 <input-file2> | -in2 <input-number1>] -of <output-file>" << STD_endl;
00071 STD_cout << " Accumulation of one data set:" << STD_endl;
00072 STD_cout << " micalc -if <input-file> -op <operation(+,*)>" << STD_endl;
00073 STD_cout << " Transformation (magnitude, logarithm, negation, inversion, or acos) of one data set:" << STD_endl;
00074 STD_cout << " micalc -if <input-file> -op <operation(abs,log,-,/,acos)> -of <output-file>" << STD_endl;
00075 STD_cout << " Statistics of one data set:" << STD_endl;
00076 STD_cout << " micalc -if <input-file> [-mean <time-mean-file>] [-stdev <time-stdev-file>] [-fluct <relative-time-stdev-file>] [-tcourse <time-course of all voxels>] [-hist <histogram> -histmin <minval> -histmax <maxval> -histslots <numofslots> -rightstairs -histfract]" << STD_endl;
00077 STD_cout << "Extra options:" << STD_endl;
00078 STD_cout << "\t-mask <Mask for data>" << STD_endl;
00079 STD_cout << "File read options:" << STD_endl;
00080 STD_cout << prot.get_cmdline_usage("\t");
00081 STD_cout << ropts.get_cmdline_usage("\t");
00082 STD_cout << "File write options:" << STD_endl;
00083 STD_cout << wopts.get_cmdline_usage("\t");
00084 STD_cout << "Filters applied to input files:" << STD_endl;
00085 STD_cout << filter.get_cmdline_usage("\t");
00086 STD_cout << "Other options:" << STD_endl;
00087 STD_cout << "\t" << LogBase::get_usage() << STD_endl;
00088 STD_cout << "\t" << helpUsage() << STD_endl;
00089 STD_cout << "Supported file extensions(formats):" << STD_endl;
00090 STD_cout << FileIO::autoformats_str("\t") << STD_endl;
00091 }
00092
00093
00094 int main(int argc, char* argv[]) {
00095 LogBase::set_log_levels(argc,argv);
00096
00097 Log<OdinData> odinlog("micalc","main");
00098
00099 Range all=Range::all();
00100
00101 FileReadOpts ropts;
00102 FileWriteOpts wopts;
00103 Protocol prot;
00104
00105 if(hasHelpOption(argc,argv)) {usage(prot); return 0;}
00106
00107
00108 prot.seqpars.set_MatrixSize(readDirection,1);
00109 prot.seqpars.set_MatrixSize(phaseDirection,1);
00110 prot.seqpars.set_MatrixSize(sliceDirection,1);
00111
00112 char optval[ODIN_MAXCHAR];
00113 STD_string operation;
00114 STD_string infile1;
00115 STD_string infile2;
00116 STD_string outfile;
00117 STD_string maskfile;
00118
00119 float innumber1=0.0;
00120 float innumber2=0.0;
00121
00122
00123 bool i1flag=false;
00124 bool i2flag=false;
00125 bool offlag=false;
00126
00127 bool ifflag=false;
00128 bool opflag=false;
00129
00130 bool binaryop=false;
00131 bool accumulate=false;
00132 bool transform=false;
00133 bool stateval=false;
00134
00135
00136 if(getCommandlineOption(argc,argv,"-op",optval,ODIN_MAXCHAR)) {operation=optval; opflag=true;}
00137
00138 if(getCommandlineOption(argc,argv,"-if1",optval,ODIN_MAXCHAR)) {infile1=optval; i1flag=true;}
00139 if(getCommandlineOption(argc,argv,"-if2",optval,ODIN_MAXCHAR)) {infile2=optval; i2flag=true;}
00140 if(getCommandlineOption(argc,argv,"-if",optval,ODIN_MAXCHAR)) {infile1=optval; ifflag=true;}
00141 if(getCommandlineOption(argc,argv,"-in1",optval,ODIN_MAXCHAR)) {innumber1=atof(optval); i1flag=true;}
00142 if(getCommandlineOption(argc,argv,"-in2",optval,ODIN_MAXCHAR)) {innumber2=atof(optval); i2flag=true;}
00143
00144 if(getCommandlineOption(argc,argv,"-of",optval,ODIN_MAXCHAR)) {outfile=optval; offlag=true;}
00145
00146 if(getCommandlineOption(argc,argv,"-mask",optval,ODIN_MAXCHAR)) {maskfile=optval;}
00147
00148 STD_string histfile;
00149 if(getCommandlineOption(argc,argv,"-hist",optval,ODIN_MAXCHAR)) histfile=optval;
00150
00151 STD_string meanfile;
00152 STD_string fluctfile;
00153 STD_string stdevfile;
00154 STD_string tcoursefile;
00155 STD_string histslots;
00156 STD_string histmin;
00157 STD_string histmax;
00158 if(getCommandlineOption(argc,argv,"-mean",optval,ODIN_MAXCHAR)) meanfile=optval;
00159 if(getCommandlineOption(argc,argv,"-fluct",optval,ODIN_MAXCHAR)) fluctfile=optval;
00160 if(getCommandlineOption(argc,argv,"-stdev",optval,ODIN_MAXCHAR)) stdevfile=optval;
00161 if(getCommandlineOption(argc,argv,"-tcourse",optval,ODIN_MAXCHAR)) tcoursefile=optval;
00162 if(getCommandlineOption(argc,argv,"-histslots",optval,ODIN_MAXCHAR)) histslots=optval;
00163 if(getCommandlineOption(argc,argv,"-histmin",optval,ODIN_MAXCHAR)) histmin=optval;
00164 if(getCommandlineOption(argc,argv,"-histmax",optval,ODIN_MAXCHAR)) histmax=optval;
00165 bool rightstairs=isCommandlineOption(argc,argv,"-rightstairs");
00166 bool histfract=isCommandlineOption(argc,argv,"-histfract");
00167
00168
00169
00170
00171 ropts.parse_cmdline_options(argc,argv);
00172 wopts.parse_cmdline_options(argc,argv);
00173 prot.parse_cmdline_options(argc,argv);
00174
00175 FilterChain filterchain(argc,argv);
00176
00177
00178 if(i1flag && i2flag && opflag) binaryop=true;
00179 if(ifflag && opflag) {
00180 if(offlag) transform=true;
00181 else accumulate=true;
00182 }
00183 if(ifflag && (!opflag)) stateval=true;
00184
00185 if( !( binaryop || accumulate || transform || stateval ) ) {usage(prot); return 0;}
00186
00187 if(accumulate || stateval) FileIO::set_trace_status(false);
00188
00189 if(binaryop && infile1=="" && infile2=="") {
00190 ODINLOG(odinlog,errorLog) << "Specifiy at least one file" << STD_endl;
00191 return -1;
00192 }
00193
00194 Data<float,4> data1;
00195 Data<float,4> data2;
00196 Data<float,4> datain;
00197
00198 if(infile1!="") {
00199 if(datain.autoread(infile1,ropts,&prot)<0) return -1;
00200 if(!filterchain.apply(prot,datain)) return -1;
00201 data1.resize(datain.shape());
00202 data1=datain;
00203 }
00204
00205 if(binaryop && infile2!="") {
00206 if(datain.autoread(infile2,ropts,&prot)<0) return -1;
00207 if(!filterchain.apply(prot,datain)) return -1;
00208 data2.resize(datain.shape());
00209 data2=datain;
00210 }
00211
00212 TinyVector<int,4> outshape;
00213
00214
00215 if(binaryop || transform) {
00216
00217 bool has_shape1=sum(data1.shape());
00218 bool has_shape2=sum(data2.shape());
00219
00220
00221 if( has_shape1 && has_shape2 && !same_shape(data1,data2) ) {
00222 ODINLOG(odinlog,errorLog) << "shape mismatch:" << STD_endl;
00223 ODINLOG(odinlog,errorLog) << "data1.shape()=" << data1.shape() << STD_endl;
00224 ODINLOG(odinlog,errorLog) << "data2.shape()=" << data2.shape() << STD_endl;
00225 return -1;
00226 } else outshape=data1.shape();
00227
00228 if(has_shape1 && !has_shape2) {
00229 outshape=data1.shape();
00230 data2.resize(outshape);
00231 data2=innumber2;
00232 }
00233 if(!has_shape1 && has_shape2) {
00234 outshape=data2.shape();
00235 data1.resize(outshape);
00236 data1=innumber1;
00237 }
00238 if(!has_shape1 && !has_shape2) {
00239 ODINLOG(odinlog,errorLog) << "Specifiy at least one valid file" << STD_endl;
00240 return -1;
00241 }
00242
00243 }
00244
00245
00246 Data<float,4>* maskptr=0;
00247 Data<float,4> maskdata;
00248 Data<float,4> inmaskdata;
00249 if(maskfile!="") {
00250 if(inmaskdata.autoread(maskfile,ropts,&prot)<0) return -1;
00251 TinyVector<int,4> datashape=data1.shape(); datashape(0)=1;
00252 TinyVector<int,4> maskshape=inmaskdata.shape(); maskshape(0)=1;
00253 if(datashape!=maskshape) {
00254 ODINLOG(odinlog,errorLog) << "Shape mismatch: datashape/maskshape" << datashape << "/" << maskshape << STD_endl;
00255 return -1;
00256 } else {
00257 maskdata.resize(data1.shape());
00258 maskdata=0.0;
00259 for(int i=0; i<data1.numElements(); i++) {
00260 TinyVector<int,4> index=data1.create_index(i);
00261 TinyVector<int,4> maskindex=index; maskindex(0)=0;
00262 if(inmaskdata(maskindex)>0.0) maskdata(index)=1.0;
00263 }
00264 maskptr=&maskdata;
00265 }
00266 }
00267
00268
00269 prot.system.set_data_type(TypeTraits::type2label(float(0)));
00270
00271 if(binaryop) {
00272 if(maskptr) {
00273 ODINLOG(odinlog,normalDebug) << "data1/data2/maskdata" << data1.shape() << "/" << data2.shape() << "/" << maskdata.shape() << STD_endl;
00274 data1=data1*maskdata;
00275 data2=data2*maskdata;
00276 }
00277
00278 if(operation=="lcorr" || operation=="kcorr") {
00279 int npts=data1.size();
00280 if(maskptr) npts=int(sum(*maskptr)+0.5);
00281 Array<float,1> vec1(npts);
00282 Array<float,1> vec2(npts);
00283 int vecindex=0;
00284 for(int i=0; i<data1.numElements(); i++) {
00285 TinyVector<int,4> index=data1.create_index(i);
00286 bool include=true;
00287 if( maskptr && (*maskptr)(index)<=0.0 ) include=false;
00288 if(include) {
00289 vec1(vecindex)=data1(index);
00290 vec2(vecindex)=data2(index);
00291 vecindex++;
00292 }
00293 }
00294 correlationResult coresult;
00295 if(operation=="lcorr") coresult=correlation(vec1, vec2);
00296 if(operation=="kcorr") coresult=kendall(vec1, vec2);
00297 STD_cout << "r/p/z(" << npts << ")=" << coresult.r << "/" << coresult.p << "/" << coresult.z << STD_endl;
00298 } else {
00299 Data<float,4> result(outshape);
00300 ODINLOG(odinlog,infoLog) << "calculating " << data1.shape() << " " << operation << " " << data2.shape() << " ..." << STD_endl;
00301 bool validop=false;
00302 if(operation=="+") {validop=true; result=data1+data2;}
00303 if(operation=="-") {validop=true; result=data1-data2;}
00304 if(operation=="*") {validop=true; result=data1*data2;}
00305 if(operation=="/") {validop=true; result=secureDivision(data1,data2);}
00306 if(!validop) {
00307 ODINLOG(odinlog,errorLog) << "invalid binary operation: " << operation << STD_endl;
00308 return -1;
00309 }
00310 if(result.autowrite(outfile,wopts,&prot)<0) return -1;
00311 }
00312 return 0;
00313 }
00314
00315 if(accumulate) {
00316 if(maskptr) {
00317 data1=data1*maskdata;
00318 }
00319 double accuresult=0.0;
00320 bool validop=false;
00321 if(operation=="+") {validop=true; accuresult=sum(data1);}
00322 if(operation=="*") {validop=true; accuresult=product(data1);}
00323 if(!validop) {
00324 ODINLOG(odinlog,errorLog) << "invalid accumulation operator: " << operation << STD_endl;
00325 return -1;
00326 }
00327 STD_cout << accuresult << STD_endl;
00328 }
00329
00330 if(transform) {
00331 if(maskptr) {
00332 data1=data1*maskdata;
00333 }
00334 Data<float,4> result(outshape);
00335 ODINLOG(odinlog,infoLog) << "Calculating " << operation << " " << data1.shape() << " ..." << STD_endl;
00336 bool validop=false;
00337 if(operation=="abs") {validop=true; result=fabs(data1);}
00338 if(operation=="log") {validop=true; result=where(Array<float,4>(data1)>0.0, Array<float,4>(log(data1)), float(0.0));}
00339 if(operation=="-") {validop=true; result=-data1;}
00340 if(operation=="/") {validop=true; result=where(Array<float,4>(data1)!=0.0, Array<float,4>(1.0/data1), float(0.0));}
00341 if(operation=="acos") {validop=true; result=where(Array<float,4>(data1)>=-1.0 && Array<float,4>(data1)<=1.0, Array<float,4>(acos(data1)), float(0.0));}
00342 if(!validop) {
00343 ODINLOG(odinlog,errorLog) << "invalid unary operation: " << operation << STD_endl;
00344 return -1;
00345 }
00346 if(result.autowrite(outfile,wopts,&prot)<0) return -1;
00347 return 0;
00348 }
00349
00350 if(stateval) {
00351
00352 if(histfile!="") {
00353
00354 int nvals=data1.numElements();
00355 if(maskptr) {
00356 nvals=0;
00357 for(int i=0; i<maskptr->numElements(); i++) {
00358 TinyVector<int,4> index=maskptr->create_index(i);
00359 if( (*maskptr)(index)>0.0 ) {
00360 nvals++;
00361 }
00362 }
00363 }
00364
00365 if(histslots=="") {usage(prot); return -1;}
00366
00367 int nslots=atoi(histslots.c_str());
00368
00369 if(histmin=="") {usage(prot); return -1;}
00370 float minval=atof(histmin.c_str());
00371 if(histmax=="") {usage(prot); return -1;}
00372 float maxval=atof(histmax.c_str());
00373 ODINLOG(odinlog,infoLog) << "calculating histogram with nslots/minval/maxval=" << nslots << "/" << minval << "/" << maxval << STD_endl;
00374
00375 if(minval>=maxval) {
00376 ODINLOG(odinlog,errorLog) << "histmax must be larger than histmin" << STD_endl;
00377 return -1;
00378 }
00379
00380 float stairoffset=0.5;
00381 if(rightstairs) stairoffset=0.0;
00382
00383 Data<float,1> xvals(nslots);
00384 Data<float,1> yvals(nslots); yvals=0.0;
00385 if(minval>0.0 && maxval>0.0) {
00386
00387 float minpow=log10(minval);
00388 float maxpow=log10(maxval);
00389 float powrange=maxpow-minpow;
00390 float step=secureDivision(powrange,nslots);
00391 ODINLOG(odinlog,normalDebug) << "minpow/maxpow/powrange/step=" << minpow << "/" << maxpow << "/" << powrange << "/" << step << STD_endl;
00392
00393 for(int i=0; i<nslots; i++) xvals(i)=pow(float(10.0),float(minpow+(float(i)+stairoffset)*step));
00394
00395 for(int i=0; i<data1.numElements(); i++) {
00396 TinyVector<int,4> index=data1.create_index(i);
00397 if( !(maskptr && (*maskptr)(index)==0.0) ) {
00398 float pow=log10(data1(index));
00399 int slot=int(secureDivision(pow-minpow,powrange)*nslots);
00400 ODINLOG(odinlog,normalDebug) << "pow/slot=" << pow << "/" << slot << STD_endl;
00401 if(slot>=0 && slot<nslots) yvals(slot)++;
00402 }
00403 }
00404
00405
00406 } else {
00407
00408 float range=maxval-minval;
00409 float step=secureDivision(range,nslots);
00410 ODINLOG(odinlog,normalDebug) << "range/step=" << range << "/" << step << STD_endl;
00411
00412 for(int i=0; i<nslots; i++) xvals(i)=minval+(float(i)+stairoffset)*step;
00413
00414 for(int i=0; i<data1.numElements(); i++) {
00415 TinyVector<int,4> index=data1.create_index(i);
00416 if( !(maskptr && (*maskptr)(index)==0.0) ) {
00417 float val=data1(index);
00418 int slot=int(secureDivision(val-minval,range)*nslots);
00419 ODINLOG(odinlog,normalDebug) << "val/slot=" << val << "/" << slot << STD_endl;
00420 if(slot>=0 && slot<nslots) yvals(slot)++;
00421 }
00422 }
00423 }
00424
00425 if(histfract) yvals/=sum(yvals);
00426 yvals.write_asc_file(histfile, xvals);
00427
00428 }
00429
00430
00431 if(meanfile!="" || fluctfile!="" || stdevfile!="" || tcoursefile!="") {
00432 outshape=data1.shape();
00433 ODINLOG(odinlog,infoLog) << "calculating mean, standard deviation and/or timecourse for shape " << outshape << " ..." << STD_endl;
00434 outshape(0)=1;
00435 Data<float,4> meanresult(outshape); meanresult=0.0;
00436 Data<float,4> fluctresult(outshape); fluctresult=0.0;
00437 Data<float,4> stdevresult(outshape); stdevresult=0.0;
00438 Data<float,1> tcourseresult(data1.extent(0)); tcourseresult=0.0;
00439 int tcoursepts=0;
00440 for(int i=0; i<meanresult.numElements(); i++) {
00441 TinyVector<int,4> index=meanresult.create_index(i);
00442 if(!index(0)) {
00443 if( !(maskptr && (*maskptr)(index)==0.0) ) {
00444 statisticResult stats=statistics(data1(all,index(1),index(2),index(3)));
00445 meanresult(index)= stats.mean;
00446 fluctresult(index)=secureDivision(stats.stdev,stats.mean);
00447 stdevresult(index)=stats.stdev;
00448 tcourseresult(all)+=data1(all,index(1),index(2),index(3));
00449 tcoursepts++;
00450 ODINLOG(odinlog,normalDebug) << "mean/stdev(" << index << ")=" << stats.mean << "/" << stats.stdev << STD_endl;
00451 }
00452 }
00453 }
00454 tcourseresult/=float(tcoursepts);
00455 if(meanfile!="") if(meanresult.autowrite(meanfile,wopts,&prot)<0) return -1;
00456 if(fluctfile!="") if(fluctresult.autowrite(fluctfile,wopts,&prot)<0) return -1;
00457 if(stdevfile!="") if(stdevresult.autowrite(stdevfile,wopts,&prot)<0) return -1;
00458 if(tcoursefile!="") if(tcourseresult.autowrite(tcoursefile,wopts,&prot)<0) return -1;
00459 } else {
00460 STD_cout << statistics(data1,maskptr) << STD_endl;
00461 }
00462
00463 }
00464
00465 return 0;
00466 }
