00001 #include <odinpara/reco.h>
00002
00003 #include <odindata/image.h>
00004 #include <odindata/data.h>
00005 #include <odindata/gridding.h>
00006 #include <odindata/fitting.h>
00007
00014 void usage() {
00015 STD_cout << STD_endl;
00016 STD_cout << "gencoil: Generates a virtual coil suitable for sequence simulation in ODIN" << STD_endl;
00017 STD_cout << STD_endl;
00018 STD_cout << "Usage and options:" << STD_endl;
00019 STD_cout << " Generate a radially inhomogenous coil:" << STD_endl;
00020 STD_cout << " gencoil -rad -n <inplane-size> -fov <FOV> -R <radial-inhomogeneity[%]> -o <Coil-file>" << STD_endl;
00021 STD_cout << " Generate an array coil by rotating pre-existing image file:" << STD_endl;
00022 STD_cout << " gencoil -rot -n <inplane-size> -fov <FOV> -nc <numof-rotated-subcoils> -i <input-image> -sl <selected-slice> -o <Coil-file>" << STD_endl;
00023 STD_cout << " Generate an array coil consisting of simple loops:" << STD_endl;
00024 STD_cout << " gencoil -arr -n <inplane-size> -nc <numof-rotated-loops> -fov <FOV> -R <array-radius> -o <Coil-file>" << STD_endl;
00025 STD_cout << " Generate coil with a B1 gradient:" << STD_endl;
00026 STD_cout << " gencoil -grad -n <inplane-size> -fov <FOV> -g <B1-gradient-endpoint> -o <Coil-file>" << STD_endl;
00027 STD_cout << "Other options:" << STD_endl;
00028 STD_cout << "\t" << LogBase::get_usage() << STD_endl;
00029 STD_cout << "\t" << helpUsage() << STD_endl;
00030 }
00031
00032 int main(int argc, char* argv[]) {
00033 LogBase::set_log_levels(argc,argv);
00034 Log<Para> odinlog("gencoil","main");
00035
00036 if(hasHelpOption(argc,argv)) {usage(); return 0;}
00037
00038 Range all=Range::all();
00039
00040 char optval[ODIN_MAXCHAR];
00041
00042 STD_string coil_fname;
00043 if(getCommandlineOption(argc,argv,"-o",optval,ODIN_MAXCHAR)) coil_fname=optval; else {usage();exit(0);}
00044
00045 CoilSensitivity sens("Coil Sensitivity");
00046
00047 bool valid_mode=false;
00048
00049
00051
00052 if(isCommandlineOption(argc,argv,"-rad")) {
00053
00054 int n;
00055 float R,fov;
00056 if(getCommandlineOption(argc,argv,"-n",optval,ODIN_MAXCHAR)) n=atoi(optval); else {usage();exit(0);}
00057 if(getCommandlineOption(argc,argv,"-R",optval,ODIN_MAXCHAR)) R=0.01*atof(optval); else {usage();exit(0);}
00058 if(getCommandlineOption(argc,argv,"-fov",optval,ODIN_MAXCHAR)) fov=atof(optval); else {usage();exit(0);}
00059
00060 carray sens_map(1,1,n,n);
00061
00062 int nx_center=n/2;
00063 int ny_center=n/2;
00064
00065 double a=-4.0*R;
00066 double c=1.0+R;
00067
00068 STD_complex b1factor;
00069
00070 ndim nn=sens_map.get_extent();
00071 ndim indexvec;
00072
00073 for(int i=0; i<nn.total(); i++) {
00074 indexvec=nn.index2extent(i);
00075
00076 int ix=indexvec[3];
00077 int iy=indexvec[2];
00078
00079 double rx=secureDivision(fabs(double(ix-nx_center)),nx_center);
00080 double ry=secureDivision(fabs(double(iy-ny_center)),ny_center);
00081 double radius=sqrt(rx*rx+ry*ry);
00082
00083 double re=a*radius*radius+c;
00084 if(re<0.0) re=0.0;
00085
00086 b1factor=STD_complex(re,0.0);
00087 sens_map(indexvec)=b1factor;
00088 }
00089
00090 sens.set_sensitivity_map(sens_map,fov,fov,fov);
00091
00092 valid_mode=true;
00093 }
00094
00096
00097 if(isCommandlineOption(argc,argv,"-rot")) {
00098
00099 int n,nc,sl;
00100 float fov;
00101 STD_string input_fname;
00102 if(getCommandlineOption(argc,argv,"-n",optval,ODIN_MAXCHAR)) n=atoi(optval); else {usage();exit(0);}
00103 if(getCommandlineOption(argc,argv,"-nc",optval,ODIN_MAXCHAR)) nc=atoi(optval); else {usage();exit(0);}
00104 if(getCommandlineOption(argc,argv,"-fov",optval,ODIN_MAXCHAR)) fov=atof(optval); else {usage();exit(0);}
00105 if(getCommandlineOption(argc,argv,"-i",optval,ODIN_MAXCHAR)) input_fname=optval; else {usage();exit(0);}
00106 if(getCommandlineOption(argc,argv,"-sl",optval,ODIN_MAXCHAR)) sl=atoi(optval); else {usage();exit(0);}
00107
00108 Image img;
00109 ODINLOG(odinlog,infoLog) << "Loading image file ... " << STD_endl;
00110 if(img.load(input_fname)<0) {
00111 ODINLOG(odinlog,errorLog) << "unable to load ODIN image file " << input_fname << STD_endl;
00112 return -1;
00113 }
00114
00115 TinyMatrix<float,2,2> rotation;
00116 TinyVector<float,2> offset;
00117 offset=0.0;
00118
00119 carray sens_map(nc,1,n,n);
00120
00121 Data<float,4> imgdata(img.get_magnitude());
00122
00123 TinyVector<int,2> newshape(n,n);
00124
00125 ODINLOG(odinlog,infoLog) << "Congridding map ... " << STD_endl;
00126 Data<float,2> onecoil=imgdata(0,sl,all,all);
00127 onecoil.congrid(newshape);
00128 float maxmap=max(onecoil);
00129 float sum=0.0;
00130 float nvals=0.0;
00131 for(int i=0; i<onecoil.numElements(); i++) {
00132 float val=onecoil(onecoil.create_index(i));
00133 if(val>0.1*maxmap) {
00134 sum+=val;
00135 nvals+=1.0;
00136 }
00137 }
00138 float mean=sum/nvals;
00139 onecoil/=(2.0*mean);
00140
00141
00142 ODINLOG(odinlog,infoLog) << "Polynomial fit ... " << STD_endl;
00143 onecoil=polyniomial_fit(onecoil,onecoil,2,6.0);
00144
00145 Data<float,2> onecoil_transform(newshape);
00146
00147
00148 for(int ic=0; ic<nc; ic++) {
00149 ODINLOG(odinlog,infoLog) << "Rotating coil " << ic << " ... " << STD_endl;
00150 float phi=2.0*PII*float(ic)/float(nc);
00151 rotation(0,0)=cos(phi); rotation(0,1)=sin(phi);
00152 rotation(1,0)=-sin(phi); rotation(1,1)=cos(phi);
00153
00154 CoordTransformation<float,2> rotate(onecoil_transform.shape(), rotation, offset, 2.0);
00155 onecoil_transform=rotate(onecoil);
00156
00157 for(int iy=0; iy<n; iy++) {
00158 for(int ix=0; ix<n; ix++) {
00159 sens_map(ic,0,iy,ix)=STD_complex(onecoil_transform(iy,ix));
00160 }
00161 }
00162 }
00163
00164 sens.set_sensitivity_map(sens_map,fov,fov,fov);
00165
00166 valid_mode=true;
00167 }
00168
00169
00171
00172
00173 if(isCommandlineOption(argc,argv,"-arr")) {
00174
00175 int n,nc;
00176 float fov,R;
00177 if(getCommandlineOption(argc,argv,"-n",optval,ODIN_MAXCHAR)) n=atoi(optval); else {usage();exit(0);}
00178 if(getCommandlineOption(argc,argv,"-nc",optval,ODIN_MAXCHAR)) nc=atoi(optval); else {usage();exit(0);}
00179 if(getCommandlineOption(argc,argv,"-fov",optval,ODIN_MAXCHAR)) fov=atof(optval); else {usage();exit(0);}
00180 if(getCommandlineOption(argc,argv,"-R",optval,ODIN_MAXCHAR)) R=atof(optval); else {usage();exit(0);}
00181
00182
00183 carray sens_map(nc,1,n,n);
00184
00185 int nx_center=n/2;
00186 int ny_center=n/2;
00187
00188 ndim nn=sens_map.get_extent();
00189 ndim indexvec;
00190
00191 float overlap=2.0;
00192 float loopR=overlap*PII*R/float(nc);
00193
00194
00195 Array<float,2> rotmat[nc];
00196 Array<float,1> offset[nc];
00197
00198 ODINLOG(odinlog,infoLog) << "Preparing rotmat/offset ... " << STD_endl;
00199 for(int icoil=0; icoil<nc; icoil++) {
00200 float ang=2.0*PII*float(icoil)/float(nc);
00201
00202 rotmat[icoil].resize(3,3);
00203 rotmat[icoil]=0.0;
00204 rotmat[icoil](0,0)=cos(ang); rotmat[icoil](0,1)=-sin(ang);
00205 rotmat[icoil](1,0)=sin(ang); rotmat[icoil](1,1)=cos(ang);
00206 rotmat[icoil](2,2)=1.0;
00207
00208 offset[icoil].resize(3);
00209 offset[icoil](0)=R*cos(ang);
00210 offset[icoil](1)=R*sin(ang);
00211 offset[icoil](2)=0.0;
00212 }
00213
00214
00215 ODINLOG(odinlog,infoLog) << "Calculating sens_map ..." << STD_endl;
00216 for(int i=0; i<nn.total(); i++) {
00217 indexvec=nn.index2extent(i);
00218
00219 int icoil=indexvec[0];
00220 int iy=indexvec[2];
00221 int ix=indexvec[3];
00222
00223 Array<float,1> x(3);
00224 x(0)=fov*secureDivision(double(ix-nx_center),nx_center);
00225 x(1)=fov*secureDivision(double(iy-ny_center),ny_center);
00226 x(2)=0.0;
00227
00228
00229 int nloopsegments=100;
00230 float dphi=2.0*PII/nloopsegments;
00231 Array<float,1> dl(3);
00232 Array<float,1> r(3);
00233 Array<float,1> r0(3);
00234 Array<float,1> dB(3);
00235 STD_complex B1(0.0);
00236 for(int iseg=0; iseg<nloopsegments; iseg++) {
00237 float phi=2.0*PII*float(iseg)/float(nloopsegments);
00238 dl(0)=0.0;
00239 dl(1)= loopR*cos(phi)*dphi;
00240 dl(2)=-loopR*sin(phi)*dphi;
00241
00242 dl=matrix_product(rotmat[icoil], dl);
00243
00244 r0(0)=0.0;
00245 r0(1)=loopR*sin(phi);
00246 r0(2)=loopR*cos(phi);
00247
00248 r0=matrix_product(rotmat[icoil], r0) + offset[icoil];
00249
00250 r=x-r0;
00251
00252 dB=vector_product(dl,r);
00253
00254 float r2=sum(r*r);
00255 if(r2>0.0) {
00256 B1+=STD_complex(dB(0),dB(1))/r2;
00257 }
00258
00259 }
00260 sens_map(indexvec)=B1/float(2.0*PII);
00261 }
00262
00263 sens.set_sensitivity_map(sens_map,fov,fov,fov);
00264
00265
00266
00267 valid_mode=true;
00268 }
00269
00270
00272
00273 if(isCommandlineOption(argc,argv,"-grad")) {
00274
00275 int n;
00276 float fov,grad;
00277 if(getCommandlineOption(argc,argv,"-n",optval,ODIN_MAXCHAR)) n=atoi(optval); else {usage();exit(0);}
00278 if(getCommandlineOption(argc,argv,"-g",optval,ODIN_MAXCHAR)) grad=atof(optval); else {usage();exit(0);}
00279 if(getCommandlineOption(argc,argv,"-fov",optval,ODIN_MAXCHAR)) fov=atof(optval); else {usage();exit(0);}
00280
00281 carray sens_map(1,1,n,n);
00282
00283 int nx_center=n/2;
00284
00285 ndim nn=sens_map.get_extent();
00286 ndim indexvec;
00287
00288 for(int i=0; i<nn.total(); i++) {
00289 indexvec=nn.index2extent(i);
00290 int ix=indexvec[3];
00291 sens_map(indexvec)=STD_complex(1.0+grad*float(ix-nx_center)/float(nx_center),0.0);
00292 }
00293
00294 sens.set_sensitivity_map(sens_map,fov,fov,fov);
00295
00296 valid_mode=true;
00297 }
00298
00299
00300
00301 if(!valid_mode) {usage();exit(0);}
00302
00303 ODINLOG(odinlog,infoLog) << "Writing coil file ... " << STD_endl;
00304 sens.write(coil_fname);
00305
00306 return 0;
00307 }
