MD++: My Documents/Codes/MD++.svn/trunk/src/xgeo.cpp Source File

00001 /*
00002   xgeo.cpp
00003   by Wei Cai  caiwei@mit.edu
00004   Last Modified : Thu Jan 11 11:58:45 2007
00005 
00006   FUNCTION  :
00007 */
00008 
00009 #include "xgeo.h"
00010 
00011 void XGeo::initparser()
00012 {
00013     bindvar("geo_box",geo_box,DOUBLE);
00014     bindvar("geo_line",geo_line,DOUBLE);
00015     bindvar("geo_plane",geo_plane,DOUBLE);
00016     bindvar("geo_xpoint",geo_xpoint,INT);
00017     bindvar("geo_xline",geo_xline,INT);
00018     bindvar("geo_vector",geo_vector,DOUBLE);
00019 
00020     bindvar("win_width",&win_width,INT);
00021     bindvar("win_height",&win_height,INT);
00022     bindvar("pointradius",&pointradius,DOUBLE);
00023     bindvar("geo_box_color",geo_box_color,STRING);
00024     bindvar("geo_line_color",geo_line_color,STRING);
00025     bindvar("geo_plane_color",geo_plane_color,STRING);
00026     bindvar("geo_xpoint_color",geo_xpoint_color,STRING);
00027     bindvar("geo_xline_color",geo_xline_color,STRING);
00028     bindvar("geo_vector_color",geo_vector_color,STRING);
00029     bindvar("backgroundcolor",backgroundcolor,STRING);
00030 }
00031 
00032 void XGeo::initvars()
00033 {
00034     initparser();
00035 }
00036 
00037 int XGeo::exec(char *name)
00038 {
00039     if(Organizer::exec(name)==0) return 0;
00040     bindcommand(name,"openwin",openwin());
00041     bindcommand(name,"plot",plot());
00042     bindcommand(name,"alloccolors",alloccolors());
00043     return -1;
00044 }
00045 
00046 #ifdef _USETCL
00047 /********************************************************************/
00048 /* Initialize Tcl parser */
00049 /********************************************************************/
00050 int XGeo::Tcl_AppInit(Tcl_Interp *interp)
00051 {
00052     /*
00053      * Tcl_Init reads init.tcl from the Tcl script library.
00054      */
00055     if (Tcl_Init(interp) == TCL_ERROR) {
00056         return TCL_ERROR;
00057     }
00058     /*
00059      * Register application-specific commands.
00060      */
00061     Tcl_CreateCommand(interp, "MD++", MD_Cmd,
00062                       (ClientData)NULL, (Tcl_CmdDeleteProc *)NULL);
00063 
00064     /*
00065       Random_Init(interp);
00066       Blob_Init(interp);
00067     */
00068     return TCL_OK;
00069 }
00070 
00071 int Tcl_Main_Wrapper(int argc, char *argv[])
00072 {
00073 #define MAXNARGS  20
00074 #define MAXARGLEN 500    
00075 //    char newargv[MAXNARGS][MAXARGLEN];
00076 
00077     if(argc>MAXNARGS)
00078         ERROR("maximum number of arguments exceeded!");
00079 
00080     Tcl_Main(argc, argv, Tcl_AppInit);
00081 //    Tcl_Main(argc, (char **)newargv, Tcl_AppInit);
00082 
00083     return TCL_OK;
00084 }
00085 #endif
00086 
00087 
00088 #include "namecolor.c"
00089 void XGeo::alloccolors()
00090 {
00091     int r,g,b;
00092     if ((win!=NULL)&&(win->alive))
00093     {
00094         Str2RGB(geo_box_color,&r,&g,&b);
00095         colors[0]=win->AllocShortRGBColor(r,g,b);
00096         Str2RGB(geo_line_color,&r,&g,&b);
00097         colors[1]=win->AllocShortRGBColor(r,g,b);
00098         Str2RGB(geo_plane_color,&r,&g,&b);
00099         colors[2]=win->AllocShortRGBColor(r,g,b);
00100         Str2RGB(geo_xpoint_color,&r,&g,&b);
00101         colors[3]=win->AllocShortRGBColor(r,g,b);
00102         Str2RGB(geo_xline_color,&r,&g,&b);
00103         colors[4]=win->AllocShortRGBColor(r,g,b);
00104         Str2RGB(geo_vector_color,&r,&g,&b);
00105         colors[5]=win->AllocShortRGBColor(r,g,b);
00106         Str2RGB(backgroundcolor,&r,&g,&b);
00107         win->bgcolor=win->AllocShortRGBColor(r,g,b);
00108     }
00109     else
00110     {
00111         WARNING("No window to allocate color for!");
00112     }
00113 }
00114 
00115 void XGeo::openwin()
00116 {
00117     openwindow(win_width,win_height,dirname);
00118     alloccolors();
00119 }
00120 
00121 int XGeo::openwindow(int w,int h,const char *n)
00122 {
00123     if(win!=NULL)
00124         if(win->alive) return -1;
00125     win=new YWindow(w,h,n);
00126     if(win==NULL) return -1;
00127     else
00128     {
00129         win->setinterval(100);
00130         win->Run();
00131         return 0;
00132     }
00133 }
00134 
00135 void XGeo::closewindow() {delete(win);}
00136 
00137 void XGeo::plot()
00138 {
00139     if((win!=NULL)&&(win->alive))
00140     {
00141         win->Lock();
00142         win->Clear();
00143         drawBox();
00144         drawLines();
00145         drawVectors();
00146         drawPlanes();
00147         drawXpoints();
00148         drawXlines();
00149     }
00150     win->Unlock();
00151     win->Refresh();
00152 }
00153 
00154 void XGeo::drawBox()
00155 {
00156     Vector3 a,b,c,e1,e2,e3,r1,r2;
00157     
00158     a.set(geo_box[0][0],geo_box[0][1],geo_box[0][2]);
00159     b.set(geo_box[1][0],geo_box[1][1],geo_box[1][2]);
00160     c.set(geo_box[2][0],geo_box[2][1],geo_box[2][2]);
00161     a*=geo_box[0][3];
00162     b*=geo_box[1][3];
00163     c*=geo_box[2][3];
00164     L=a.norm();
00165     if(L<b.norm()) L=b.norm();
00166     if(L<c.norm()) L=c.norm();
00167 
00168     L/=2;
00169 //    a/=L; b/=L; c/=L;
00170     _H.setcol(a,b,c);
00171 
00172     e1=a/a.norm();
00173     e2=b/b.norm(); e2=schmidt(e2,e1);
00174     e3=cross(e1,e2);
00175     _RV.setcol(e1,e2,e3); _RV=_RV.inv(); _RV/=L;
00176     
00177 #define _Geo_DoDraw(_1,_2,_3,_4,_5,_6) \
00178     r1=(a*(_1 1)+b*(_2 1)+c*(_3 1))*0.5; \
00179     r2=(a*(_4 1)+b*(_5 1)+c*(_6 1))*0.5; \
00180     r1=_RV*r1; r2=_RV*r2; \
00181     win->DrawLine(r1.x,r1.y,r1.z,r2.x,r2.y,r2.z,colors[0]);
00182 
00183     _Geo_DoDraw(+,-,+, +,+,+);
00184     _Geo_DoDraw(-,-,+, +,-,+);
00185     _Geo_DoDraw(-,+,+, -,-,+);
00186     _Geo_DoDraw(+,+,+, -,+,+);
00187     _Geo_DoDraw(+,-,-, +,+,-);
00188     _Geo_DoDraw(-,-,-, +,-,-);
00189     _Geo_DoDraw(-,+,-, -,-,-);
00190     _Geo_DoDraw(+,+,-, -,+,-);
00191     _Geo_DoDraw(+,+,+, +,+,-);
00192     _Geo_DoDraw(+,-,+, +,-,-);
00193     _Geo_DoDraw(-,-,+, -,-,-);
00194     _Geo_DoDraw(-,+,+, -,+,-);
00195 }
00196 
00197 void XGeo::drawLines()
00198 {
00199     int i, nlines;
00200     Vector3 dl, sr0, r1, r2;
00201     
00202     nlines=(int)geo_line[0];
00203     if(nlines<=0) return;
00204 
00205     for(i=0;i<nlines;i++)
00206     {
00207         dl.set(geo_line[i*7+2],geo_line[i*7+3],geo_line[i*7+4]);
00208         sr0.set(geo_line[i*7+5],geo_line[i*7+6],geo_line[i*7+7]);
00209         getlineendpoints(dl,sr0,r1,r2);
00210         r1=_RV*r1; r2=_RV*r2;
00211         win->DrawLine(r1.x,r1.y,r1.z,r2.x,r2.y,r2.z,colors[1],0.015);
00212     }
00213 }
00214 
00215 void XGeo::drawVectors()
00216 {
00217     int i, nvectors;
00218     Vector3 dl, sr0, r1, r2;
00219     double scale;
00220     
00221     nvectors=(int)geo_vector[0];
00222     if(nvectors<=0) return;
00223 
00224     for(i=0;i<nvectors;i++)
00225     {
00226         dl.set(geo_vector[i*8+2],geo_vector[i*8+3],geo_vector[i*8+4]);
00227         sr0.set(geo_vector[i*8+5],geo_vector[i*8+6],geo_vector[i*8+7]);
00228         scale=geo_vector[i*8+8];
00229         r1=_H*sr0;
00230         r2=r1+dl*scale/L/dl.norm();
00231         r1=_RV*r1; r2=_RV*r2;
00232         win->DrawLine(r1.x,r1.y,r1.z,r2.x,r2.y,r2.z,colors[5],0.02);
00233     }
00234 }
00235 
00236 bool XGeo::getlineendpoints(Vector3 &dl, Vector3 &sr0,
00237                             Vector3 &r1, Vector3 &r2)
00238 {
00239     Vector3 sl, sr[6], dr; Matrix33 hinv;
00240     double x, y, z, a;
00241     int i, n; bool succ;
00242 
00243 //    INFO("dl ="<<dl);
00244 //    INFO("sr0="<<sr0);
00245     hinv=_H.inv();
00246     sl=hinv*dl;
00247     x=y=z=0;
00248 
00249 #define testendpoint(_1,_2,_3)                             \
00250     if(sl._1!=0)                                           \
00251     {                                                      \
00252         a=(0.5-sr0._1)/sl._1;                              \
00253         _1=0.5;                                            \
00254         _2=a*sl._2+sr0._2;                                 \
00255         _3=a*sl._3+sr0._3;                                 \
00256         if((fabs(_2) <= 0.5001)&&(fabs(_3) <= 0.5001))     \
00257         {                                                  \
00258             sr[n].set(x,y,z); n++;                         \
00259         }                                                  \
00260         a=(-0.5-sr0._1)/sl._1;                             \
00261         _1=-0.5;                                           \
00262         _2=a*sl._2+sr0._2;                                 \
00263         _3=a*sl._3+sr0._3;                                 \
00264         if((fabs(_2) <= 0.5001)&&(fabs(_3) <= 0.5001))     \
00265         {                                                  \
00266             sr[n].set(x,y,z); n++;                         \
00267         }                                                  \
00268     }                                                      
00269 
00270     n=0;
00271     /* x=+-0.5 */
00272     testendpoint(x,y,z);
00273     
00274     /* y=+-0.5 */
00275     testendpoint(y,x,z);
00276 
00277     /* z=+-0.5 */
00278     testendpoint(z,x,y);
00279 
00280     if(n<2)
00281     {
00282         WARNING("getlineendpoints n="<<n<<" !=2");
00283         r1.clear(); r2.clear();
00284         return false;
00285     }
00286     
00287     succ=false;
00288     for(i=1;i<n;i++)
00289     {
00290         dr=sr[i]-sr[0];
00291         if(dr.norm2()>1e-6)
00292         {
00293             succ=true;
00294             break;
00295         }
00296     }
00297     sr[1]=sr[i];
00298     
00299     if(!succ)
00300     {
00301         WARNING("getlineendpoints failed");
00302         r1.clear(); r2.clear();
00303         return false;
00304     }
00305     else
00306     {
00307         r1=_H*sr[0];
00308         r2=_H*sr[1];
00309 //        INFO("sr[0]="<<sr[0]);
00310 //        INFO("sr[1]="<<sr[1]);
00311 //        INFO("r1="<<r1);
00312 //        INFO("r2="<<r2);
00313         return true;
00314     }
00315 }
00316 
00317 bool XGeo::getxline(Vector3 &dn1, Vector3 &sr1,
00318                     Vector3 &dn2, Vector3 &sr2,
00319                     Vector3 &dl,  Vector3 &sr)
00320 {
00321     Vector3 r1, r2, r;
00322     Vector3 e1, e2;
00323     Vector3 x, b;
00324     Matrix33 A, Ainv, hinv;
00325     
00326     hinv=_H.inv();
00327     r1=_H*sr1;
00328     r2=_H*sr2;            
00329     dl=cross(dn1, dn2);
00330 //    INFO("dl="<<dl);               
00331     if(dl.norm2()>1e-6)
00332     {                             
00333         e2=cross(dn2,dl); e1=cross(dn1,dl);               
00334         A.setcol(e1,e2,dl); Ainv=A.inv();                
00335         b=r2-r1;                                         
00336         x=Ainv*b;                                        
00337         r=r1+e1*x[0]; sr=hinv*r;
00338 //        INFO("getxline: sr="<<sr);
00339         return true;
00340     }
00341     else
00342     {
00343         sr.clear();
00344         return false;
00345     }
00346 }
00347 
00348 void XGeo::drawPlanes()
00349 {
00350     int i, nplanes;
00351     Vector3 dn, sr1, ra, rb;
00352     Vector3 dn2, sr2;
00353     Vector3 dl, srl;
00354     unsigned c; double alpha;
00355     
00356     nplanes=(int)geo_plane[0];
00357     if(nplanes<=0) return;
00358 
00359     for(i=0;i<nplanes;i++)
00360     {
00361         dn.set(geo_plane[i*7+2],geo_plane[i*7+3],geo_plane[i*7+4]);
00362         sr1.set(geo_plane[i*7+5],geo_plane[i*7+6],geo_plane[i*7+7]);
00363 
00364         alpha=1.0*i/nplanes/2;
00365         c=(unsigned)(colors[2]*(1-alpha)+colors[5]*alpha);
00366         
00367 #define testendline(_1,_2)                                \
00368         dn2.clear(); dn2._1=1;                            \
00369         sr2.clear(); sr2._1=_2 0.5;                       \
00370         if(getxline(dn,sr1,dn2,sr2,dl,srl)){              \
00371         if(getlineendpoints(dl,srl,ra,rb)){               \
00372         ra=_RV*ra; rb=_RV*rb;                             \
00373         win->DrawLine(ra.x,ra.y,ra.z,rb.x,rb.y,rb.z,c);}}
00374 
00375         /* x=+-0.5 */
00376         testendline(x,+);
00377         testendline(x,-);
00378         /* y=+-0.5 */
00379         testendline(y,+);
00380         testendline(y,-);
00381         /* z=+-0.5 */
00382         testendline(z,+);
00383         testendline(z,-);
00384     }
00385 }
00386 
00387 bool XGeo::getxpoint(Vector3 &dl1, Vector3 &sr1,
00388                      Vector3 &dn2, Vector3 &sr2, Vector3 &r)
00389 {
00390     Vector3 r1, r2, dr;
00391     double a, t;
00392     
00393     r1=_H*sr1;
00394     r2=_H*sr2;
00395 
00396 //    INFO("getxpoint: dl1="<<dl1<<"  r1="<<r1);
00397 //    INFO("getxpoint: dn2="<<dn2<<"  r2="<<r2);
00398 
00399     dr=r2-r1;
00400     t=dl1*dn2;
00401     if(fabs(t)<1e-6)
00402     {
00403         r.clear();
00404         return false;
00405     }
00406     else
00407     {
00408         a=(dr*dn2)/t;
00409         r=r1+dl1*a;
00410         return true;
00411     }
00412 }
00413 
00414 void XGeo::drawXpoints()
00415 {
00416     int i, nxpoints, l1, p2;
00417     Vector3 dl1, sr1;
00418     Vector3 dn2, sr2;
00419     Vector3 r;
00420     
00421     nxpoints=(int)geo_xpoint[0];
00422 //    INFO("nxpoints="<<nxpoints);
00423     if(nxpoints<=0) return;
00424 
00425     for(i=0;i<nxpoints;i++)
00426     {
00427         l1=geo_xpoint[i*3+2]; l1--;
00428         p2=geo_xpoint[i*3+3]; p2--;
00429 
00430 //        INFO("drawXlines: l1="<<l1<<"  p2="<<p2);
00431         
00432         dl1.set(geo_line[l1*7+2],geo_line[l1*7+3],geo_line[l1*7+4]);
00433         sr1.set(geo_line[l1*7+5],geo_line[l1*7+6],geo_line[l1*7+7]);
00434         dn2.set(geo_plane[p2*7+2],geo_plane[p2*7+3],geo_plane[p2*7+4]);
00435         sr2.set(geo_plane[p2*7+5],geo_plane[p2*7+6],geo_plane[p2*7+7]);
00436 
00437         if(getxpoint(dl1,sr1,dn2,sr2,r))
00438         {
00439 //            INFO("xpoint: r="<<r);
00440             r=_RV*r;
00441             win->DrawPoint(r.x,r.y,r.z,pointradius/L,colors[3]);
00442         }
00443     }
00444 }
00445 
00446 void XGeo::drawXlines()
00447 {
00448     int i, nxlines, p1, p2;
00449     Vector3 dn1, sr1, ra, rb;
00450     Vector3 dn2, sr2;
00451     Vector3 dl, sr;
00452     
00453     nxlines=(int)geo_xline[0];
00454 //    INFO("nxlines="<<nxlines);
00455     if(nxlines<=0) return;
00456 
00457     for(i=0;i<nxlines;i++)
00458     {
00459         p1=geo_xline[i*3+2]; p1--;
00460         p2=geo_xline[i*3+3]; p2--;
00461 
00462 //        INFO("drawXlines: p1="<<p1<<"  p2="<<p2);
00463         
00464         dn1.set(geo_plane[p1*7+2],geo_plane[p1*7+3],geo_plane[p1*7+4]);
00465         sr1.set(geo_plane[p1*7+5],geo_plane[p1*7+6],geo_plane[p1*7+7]);
00466         dn2.set(geo_plane[p2*7+2],geo_plane[p2*7+3],geo_plane[p2*7+4]);
00467         sr2.set(geo_plane[p2*7+5],geo_plane[p2*7+6],geo_plane[p2*7+7]);
00468 
00469         if(getxline(dn1,sr1,dn2,sr2,dl,sr))
00470         {
00471             if(getlineendpoints(dl,sr,ra,rb))
00472             {
00473                 ra=_RV*ra; rb=_RV*rb;                             
00474                 win->DrawLine(ra.x,ra.y,ra.z,rb.x,rb.y,rb.z,colors[4]);
00475             }
00476         }
00477     }
00478 }
00479 
00480 class XGeo sim;
00481 
00482 /* The main program is defined here */
00483 #include "main.cpp"
00484