00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "symmetry.h"
00033 #include "transform.h"
00034 #include "vec3.h"
00035 #include "exception.h"
00036 #include "util.h"
00037
00038 using namespace EMAN;
00039
00040 const string CSym::NAME = "c";
00041 const string DSym::NAME = "d";
00042 const string HSym::NAME = "h";
00043 const string TetrahedralSym::NAME = "tet";
00044 const string OctahedralSym::NAME = "oct";
00045 const string IcosahedralSym::NAME = "icos";
00046 const string Icosahedral2Sym::NAME = "icos2";
00047 const string EmanOrientationGenerator::NAME = "eman";
00048 const string SaffOrientationGenerator::NAME = "saff";
00049 const string EvenOrientationGenerator::NAME = "even";
00050 const string RandomOrientationGenerator::NAME = "rand";
00051 const string OptimumOrientationGenerator::NAME = "opt";
00052
00053 template <> Factory < Symmetry3D >::Factory()
00054 {
00055 force_add<CSym>();
00056 force_add<DSym>();
00057 force_add<HSym>();
00058 force_add<TetrahedralSym>();
00059 force_add<OctahedralSym>();
00060 force_add<IcosahedralSym>();
00061 force_add<Icosahedral2Sym>();
00062 }
00063
00064 void EMAN::dump_symmetries()
00065 {
00066 dump_factory < Symmetry3D > ();
00067 }
00068
00069 map<string, vector<string> > EMAN::dump_symmetries_list()
00070 {
00071 return dump_factory_list < Symmetry3D > ();
00072 }
00073
00074 template <>
00075 Symmetry3D* Factory < Symmetry3D >::get(const string & instancename_)
00076 {
00077 init();
00078
00079 string instancename = Util::str_to_lower(instancename_);
00080
00081 unsigned int n = instancename.size();
00082 if ( n == 0 ) throw NotExistingObjectException(instancename, "Empty instance name!");
00083
00084 char leadingchar = instancename[0];
00085 if (leadingchar == 'c' || leadingchar == 'd' || leadingchar == 'h' ) {
00086 Dict parms;
00087 if (n > 1) {
00088 int nsym = atoi(instancename.c_str() + 1);
00089 parms["nsym"] = nsym;
00090 }
00091
00092 if (leadingchar == 'c') {
00093 return get("c",parms);
00094 }
00095 if (leadingchar == 'd') {
00096 return get("d",parms);
00097 }
00098 if (leadingchar == 'h') {
00099 int nstart=1,nsym=1;
00100 float daz=5.,tz=5.,maxtilt=1.0;
00101 string temp;
00102 temp=instancename;
00103 temp.erase(0,1);
00104 sscanf(temp.c_str(),"%d,%d,%f,%f,%f",&nsym,&nstart,&daz,&tz,&maxtilt);
00105 parms["nstart"]=nstart;
00106 parms["nsym"]=nsym;
00107 parms["daz"]=daz;
00108 parms["tz"]=tz;
00109 parms["maxtilt"]=maxtilt;
00110 return get("h",parms);
00111 }
00112
00113
00114 }
00115 else if ( instancename == "icos" || instancename == "oct" || instancename == "tet" || instancename == "icos2" || instancename == "icos5" )
00116 {
00117 if (instancename == "icos5") {
00118 instancename = "icos";
00119 }
00120
00121 map < string, InstanceType >::iterator fi =
00122 my_instance->my_dict.find(instancename);
00123 if (fi != my_instance->my_dict.end()) {
00124 return my_instance->my_dict[instancename] ();
00125 }
00126
00127 string lower = instancename;
00128 for (unsigned int i=0; i<lower.length(); i++) lower[i]=tolower(lower[i]);
00129
00130 fi = my_instance->my_dict.find(lower);
00131 if (fi != my_instance->my_dict.end()) {
00132 return my_instance->my_dict[lower] ();
00133 }
00134
00135 throw NotExistingObjectException(instancename, "No such an instance existing");
00136 }
00137 else throw NotExistingObjectException(instancename, "No such an instance existing");
00138
00139 throw NotExistingObjectException(instancename, "No such an instance existing");
00140 }
00141
00142 template <> Factory < OrientationGenerator >::Factory()
00143 {
00144 force_add<EmanOrientationGenerator>();
00145 force_add<RandomOrientationGenerator>();
00146 force_add<EvenOrientationGenerator>();
00147 force_add<SaffOrientationGenerator>();
00148 force_add<OptimumOrientationGenerator>();
00149 }
00150
00151
00152
00153 void EMAN::dump_orientgens()
00154 {
00155 dump_factory < OrientationGenerator > ();
00156 }
00157
00158 map<string, vector<string> > EMAN::dump_orientgens_list()
00159 {
00160 return dump_factory_list < OrientationGenerator > ();
00161 }
00162
00163 vector<Transform> Symmetry3D::gen_orientations(const string& generatorname, const Dict& parms)
00164 {
00165 ENTERFUNC;
00166 vector<Transform> ret;
00167 OrientationGenerator *g = Factory < OrientationGenerator >::get(Util::str_to_lower(generatorname), parms);
00168 if (g) {
00169 ret = g->gen_orientations(this);
00170 if( g )
00171 {
00172 delete g;
00173 g = 0;
00174 }
00175 }
00176 else throw;
00177
00178 EXITFUNC;
00179
00180 return ret;
00181 }
00182
00183 void OrientationGenerator::get_az_max(const Symmetry3D* const sym, const float& altmax, const bool inc_mirror, const float& alt_iterator,const float& h,bool& d_odd_mirror_flag, float& azmax_adjusted) const
00184 {
00185
00186 if ( sym->is_d_sym() && alt_iterator == altmax && ( (sym->get_nsym())/2 % 2 == 1 )) {
00187 if (inc_mirror) {
00188 azmax_adjusted /= 4.0f;
00189 d_odd_mirror_flag = true;
00190 }
00191 else azmax_adjusted /= 2.0f;
00192 }
00193 else if (sym->is_d_sym() && alt_iterator == altmax && ( (sym->get_nsym())/2 % 2 == 0 ) && inc_mirror) {
00194 azmax_adjusted /= 2.0f;
00195 }
00196
00197
00198 else if (sym->is_c_sym() && !inc_mirror && alt_iterator == altmax && (sym->get_nsym() % 2 == 1 ) ){
00199 azmax_adjusted /= 2.0f;
00200 }
00201
00202
00203 else if (sym->is_c_sym() || sym->is_tet_sym() ) {
00204 azmax_adjusted -= h/4.0f;
00205 }
00206
00207
00208 else if (inc_mirror && ( sym->is_d_sym() || sym->is_platonic_sym() ) ) {
00209 azmax_adjusted -= h/4.0f;
00210 }
00211
00212
00213
00214
00215
00216 }
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232 float OrientationGenerator::get_optimal_delta(const Symmetry3D* const sym, const int& n) const
00233 {
00234
00235
00236 float delta_soln = 180.0f;
00237 float delta_upper_bound = delta_soln;
00238 float delta_lower_bound = 0.0f;
00239
00240 int prev_tally = -1;
00241
00242
00243
00244 bool soln_found = false;
00245
00246 while ( soln_found == false ) {
00247 int tally = get_orientations_tally(sym,delta_soln);
00248 if ( tally == n ) soln_found = true;
00249 else if ( (delta_upper_bound - delta_lower_bound) < 0.0001 ) {
00250
00251
00252 soln_found = true;
00253 delta_soln = (delta_upper_bound+delta_lower_bound)/2.0f;
00254 }
00255 else if (tally < n) {
00256 delta_upper_bound = delta_soln;
00257 delta_soln = delta_soln - (delta_soln-delta_lower_bound)/2.0f;
00258 }
00259 else {
00260 delta_lower_bound = delta_soln;
00261 delta_soln = delta_soln + (delta_upper_bound-delta_soln)/2.0f;
00262 }
00263 prev_tally = tally;
00264 }
00265
00266 return delta_soln;
00267 }
00268
00269 bool OrientationGenerator::add_orientation(vector<Transform>& v, const float& az, const float& alt) const
00270 {
00271 bool randphi = params.set_default("random_phi",false);
00272 float phi = 0.0f;
00273 if (randphi) phi = Util::get_frand(0.0f,359.99999f);
00274 float phitoo = params.set_default("phitoo",0.0f);
00275 if ( phitoo < 0 ) throw InvalidValueException(phitoo, "Error, if you specify phitoo is must be positive");
00276 Dict d;
00277 d["type"] = "eman";
00278 d["az"] = az;
00279 d["alt"] = alt;
00280 d["phi"] = phi;
00281 Transform t(d);
00282 v.push_back(t);
00283 if ( phitoo != 0 ) {
00284 if (phitoo < 0) return false;
00285 else {
00286 for ( float p = phitoo; p <= 360.0f-phitoo; p+= phitoo )
00287 {
00288 d["phi"] = fmod(phi+p,360);
00289 Transform t(d);
00290 v.push_back(t);
00291 }
00292 }
00293 }
00294 return true;
00295 }
00296
00297 float EmanOrientationGenerator::get_az_delta(const float& delta,const float& altitude, const int) const
00298 {
00299
00300 float tmp = (float)(EMConsts::deg2rad * altitude);
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314 return altitude==0?360.0f:delta/sin(tmp);
00315
00316 }
00317
00318
00319 int EmanOrientationGenerator::get_orientations_tally(const Symmetry3D* const sym, const float& delta) const
00320 {
00321
00322
00323 bool inc_mirror = params.set_default("inc_mirror",false);
00324 bool breaksym = params.set_default("breaksym",false);
00325 Dict delimiters = sym->get_delimiters(inc_mirror);
00326 float altmax = delimiters["alt_max"];
00327 float azmax = delimiters["az_max"];
00328
00329 float paltmin = params.set_default("alt_min",0.0f);
00330 float paltmax = params.set_default("alt_max",180.0f);
00331 if (altmax>paltmax) altmax=paltmax;
00332
00333 float alt_iterator = 0.0f;
00334
00335
00336
00337 if (sym->is_h_sym()) alt_iterator = delimiters["alt_min"];
00338
00339 int tally = 0;
00340 while ( alt_iterator <= altmax ) {
00341 float h = get_az_delta(delta,alt_iterator, sym->get_max_csym() );
00342
00343
00344 if ( (alt_iterator > 0) && ( (azmax/h) < 2.8f) ) h = azmax / 2.1f;
00345 else if (alt_iterator == 0) h = azmax;
00346
00347 float az_iterator = 0.0f;
00348
00349 float azmax_adjusted = azmax;
00350 bool d_odd_mirror_flag = false;
00351 get_az_max(sym,altmax, inc_mirror,alt_iterator, h,d_odd_mirror_flag, azmax_adjusted);
00352 if (alt_iterator<paltmin) { alt_iterator += delta; continue; }
00353
00354 while ( az_iterator <= azmax_adjusted ) {
00355
00356
00357
00358
00359
00360
00361 if (sym->is_platonic_sym()) {
00362 if ( sym->is_in_asym_unit(alt_iterator, az_iterator,inc_mirror) == false ) {
00363 az_iterator += h;
00364 continue;
00365 }
00366 }
00367
00368 tally++;
00369 if ( sym->is_h_sym() && inc_mirror && alt_iterator != (float) delimiters["alt_min"] ) {
00370 tally++;
00371 }
00372 az_iterator += h;
00373 if ( (az_iterator > azmax_adjusted) && d_odd_mirror_flag) {
00374 azmax_adjusted = azmax;
00375 az_iterator += azmax/2.0f;
00376 }
00377 }
00378 alt_iterator += delta;
00379 }
00380
00381 if (breaksym) return tally*sym->get_nsym();
00382 return tally;
00383 }
00384
00385 vector<Transform> EmanOrientationGenerator::gen_orientations(const Symmetry3D* const sym) const
00386 {
00387 float delta = params.set_default("delta", 0.0f);
00388 int n = params.set_default("n", 0);
00389 bool breaksym = params.set_default("breaksym",false);
00390
00391 if ( delta <= 0 && n <= 0 ) throw InvalidParameterException("Error, you must specify a positive non-zero delta or n");
00392 if ( delta > 0 && n > 0 ) throw InvalidParameterException("Error, the delta and the n arguments are mutually exclusive");
00393
00394 if ( n > 0 ) {
00395 delta = get_optimal_delta(sym,n);
00396 }
00397
00398 bool inc_mirror = params.set_default("inc_mirror",false);
00399 bool inc_mirror_real = inc_mirror;
00400 if (breaksym) inc_mirror=true;
00401 Dict delimiters = sym->get_delimiters(inc_mirror);
00402 float altmax = delimiters["alt_max"];
00403 float azmax = delimiters["az_max"];
00404
00405 float paltmin = params.set_default("alt_min",0.0f);
00406 float paltmax = params.set_default("alt_max",180.0f);
00407 if (altmax>paltmax) altmax=paltmax;
00408
00409 bool perturb = params.set_default("perturb",true);
00410
00411 float alt_iterator = 0.0f;
00412
00413
00414
00415 if (sym->is_h_sym()) alt_iterator = delimiters["alt_min"];
00416
00417 vector<Transform> ret;
00418 while ( alt_iterator <= altmax ) {
00419 float h = get_az_delta(delta,alt_iterator, sym->get_max_csym() );
00420
00421
00422 if ( (alt_iterator > 0) && ( (azmax/h) < 2.8f) ) h = azmax / 2.1f;
00423 else if (alt_iterator == 0) h = azmax;
00424
00425 float az_iterator = 0.0f;
00426
00427 float azmax_adjusted = azmax;
00428
00429 bool d_odd_mirror_flag = false;
00430 get_az_max(sym,altmax, inc_mirror,alt_iterator, h,d_odd_mirror_flag, azmax_adjusted);
00431 if (alt_iterator<paltmin) { alt_iterator += delta; continue; }
00432
00433
00434 while ( az_iterator <= azmax_adjusted ) {
00435
00436
00437
00438
00439
00440
00441
00442
00443 if (alt_iterator == 0 && az_iterator > 0){
00444 az_iterator += h;
00445 continue;
00446 }
00447
00448
00449 float alt_soln = alt_iterator;
00450 float az_soln = az_iterator;
00451
00452 if (sym->is_platonic_sym()) {
00453 if ( sym->is_in_asym_unit(alt_soln, az_soln,inc_mirror) == false ) {
00454 az_iterator += h;
00455 continue;
00456 }
00457
00458 az_soln += sym->get_az_alignment_offset();
00459 }
00460
00461 if ( perturb && alt_soln != 0 ) {
00462 alt_soln += Util::get_gauss_rand(0.0f,.125f*delta);
00463 az_soln += Util::get_gauss_rand(0.0f,h*.125f);
00464 }
00465
00466 add_orientation(ret,az_soln,alt_soln);
00467
00468
00469
00470 if ( sym->is_h_sym() && inc_mirror && alt_iterator != (float) delimiters["alt_min"] ) {
00471 add_orientation(ret, az_soln,2.0f*(float)delimiters["alt_min"]-alt_soln);
00472 }
00473
00474 az_iterator += h;
00475 if ( (az_iterator > azmax_adjusted) && d_odd_mirror_flag) {
00476 azmax_adjusted = azmax;
00477 az_iterator += azmax/2.0f;
00478 }
00479 }
00480 alt_iterator += delta;
00481 }
00482
00483
00484
00485 if (breaksym) {
00486
00487 int nwithsym=ret.size();
00488 int nsym=sym->get_nsym();
00489 for (int j=1; j<nsym; j++) {
00490 Transform t=sym->get_sym(j);
00491 for (int i=0; i<nwithsym; i++) {
00492 ret.push_back(ret[i]*t);
00493 }
00494 }
00495
00496
00497 if (!inc_mirror_real) {
00498 vector<Transform> ret2;
00499 for (vector<Transform>::iterator t=ret.begin(); t!=ret.end(); ++t) {
00500 if ((*t)[2][2]>=0) ret2.push_back(*t);
00501
00502 }
00503 return ret2;
00504 }
00505 }
00506
00507 return ret;
00508 }
00509
00510 vector<Transform> RandomOrientationGenerator::gen_orientations(const Symmetry3D* const sym) const
00511 {
00512 int n = params.set_default("n", 0);
00513
00514 if ( n <= 0 ) throw InvalidParameterException("You must specify a positive, non zero n for the Random Orientation Generator");
00515
00516 bool phitoo = params.set_default("phitoo", false);
00517 bool inc_mirror = params.set_default("inc_mirror", false);
00518
00519 vector<Transform> ret;
00520
00521 int i = 0;
00522 Dict d("type","eman");
00523 while ( i < n ){
00524 float u1 = Util::get_frand(-1.0f,1.0f);
00525 float u2 = Util::get_frand(-1.0f,1.0f);
00526 float s = u1*u1 + u2*u2;
00527 if ( s > 1.0f ) continue;
00528 float alpha = 2.0f*sqrtf(1.0f-s);
00529 float x = alpha * u1;
00530 float y = alpha * u2;
00531 float z = 2.0f*s-1.0f;
00532
00533 float altitude = (float)EMConsts::rad2deg*acos(z);
00534 float azimuth = (float)EMConsts::rad2deg*atan2(y,x);
00535
00536 float phi = 0.0f;
00537 if ( phitoo ) phi = Util::get_frand(0.0f,359.9999f);
00538
00539 d["az"] = azimuth; d["phi"] = phi; d["alt"] = altitude;
00540 Transform t(d);
00541
00542 if ( !(sym->is_c_sym() && sym->get_nsym() == 1)) t = sym->reduce(t);
00543
00544 if ( !sym->is_in_asym_unit(altitude,azimuth,inc_mirror) ){
00545
00546
00547
00548 }
00549 ret.push_back(t);
00550 i++;
00551 }
00552 return ret;
00553 }
00554
00555 int EvenOrientationGenerator::get_orientations_tally(const Symmetry3D* const sym, const float& delta) const
00556 {
00557 bool inc_mirror = params.set_default("inc_mirror",false);
00558 Dict delimiters = sym->get_delimiters(inc_mirror);
00559 float altmax = delimiters["alt_max"];
00560 float azmax = delimiters["az_max"];
00561
00562 float altmin = 0.0f;
00563
00564
00565 if (sym->is_h_sym()) altmin = delimiters["alt_min"];
00566
00567 int tally = 0;
00568
00569 for (float alt = altmin; alt <= altmax; alt += delta) {
00570 float detaz;
00571 int lt;
00572 if ((0.0f == alt)||(180.0f == alt)) {
00573 detaz = 360.0f;
00574 lt = 1;
00575 } else {
00576 detaz = delta/(float)sin(alt*EMConsts::deg2rad);
00577 lt = int(azmax/detaz)-1;
00578 if (lt < 1) lt = 1;
00579 detaz = azmax/(float)lt;
00580 }
00581
00582
00583 for (int i = 0; i < lt; i++) {
00584 float az = (float)i*detaz;
00585 if (sym->is_platonic_sym()) {
00586 if ( sym->is_in_asym_unit(alt, az,inc_mirror) == false ) continue;
00587 }
00588 tally++;
00589 if ( sym->is_h_sym() && inc_mirror && alt != altmin ) {
00590 tally++;
00591 }
00592 }
00593 }
00594
00595 return tally;
00596 }
00597
00598 vector<Transform> EvenOrientationGenerator::gen_orientations(const Symmetry3D* const sym) const
00599 {
00600 float delta = params.set_default("delta", 0.0f);
00601 int n = params.set_default("n", 0);
00602
00603 if ( delta <= 0 && n <= 0 ) throw InvalidParameterException("Error, you must specify a positive non-zero delta or n");
00604 if ( delta > 0 && n > 0 ) throw InvalidParameterException("Error, the delta and the n arguments are mutually exzclusive");
00605
00606 if ( n > 0 ) {
00607 delta = get_optimal_delta(sym,n);
00608 }
00609
00610 bool inc_mirror = params.set_default("inc_mirror",false);
00611 Dict delimiters = sym->get_delimiters(inc_mirror);
00612 float altmax = delimiters["alt_max"];
00613 float azmax = delimiters["az_max"];
00614
00615 float altmin = 0.0f;
00616
00617
00618 if (sym->is_h_sym()) altmin = delimiters["alt_min"];
00619
00620 vector<Transform> ret;
00621
00622 for (float alt = altmin; alt <= altmax; alt += delta) {
00623 float detaz;
00624 int lt;
00625 if ((0.0f == alt)||(180.0f == alt)) {
00626 detaz = 360.0f;
00627 lt = 1;
00628 } else {
00629 detaz = delta/(float)sin(alt*EMConsts::deg2rad);
00630 lt = int(azmax/detaz)-1;
00631 if (lt < 1) lt = 1;
00632 detaz = azmax/(float)lt;
00633 }
00634
00635
00636
00637 for (int i = 0; i < lt; i++) {
00638 float az = (float)i*detaz;
00639 if (sym->is_platonic_sym()) {
00640 if ( sym->is_in_asym_unit(alt, az,inc_mirror) == false ) continue;
00641 else {
00642 az += sym->get_az_alignment_offset();
00643 }
00644 }
00645 add_orientation(ret,az,alt);
00646 if ( sym->is_h_sym() && inc_mirror && alt != altmin ) {
00647 add_orientation(ret,az,2.0f*altmin-alt);
00648 }
00649 }
00650 }
00651
00652 return ret;
00653 }
00654
00655 int SaffOrientationGenerator::get_orientations_tally(const Symmetry3D* const sym, const float& delta) const
00656 {
00657 bool inc_mirror = params.set_default("inc_mirror",false);
00658 Dict delimiters = sym->get_delimiters(inc_mirror);
00659 float altmax = delimiters["alt_max"];
00660 float azmax = delimiters["az_max"];
00661
00662 float altmin = 0.0f;
00663
00664
00665 if (sym->is_h_sym()){
00666 altmin = delimiters["alt_min"];
00667 if (inc_mirror) {
00668 altmin -= (float) sym->get_params()["maxtilt"];
00669 }
00670 }
00671
00672 float Deltaz = (float)(cos(altmax*EMConsts::deg2rad)-cos(altmin*EMConsts::deg2rad));
00673 float s = delta*M_PI/180.0f;
00674 float NFactor = 3.6f/s;
00675 float wedgeFactor = fabs( Deltaz*(azmax)/720.0f) ;
00676 int NumPoints = static_cast<int> (NFactor*NFactor*wedgeFactor);
00677
00678 int tally = 0;
00679 if (!sym->is_h_sym()) ++tally;
00680 float az = 0.0f;
00681 float dz = (float)cos(altmin*EMConsts::deg2rad);
00682 for(int i = 1; i < NumPoints; ++i ){
00683 float z = dz + Deltaz* (float)i/ float(NumPoints-1);
00684 float r= sqrt(1.0f-z*z);
00685 az = fmod(az + delta/r,azmax);
00686 float alt = (float)(acos(z)*EMConsts::rad2deg);
00687 if (sym->is_platonic_sym()) {
00688 if ( sym->is_in_asym_unit(alt,az,inc_mirror) == false ) continue;
00689 }
00690 tally++;
00691 }
00692
00693 return tally;
00694 }
00695
00696 vector<Transform> SaffOrientationGenerator::gen_orientations(const Symmetry3D* const sym) const
00697 {
00698 float delta = params.set_default("delta", 0.0f);
00699 int n = params.set_default("n", 0);
00700
00701 if ( delta <= 0 && n <= 0 ) throw InvalidParameterException("Error, you must specify a positive non-zero delta or n");
00702 if ( delta > 0 && n > 0 ) throw InvalidParameterException("Error, the delta and the n arguments are mutually exclusive");
00703
00704 if ( n > 0 ) {
00705 delta = get_optimal_delta(sym,n);
00706 }
00707
00708
00709
00710 bool inc_mirror = params.set_default("inc_mirror",false);
00711 Dict delimiters = sym->get_delimiters(inc_mirror);
00712 float altmax = delimiters["alt_max"];
00713 float azmax = delimiters["az_max"];
00714
00715 float altmin = 0.0f;
00716
00717
00718 if (sym->is_h_sym()){
00719 altmin = delimiters["alt_min"];
00720 if (inc_mirror) {
00721 altmin -= (float) sym->get_params()["maxtilt"];
00722 }
00723 }
00724
00725 float Deltaz = (float)(cos(altmax*EMConsts::deg2rad)-cos(altmin*EMConsts::deg2rad));
00726 float s = delta*M_PI/180.0f;
00727 float NFactor = 3.6f/s;
00728 float wedgeFactor = fabs( Deltaz*(azmax)/720.0f) ;
00729 int NumPoints = static_cast<int> (NFactor*NFactor*wedgeFactor);
00730
00731 vector<Transform> ret;
00732
00733 if (!sym->is_h_sym()) add_orientation(ret,0,0);
00734 float az = 0.0f;
00735 float dz = (float)cos(altmin*EMConsts::deg2rad);
00736 for(int i = 1; i < NumPoints; ++i ){
00737 float z = dz + Deltaz* (float)i/ float(NumPoints-1);
00738 float r= sqrt(1.0f-z*z);
00739 az = fmod(az + delta/r,azmax);
00740 float alt = (float)(acos(z)*EMConsts::rad2deg);
00741 if (sym->is_platonic_sym()) {
00742 if ( sym->is_in_asym_unit(alt,az,inc_mirror) == false ) continue;
00743 else {
00744 az += sym->get_az_alignment_offset();
00745 }
00746 }
00747 add_orientation(ret,az,alt);
00748 }
00749
00750 return ret;
00751 }
00752
00753 int OptimumOrientationGenerator::get_orientations_tally(const Symmetry3D* const sym, const float& delta) const
00754 {
00755 string deltaoptname = params.set_default("use","saff");
00756 Dict a;
00757 a["inc_mirror"] = (bool)params.set_default("inc_mirror",false);
00758 OrientationGenerator *g = Factory < OrientationGenerator >::get(deltaoptname,a);
00759 if (g) {
00760 int tally = g->get_orientations_tally(sym,delta);
00761 delete g;
00762 g = 0;
00763 return tally;
00764 }
00765 else throw;
00766 }
00767
00768 vector<Transform> OptimumOrientationGenerator::gen_orientations(const Symmetry3D* const sym) const
00769 {
00770 float delta = params.set_default("delta", 0.0f);
00771 int n = params.set_default("n", 0);
00772
00773 bool inc_mirror = params.set_default("inc_mirror",false);
00774
00775 if ( delta <= 0 && n <= 0 ) throw InvalidParameterException("Error, you must specify a positive non-zero delta or n");
00776 if ( delta > 0 && n > 0 ) throw InvalidParameterException("Error, the delta and the n arguments are mutually exclusive");
00777
00778 string generatorname = params.set_default("use","saff");
00779
00780 if ( n > 0 && generatorname != RandomOrientationGenerator::NAME ) {
00781 params["delta"] = get_optimal_delta(sym,n);
00782 params["n"] = (int)0;
00783 }
00784
00785
00786
00787
00788 params["inc_mirror"] = true;
00789 OrientationGenerator* g = Factory < OrientationGenerator >::get(generatorname);
00790 g->set_params(copy_relevant_params(g));
00791
00792
00793
00794 CSym* unit_sphere = new CSym();
00795 Dict nsym; nsym["nsym"] = 1; unit_sphere->set_params(nsym);
00796
00797 vector<Transform> unitsphereorientations = g->gen_orientations(unit_sphere);
00798 delete g; g = 0;
00799 delete unit_sphere; unit_sphere = 0;
00800
00801 vector<Vec3f> angles = optimize_distances(unitsphereorientations);
00802
00803 vector<Transform> ret;
00804 for (vector<Vec3f>::const_iterator it = angles.begin(); it != angles.end(); ++it ) {
00805 if ( sym->is_in_asym_unit((*it)[1],(*it)[0],inc_mirror) ) {
00806 add_orientation(ret,(*it)[0],(*it)[1]);
00807 }
00808 }
00809
00810
00811 params["inc_mirror"] = inc_mirror;
00812 params["delta"] = delta;
00813 params["n"] = n;
00814
00815 return ret;
00816 }
00817
00818 vector<Vec3f> OptimumOrientationGenerator::optimize_distances(const vector<Transform>& v) const
00819 {
00820 vector<Vec3f> points;
00821
00822 for (vector<Transform>::const_iterator it = v.begin(); it != v.end(); ++it ) {
00823 points.push_back(Vec3f(0,0,1)*(*it));
00824 }
00825
00826 if ( points.size() >= 2 ) {
00827 int max_it = 100;
00828 float percentage = 0.01f;
00829
00830 for ( int i = 0; i < max_it; ++i ){
00831 unsigned int p1 = 0;
00832 unsigned int p2 = 1;
00833
00834 float distsquared = (points[p1]-points[p2]).squared_length();
00835
00836
00837 for(unsigned int j = 0; j < points.size(); ++j) {
00838 for(unsigned int k = j+1; k < points.size(); ++k) {
00839 float d = (points[j]-points[k]).squared_length();
00840 if ( d < distsquared ) {
00841 distsquared = d;
00842 p1 = j;
00843 p2 = k;
00844 }
00845 }
00846 }
00847
00848
00849 Vec3f delta = percentage*(points[p2]-points[p1]);
00850
00851 points[p2] += delta;
00852 points[p2].normalize();
00853 points[p1] -= delta;
00854 points[p1].normalize();
00855 }
00856 }
00857
00858 vector<Vec3f> ret;
00859 for (vector<Vec3f>::const_iterator it = points.begin(); it != points.end(); ++it ) {
00860 float altitude = (float)(EMConsts::rad2deg*acos((*it)[2]));
00861 float azimuth = (float)(EMConsts::rad2deg*atan2((*it)[1],(*it)[0]));
00862 ret.push_back(Vec3f(90.0f+azimuth,altitude,0));
00863 }
00864
00865 return ret;
00866 }
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936 Symmetry3D::Symmetry3D() : cached_au_planes(0),cache_size(0),num_triangles(0),au_sym_triangles() {}
00937 Symmetry3D::~Symmetry3D() {
00938 if (cached_au_planes != 0 ) {
00939 delete_au_planes();
00940 }
00941 }
00942
00943 void verify(const Vec3f& tmp, float * plane, const string& message )
00944 {
00945 cout << message << " residual " << plane[0]*tmp[0]+plane[1]*tmp[1]+plane[2]*tmp[2] + plane[3] << endl;
00946 }
00947
00948 Transform Symmetry3D::reduce(const Transform& t, int n) const
00949 {
00950
00951 int soln = in_which_asym_unit(t);
00952
00953
00954 if ( soln == -1 ) {
00955 cout << "error, no solution found!" << endl;
00956
00957 return t;
00958 }
00959
00960
00961 Transform nt = get_sym(soln);
00962
00963 nt.invert();
00964
00965 nt = t*nt;
00966
00967 if ( n != 0 ) {
00968 nt = nt*get_sym(n);
00969 }
00970
00971 return nt;
00972
00973 }
00974
00975 int Symmetry3D::in_which_asym_unit(const Transform& t3d) const
00976 {
00977
00978
00979
00980 Vec3f p(0,0,1);
00981
00982 Transform o(t3d);
00983
00984
00985 o.invert();
00986
00987
00988 p = o*p;
00989
00990 return point_in_which_asym_unit(p);
00991 }
00992
00993
00994 void Symmetry3D::cache_au_planes() const {
00995 if (cached_au_planes == 0 ) {
00996 vector< vector<Vec3f> > au_triangles = get_asym_unit_triangles(true);
00997 num_triangles = au_triangles.size();
00998 cache_size = get_nsym()*au_triangles.size();
00999
01000 cached_au_planes = new float*[cache_size];
01001 float** fit = cached_au_planes;
01002 for(int i =0; i < cache_size; ++i,++fit) {
01003 float *t = new float[4];
01004 *fit = t;
01005 }
01006
01007
01008 int k = 0;
01009 for(int i = 0; i < get_nsym(); ++i) {
01010
01011 for( ncit it = au_triangles.begin(); it != au_triangles.end(); ++it, ++k)
01012 {
01013
01014 vector<Vec3f> points = *it;
01015 if ( i != 0 ) {
01016 for (vector<Vec3f>::iterator iit = points.begin(); iit != points.end(); ++iit ) {
01017
01018
01019 *iit = (*iit)*get_sym(i);
01020 }
01021 }
01022
01023 au_sym_triangles.push_back(points);
01024
01025
01026 Util::equation_of_plane(points[0],points[2],points[1],cached_au_planes[k]);
01027 }
01028 }
01029 }
01030 else throw UnexpectedBehaviorException("Attempt to generate a cache when cache exists");
01031 }
01032
01033 void Symmetry3D::delete_au_planes() {
01034 if (cached_au_planes == 0 ) throw UnexpectedBehaviorException("Attempt to delete a cache that does not exist");
01035 float** fit = cached_au_planes;
01036 for(int i =0; i < cache_size; ++i,++fit) {
01037 if (*fit == 0) throw UnexpectedBehaviorException("Attempt to delete a cache that does not exist");
01038 delete [] *fit;
01039 *fit = 0;
01040 }
01041
01042 delete [] cached_au_planes;
01043 cached_au_planes = 0;
01044 }
01045
01046 int Symmetry3D::point_in_which_asym_unit(const Vec3f& p) const
01047 {
01048 if (cached_au_planes == 0) {
01049 cache_au_planes();
01050 }
01051
01052 float epsNow=0.01f;
01053 int k = 0;
01054 for(int i = 0; i < get_nsym(); ++i) {
01055 for( int j = 0; j < num_triangles; ++j,++k) {
01056 vector<Vec3f> points = au_sym_triangles[k];
01057
01058 float* plane = cached_au_planes[k];
01059 Vec3f tmp = p;
01060
01061
01062 float scale = plane[0]*tmp[0]+plane[1]*tmp[1]+plane[2]*tmp[2];
01063 if ( scale != 0 )
01064 scale = -plane[3]/scale;
01065 else {
01066
01067 continue;
01068 }
01069
01070
01071 if (scale <= 0) continue;
01072
01073
01074 Vec3f pp = tmp*scale;
01075
01076
01077
01078
01079
01080 Vec3f v = points[2]-points[0];
01081 Vec3f u = points[1]-points[0];
01082 Vec3f w = pp - points[0];
01083
01084 float udotu = u.dot(u);
01085 float udotv = u.dot(v);
01086 float udotw = u.dot(w);
01087 float vdotv = v.dot(v);
01088 float vdotw = v.dot(w);
01089
01090 float d = 1.0f/(udotv*udotv - udotu*vdotv);
01091 float s = udotv*vdotw - vdotv*udotw;
01092 s *= d;
01093
01094 float t = udotv*udotw - udotu*vdotw;
01095 t *= d;
01096
01097
01098
01099 if (fabs(s) < Transform::ERR_LIMIT ) s = 0;
01100 if (fabs(t) < Transform::ERR_LIMIT ) t = 0;
01101
01102 if ( fabs((fabs(s)-1.0)) < Transform::ERR_LIMIT ) s = 1;
01103 if ( fabs((fabs(t)-1.0)) < Transform::ERR_LIMIT ) t = 1;
01104
01105
01106 if ( s >= -epsNow && t >= -epsNow && (s+t) <= 1+epsNow ) {
01107
01108 return i;
01109 }
01110 }
01111 }
01112
01113 return -1;
01114 }
01115 vector<Transform> Symmetry3D::get_touching_au_transforms(bool inc_mirror) const
01116 {
01117 vector<Transform> ret;
01118 vector<int> hit_cache;
01119
01120 vector<Vec3f> points = get_asym_unit_points(inc_mirror);
01121
01122
01123
01124
01125 if (inc_mirror && is_d_sym() && (get_nsym()/2 % 2 == 0)) {
01126 Dict delim = get_delimiters(false);
01127 float angle = (float)(EMConsts::deg2rad*float(delim["az_max"]));
01128 float y = -cos(angle);
01129 float x = sin(angle);
01130 points.push_back(Vec3f(x,y,0));
01131 }
01132 else if ( is_d_sym() && (get_nsym()/2 % 2 == 1)) {
01133 Dict delim = get_delimiters(false);
01134 float angle = float(delim["az_max"])/2.0f;
01135
01136 angle *= (float)EMConsts::deg2rad;
01137 float y = -cos(angle);
01138 float x = sin(angle);
01139 points.push_back(Vec3f(x,y,0));
01140
01141 if ( inc_mirror ) {
01142 angle = 3.0f*(float(delim["az_max"]))/2.0f;
01143 angle *= (float)EMConsts::deg2rad;
01144 float y = -cos(angle);
01145 float x = sin(angle);
01146 points.push_back(Vec3f(x,y,0));
01147 }
01148 }
01149
01150 typedef vector<Vec3f>::const_iterator const_point_it;
01151 for(const_point_it point = points.begin(); point != points.end(); ++point ) {
01152
01153 for(int i = 1; i < get_nsym(); ++i) {
01154
01155 if ( find(hit_cache.begin(),hit_cache.end(),i) != hit_cache.end() ) continue;
01156 Transform t = get_sym(i);
01157 Vec3f result = (*point)*t;
01158
01159 if (is_platonic_sym()) {
01160 for(const_point_it tmp = points.begin(); tmp != points.end(); ++tmp ) {
01161 Vec3f tt = result-(*tmp);
01162 if (tt.squared_length() < 0.01f) {
01163 hit_cache.push_back(i);
01164 ret.push_back(t);
01165 }
01166
01167 }
01168 }
01169 else {
01170 result -= *point;
01171 if (result.squared_length() < 0.05f) {
01172 hit_cache.push_back(i);
01173 ret.push_back(t);
01174 }
01175 }
01176 }
01177
01178 }
01179
01180 return ret;
01181 }
01182
01183
01184 vector<Transform> Symmetry3D::get_syms() const
01185 {
01186
01187 vector<Transform> ret;
01188
01189 for(int i = 0; i < get_nsym(); ++i ) {
01190 ret.push_back(get_sym(i));
01191 }
01192
01193
01194
01195
01196
01197 return ret;
01198 }
01199
01200 vector<Transform> Symmetry3D::get_symmetries(const string& symmetry)
01201 {
01202 Symmetry3D* sym = Factory<Symmetry3D>::get(Util::str_to_lower(symmetry));
01203 vector<Transform> ret = sym->get_syms();
01204 delete sym;
01205 return ret;
01206 }
01207
01208
01209 Dict CSym::get_delimiters(const bool inc_mirror) const {
01210 Dict returnDict;
01211
01212 int nsym = params.set_default("nsym",0);
01213 if ( nsym <= 0 ) throw InvalidValueException(nsym,"Error, you must specify a positive non zero nsym");
01214
01215 if ( inc_mirror ) returnDict["alt_max"] = 180.0f;
01216 else returnDict["alt_max"] = 90.0f;
01217
01218 returnDict["az_max"] = 360.0f/(float)nsym;
01219
01220 return returnDict;
01221 }
01222
01223 bool CSym::is_in_asym_unit(const float& altitude, const float& azimuth, const bool inc_mirror = false) const
01224 {
01225 Dict d = get_delimiters(inc_mirror);
01226 float alt_max = d["alt_max"];
01227 float az_max = d["az_max"];
01228
01229 int nsym = params.set_default("nsym",0);
01230 if ( nsym != 1 && azimuth < 0) return false;
01231 if ( altitude <= alt_max && azimuth <= az_max ) return true;
01232 return false;
01233 }
01234
01235 vector<vector<Vec3f> > CSym::get_asym_unit_triangles(bool inc_mirror ) const{
01236 vector<Vec3f> v = get_asym_unit_points(inc_mirror);
01237 int nsym = params.set_default("nsym",0);
01238
01239 vector<vector<Vec3f> > ret;
01240 if (v.size() == 0) return ret;
01241 if (nsym == 1 && !inc_mirror) {
01242 Vec3f z(0,0,1);
01243 vector<Vec3f> tmp;
01244 tmp.push_back(z);
01245 tmp.push_back(v[1]);
01246 tmp.push_back(v[0]);
01247 ret.push_back(tmp);
01248
01249 vector<Vec3f> tmp2;
01250 tmp2.push_back(z);
01251 tmp2.push_back(v[2]);
01252 tmp2.push_back(v[1]);
01253 ret.push_back(tmp2);
01254
01255 vector<Vec3f> tmp3;
01256 tmp3.push_back(z);
01257 tmp3.push_back(v[3]);
01258 tmp3.push_back(v[2]);
01259 ret.push_back(tmp3);
01260
01261 vector<Vec3f> tmp4;
01262 tmp4.push_back(z);
01263 tmp4.push_back(v[0]);
01264 tmp4.push_back(v[3]);
01265 ret.push_back(tmp4);
01266 }
01267 else if (nsym == 2 && inc_mirror) {
01268 Vec3f x(1,0,0);
01269 vector<Vec3f> tmp;
01270 tmp.push_back(v[1]);
01271 tmp.push_back(v[0]);
01272 tmp.push_back(x);
01273 ret.push_back(tmp);
01274
01275 vector<Vec3f> tmp2;
01276 tmp2.push_back(v[2]);
01277 tmp2.push_back(v[1]);
01278 tmp2.push_back(x);
01279 ret.push_back(tmp2);
01280
01281 vector<Vec3f> tmp3;
01282 tmp3.push_back(v[3]);
01283 tmp3.push_back(v[2]);
01284 tmp3.push_back(x);
01285 ret.push_back(tmp3);
01286
01287 vector<Vec3f> tmp4;
01288 tmp4.push_back(v[0]);
01289 tmp4.push_back(v[3]);
01290 tmp4.push_back(x);
01291 ret.push_back(tmp4);
01292 }
01293 else if (nsym == 2 && !inc_mirror) {
01294 vector<Vec3f> tmp;
01295 tmp.push_back(v[0]);
01296 tmp.push_back(v[2]);
01297 tmp.push_back(v[1]);
01298 ret.push_back(tmp);
01299
01300 vector<Vec3f> tmp2;
01301 tmp2.push_back(v[2]);
01302 tmp2.push_back(v[0]);
01303 tmp2.push_back(v[3]);
01304 ret.push_back(tmp2);
01305 }
01306 else if (v.size() == 3) {
01307 vector<Vec3f> tmp;
01308 tmp.push_back(v[0]);
01309 tmp.push_back(v[2]);
01310 tmp.push_back(v[1]);
01311 ret.push_back(tmp);
01312 }
01313 else if (v.size() == 4) {
01314 vector<Vec3f> tmp;
01315 tmp.push_back(v[0]);
01316 tmp.push_back(v[3]);
01317 tmp.push_back(v[1]);
01318 ret.push_back(tmp);
01319
01320 vector<Vec3f> tmp2;
01321 tmp2.push_back(v[1]);
01322 tmp2.push_back(v[3]);
01323 tmp2.push_back(v[2]);
01324 ret.push_back(tmp2);
01325 }
01326
01327 return ret;
01328 }
01329
01330 vector<Vec3f> CSym::get_asym_unit_points(bool inc_mirror) const
01331 {
01332 Dict delim = get_delimiters(inc_mirror);
01333 int nsym = params.set_default("nsym",0);
01334 vector<Vec3f> ret;
01335
01336 if ( nsym == 1 ) {
01337 if (inc_mirror == false ) {
01338 ret.push_back(Vec3f(0,-1,0));
01339 ret.push_back(Vec3f(1,0,0));
01340 ret.push_back(Vec3f(0,1,0));
01341 ret.push_back(Vec3f(-1,0,0));
01342 }
01343
01344 }
01345 else if (nsym == 2 && !inc_mirror) {
01346 ret.push_back(Vec3f(0,0,1));
01347 ret.push_back(Vec3f(0,-1,0));
01348 ret.push_back(Vec3f(1,0,0));
01349 ret.push_back(Vec3f(0,1,0));
01350 }
01351 else {
01352 ret.push_back(Vec3f(0,0,1));
01353 ret.push_back(Vec3f(0,-1,0));
01354 if (inc_mirror == true) {
01355 ret.push_back(Vec3f(0,0,-1));
01356 }
01357 float angle = (float)(EMConsts::deg2rad*float(delim["az_max"]));
01358 float y = -cos(angle);
01359 float x = sin(angle);
01360 ret.push_back(Vec3f(x,y,0));
01361 }
01362
01363 return ret;
01364
01365 }
01366
01367 Transform CSym::get_sym(const int n) const {
01368 int nsym = params.set_default("nsym",0);
01369 if ( nsym <= 0 ) throw InvalidValueException(n,"Error, you must specify a positive non zero nsym");
01370
01371 Dict d("type","eman");
01372
01373 d["az"] = (n%nsym) * 360.0f / nsym;
01374 d["alt"] = 0.0f;
01375 d["phi"] = 0.0f;
01376 return Transform(d);
01377 }
01378
01379
01380 Dict DSym::get_delimiters(const bool inc_mirror) const {
01381 Dict returnDict;
01382
01383
01384 int nsym = params.set_default("nsym",0);
01385 if ( nsym <= 0 ) throw InvalidValueException(nsym,"Error, you must specify a positive non zero nsym");
01386
01387 returnDict["alt_max"] = 90.0f;
01388
01389 if ( inc_mirror ) returnDict["az_max"] = 360.0f/(float)nsym;
01390 else returnDict["az_max"] = 180.0f/(float)nsym;
01391
01392 return returnDict;
01393 }
01394
01395 bool DSym::is_in_asym_unit(const float& altitude, const float& azimuth, const bool inc_mirror = false) const
01396 {
01397 Dict d = get_delimiters(inc_mirror);
01398 float alt_max = d["alt_max"];
01399 float az_max = d["az_max"];
01400
01401 int nsym = params.set_default("nsym",0);
01402
01403 if ( nsym == 1 && inc_mirror ) {
01404 if (altitude >= 0 && altitude <= alt_max && azimuth <= az_max ) return true;
01405 }
01406 else {
01407 if ( altitude >= 0 && altitude <= alt_max && azimuth <= az_max && azimuth >= 0 ) return true;
01408 }
01409 return false;
01410 }
01411
01412 Transform DSym::get_sym(const int n) const
01413 {
01414 int nsym = 2*params.set_default("nsym",0);
01415 if ( nsym <= 0 ) throw InvalidValueException(n,"Error, you must specify a positive non zero nsym");
01416
01417 Dict d("type","eman");
01418
01419 if (n >= nsym / 2) {
01420 d["az"] = ( (n%nsym) - nsym/2) * 360.0f / (nsym / 2);
01421 d["alt"] = 180.0f;
01422 d["phi"] = 0.0f;
01423 }
01424 else {
01425 d["az"] = (n%nsym) * 360.0f / (nsym / 2);
01426 d["alt"] = 0.0f;
01427 d["phi"] = 0.0f;
01428 }
01429 return Transform(d);
01430 }
01431
01432 vector<vector<Vec3f> > DSym::get_asym_unit_triangles(bool inc_mirror) const{
01433 vector<Vec3f> v = get_asym_unit_points(inc_mirror);
01434 int nsym = params.set_default("nsym",0);
01435 vector<vector<Vec3f> > ret;
01436 if ( (nsym == 1 && inc_mirror == false) || (nsym == 2 && inc_mirror)) {
01437 vector<Vec3f> tmp;
01438 tmp.push_back(v[0]);
01439 tmp.push_back(v[2]);
01440 tmp.push_back(v[1]);
01441 ret.push_back(tmp);
01442
01443 vector<Vec3f> tmp2;
01444 tmp2.push_back(v[2]);
01445 tmp2.push_back(v[0]);
01446 tmp2.push_back(v[3]);
01447 ret.push_back(tmp2);
01448 }
01449 else if (nsym == 1) {
01450 Vec3f z(0,0,1);
01451 vector<Vec3f> tmp;
01452 tmp.push_back(z);
01453 tmp.push_back(v[1]);
01454 tmp.push_back(v[0]);
01455 ret.push_back(tmp);
01456
01457 vector<Vec3f> tmp2;
01458 tmp2.push_back(z);
01459 tmp2.push_back(v[2]);
01460 tmp2.push_back(v[1]);
01461 ret.push_back(tmp2);
01462
01463 vector<Vec3f> tmp3;
01464 tmp3.push_back(z);
01465 tmp3.push_back(v[3]);
01466 tmp3.push_back(v[2]);
01467 ret.push_back(tmp3);
01468
01469 vector<Vec3f> tmp4;
01470 tmp4.push_back(z);
01471 tmp4.push_back(v[0]);
01472 tmp4.push_back(v[3]);
01473 ret.push_back(tmp4);
01474 }
01475 else {
01476
01477 vector<Vec3f> tmp;
01478 tmp.push_back(v[0]);
01479 tmp.push_back(v[2]);
01480 tmp.push_back(v[1]);
01481 ret.push_back(tmp);
01482 }
01483
01484 return ret;
01485 }
01486
01487 vector<Vec3f> DSym::get_asym_unit_points(bool inc_mirror) const
01488 {
01489 Dict delim = get_delimiters(inc_mirror);
01490
01491 vector<Vec3f> ret;
01492 int nsym = params.set_default("nsym",0);
01493 if ( nsym == 1 ) {
01494 if (inc_mirror == false ) {
01495 ret.push_back(Vec3f(0,0,1));
01496 ret.push_back(Vec3f(0,-1,0));
01497 ret.push_back(Vec3f(1,0,0));
01498 ret.push_back(Vec3f(0,1,0));
01499 }
01500 else {
01501 ret.push_back(Vec3f(0,-1,0));
01502 ret.push_back(Vec3f(1,0,0));
01503 ret.push_back(Vec3f(0,1,0));
01504 ret.push_back(Vec3f(-1,0,0));
01505 }
01506 }
01507 else if ( nsym == 2 && inc_mirror ) {
01508 ret.push_back(Vec3f(0,0,1));
01509 ret.push_back(Vec3f(0,-1,0));
01510 ret.push_back(Vec3f(1,0,0));
01511 ret.push_back(Vec3f(0,1,0));
01512 }
01513 else {
01514 float angle = (float)(EMConsts::deg2rad*float(delim["az_max"]));
01515 ret.push_back(Vec3f(0,0,1));
01516 ret.push_back(Vec3f(0,-1,0));
01517 float y = -cos(angle);
01518 float x = sin(angle);
01519 ret.push_back(Vec3f(x,y,0));
01520 }
01521
01522 return ret;
01523
01524 }
01525
01526
01527 Dict HSym::get_delimiters(const bool) const {
01528 Dict returnDict;
01529
01530
01531 int nsym = params.set_default("nsym",0);
01532 if ( nsym <= 0 ) throw InvalidValueException(nsym,"Error, you must specify a positive non zero nsym");
01533
01534 float maxtilt = params.set_default("maxtilt",5.0f);
01535
01536 returnDict["alt_max"] = 90.0f;
01537
01538 returnDict["alt_min"] = 90.0f - maxtilt;
01539
01540 returnDict["az_max"] = 360.0f;
01541
01542 return returnDict;
01543 }
01544
01545 bool HSym::is_in_asym_unit(const float& altitude, const float& azimuth, const bool inc_mirror = false) const
01546 {
01547 Dict d = get_delimiters(inc_mirror);
01548 float alt_max = d["alt_max"];
01549 float alt_min = d["alt_min"];
01550
01551 if (inc_mirror) {
01552 float e = params.set_default("maxtilt",5.0f);
01553 alt_min -= e;
01554 }
01555
01556 float az_max = d["az_max"];
01557
01558 if ( altitude >=alt_min && altitude <= alt_max && azimuth <= az_max && azimuth >= 0 ) return true;
01559 return false;
01560 }
01561
01562 vector<vector<Vec3f> > HSym::get_asym_unit_triangles(bool ) const{
01563
01564 vector<vector<Vec3f> > ret;
01565 return ret;
01566 }
01567
01568 vector<Vec3f> HSym::get_asym_unit_points(bool inc_mirror) const
01569 {
01570 vector<Vec3f> ret;
01571
01572 Dict delim = get_delimiters(inc_mirror);
01573 int nsym = params.set_default("nsym",1);
01574 float az = -(float)delim["az_max"];
01575
01576
01577 bool tracing_arcs = false;
01578
01579
01580 if ( !tracing_arcs) {
01581 Vec3f a(0,-1,0);
01582 ret.push_back(a);
01583
01584 if ( nsym > 2 ) {
01585 Dict d("type","eman");
01586 d["phi"] = 0.0f;
01587 d["alt"] = 0.0f;
01588 d["az"] = az;
01589 Vec3f b = Transform(d)*a;
01590 ret.push_back(b);
01591 }
01592 else
01593 {
01594 ret.push_back(Vec3f(1,0,0));
01595
01596 ret.push_back(Vec3f(0,1,0));
01597
01598 if ( nsym == 1 ) {
01599 ret.push_back(Vec3f(-1,0,0));
01600 ret.push_back(a);
01601 }
01602 }
01603 }
01604 return ret;
01605
01606 }
01607
01608 Transform HSym::get_sym(const int n) const
01609 {
01610 int nstart=params["nstart"];
01611
01612 float apix = params.set_default("apix",1.0f);
01613 float daz= params["daz"];
01614 float tz=params["tz"];
01615 float dz=tz/apix;
01616 Dict d("type","eman");
01617
01618
01619
01620
01621
01622 d["az"]=(n%nstart)*(360.0/nstart)+floor(float(n)/nstart)*daz;
01623 d["alt"] = 0.0f;
01624 d["phi"] = 0.0f;
01625 Transform ret(d);
01626 ret.set_trans(0,0,(n/nstart)*dz);
01627 return ret;
01628 }
01629
01630
01631 void PlatonicSym::init()
01632 {
01633
01634
01635
01636 float cap_sig = 2.0f*M_PI/ get_max_csym();
01637
01638 platonic_params["az_max"] = cap_sig;
01639
01640
01641
01642 float alpha = acos(1.0f/(sqrtf(3.0f)*tan(cap_sig/2.0f)));
01643
01644 platonic_params["alt_max"] = alpha;
01645
01646
01647 platonic_params["theta_c_on_two"] = 1.0f/2.0f*acos( cos(cap_sig)/(1.0f-cos(cap_sig)));
01648
01649 }
01650
01651
01652 Dict PlatonicSym::get_delimiters(const bool inc_mirror) const
01653 {
01654 Dict ret;
01655 ret["az_max"] = EMConsts::rad2deg * (float) platonic_params["az_max"];
01656
01657 if ( inc_mirror == false )
01658 if ( get_name() == IcosahedralSym::NAME || get_name() == OctahedralSym::NAME )
01659 ret["az_max"] = 0.5f*EMConsts::rad2deg * (float) platonic_params["az_max"];
01660
01661
01662
01663
01664 ret["alt_max"] = (float)(EMConsts::rad2deg * (float) platonic_params["alt_max"]);
01665 return ret;
01666 }
01667
01668
01669 bool PlatonicSym::is_in_asym_unit(const float& altitude, const float& azimuth, const bool inc_mirror) const
01670 {
01671 Dict d = get_delimiters(inc_mirror);
01672 float alt_max = d["alt_max"];
01673 float az_max = d["az_max"];
01674
01675 if ( altitude >= 0 && altitude <= alt_max && azimuth <= az_max && azimuth >= 0) {
01676
01677
01678 float tmpaz = (float)(EMConsts::deg2rad * azimuth);
01679
01680 float cap_sig = platonic_params["az_max"];
01681 float alt_max = platonic_params["alt_max"];
01682 if ( tmpaz > ( cap_sig/2.0f ) )tmpaz = cap_sig - tmpaz;
01683
01684 float lower_alt_bound = platonic_alt_lower_bound(tmpaz, alt_max );
01685
01686
01687 float tmpalt = (float)(EMConsts::deg2rad * altitude);
01688 if ( lower_alt_bound > tmpalt ) {
01689 if ( inc_mirror == false )
01690 {
01691 if ( cap_sig/2.0f < tmpaz ) return false;
01692 else return true;
01693 }
01694 else return true;
01695 }
01696 return false;
01697 }
01698 return false;
01699 }
01700
01701 float PlatonicSym::platonic_alt_lower_bound(const float& azimuth, const float& alpha) const
01702 {
01703 float cap_sig = platonic_params["az_max"];
01704 float theta_c_on_two = platonic_params["theta_c_on_two"];
01705
01706 float baldwin_lower_alt_bound = sin(cap_sig/2.0f-azimuth)/tan(theta_c_on_two);
01707 baldwin_lower_alt_bound += sin(azimuth)/tan(alpha);
01708 baldwin_lower_alt_bound *= 1/sin(cap_sig/2.0f);
01709 baldwin_lower_alt_bound = atan(1/baldwin_lower_alt_bound);
01710
01711 return baldwin_lower_alt_bound;
01712 }
01713
01714 vector<vector<Vec3f> > PlatonicSym::get_asym_unit_triangles(bool inc_mirror) const{
01715 vector<Vec3f> v = get_asym_unit_points(inc_mirror);
01716 vector<vector<Vec3f> > ret;
01717 if (v.size() == 3) {
01718 vector<Vec3f> tmp;
01719 tmp.push_back(v[0]);
01720 tmp.push_back(v[2]);
01721 tmp.push_back(v[1]);
01722 ret.push_back(tmp);
01723 }
01724 else {
01725 vector<Vec3f> tmp;
01726 tmp.push_back(v[0]);
01727 tmp.push_back(v[2]);
01728 tmp.push_back(v[1]);
01729 ret.push_back(tmp);
01730
01731 vector<Vec3f> tmp2;
01732 tmp2.push_back(v[0]);
01733 tmp2.push_back(v[3]);
01734 tmp2.push_back(v[2]);
01735 ret.push_back(tmp2);
01736 }
01737
01738 return ret;
01739 }
01740
01741 vector<Vec3f> PlatonicSym::get_asym_unit_points(bool inc_mirror) const
01742 {
01743 vector<Vec3f> ret;
01744
01745 Vec3f b = Vec3f(0,0,1);
01746 ret.push_back(b);
01747 float theta_c_on_two = (float)platonic_params["theta_c_on_two"];
01748 float theta_c = 2*theta_c_on_two;
01749
01750 Vec3f c_on_two = Vec3f(0,-sin(theta_c_on_two),cos(theta_c_on_two));
01751 Vec3f c = Vec3f(0,-sin(theta_c),cos(theta_c));
01752 ret.push_back(c_on_two);
01753
01754 float cap_sig = platonic_params["az_max"];
01755 Vec3f a = Vec3f(sin(theta_c)*sin(cap_sig),-sin(theta_c)*cos(cap_sig),cos(theta_c));
01756
01757 Vec3f f = a+b+c;
01758 f.normalize();
01759
01760 ret.push_back(f);
01761
01762 if ( inc_mirror ) {
01763 Vec3f a_on_two = Vec3f(sin(theta_c_on_two)*sin(cap_sig),-sin(theta_c_on_two)*cos(cap_sig),cos(theta_c_on_two));
01764 ret.push_back(a_on_two);
01765 }
01766
01767 if ( get_az_alignment_offset() != 0 ) {
01768 Dict d("type","eman");
01769 d["az"] = get_az_alignment_offset();
01770 d["phi"] = 0.0f;
01771 d["alt"] = 0.0f;
01772 Transform t(d);
01773 for (vector<Vec3f>::iterator it = ret.begin(); it != ret.end(); ++it )
01774 {
01775 *it = (*it)*t;
01776 }
01777 }
01778
01779 return ret;
01780
01781 }
01782
01783 float IcosahedralSym::get_az_alignment_offset() const { return 234.0; }
01784
01785 Transform IcosahedralSym::get_sym(const int n) const
01786 {
01787
01788 static double lvl0=0.;
01789 static double lvl1= 63.4349;
01790 static double lvl2=116.5651;
01791 static double lvl3=180.0;
01792
01793 static double ICOS[180] = {
01794 0,lvl0,0, 0,lvl0,288, 0,lvl0,216, 0,lvl0,144, 0,lvl0,72,
01795 0,lvl1,36, 0,lvl1,324, 0,lvl1,252, 0,lvl1,180, 0,lvl1,108,
01796 72,lvl1,36, 72,lvl1,324, 72,lvl1,252, 72,lvl1,180, 72,lvl1,108,
01797 144,lvl1,36, 144,lvl1,324, 144,lvl1,252, 144,lvl1,180, 144,lvl1,108,
01798 216,lvl1,36, 216,lvl1,324, 216,lvl1,252, 216,lvl1,180, 216,lvl1,108,
01799 288,lvl1,36, 288,lvl1,324, 288,lvl1,252, 288,lvl1,180, 288,lvl1,108,
01800 36,lvl2,0, 36,lvl2,288, 36,lvl2,216, 36,lvl2,144, 36,lvl2,72,
01801 108,lvl2,0, 108,lvl2,288, 108,lvl2,216, 108,lvl2,144, 108,lvl2,72,
01802 180,lvl2,0, 180,lvl2,288, 180,lvl2,216, 180,lvl2,144, 180,lvl2,72,
01803 252,lvl2,0, 252,lvl2,288, 252,lvl2,216, 252,lvl2,144, 252,lvl2,72,
01804 324,lvl2,0, 324,lvl2,288, 324,lvl2,216, 324,lvl2,144, 324,lvl2,72,
01805 0,lvl3,0, 0,lvl3,288, 0,lvl3,216, 0,lvl3,144, 0,lvl3,72
01806 };
01807
01808 int idx = n % 60;
01809 Dict d("type","eman");
01810
01811 if (get_az_alignment_offset() == 234.0) {
01812 d["az"] =(float)ICOS[idx * 3 ]+90;
01813 d["alt"] = (float)ICOS[idx * 3 + 1];
01814 d["phi"] = (float)ICOS[idx * 3 + 2]-90;
01815
01816 }
01817 else {
01818 d["az"] =(float)(float)ICOS[idx * 3 ];
01819 d["alt"] = (float)ICOS[idx * 3 + 1];
01820 d["phi"] = (float)ICOS[idx * 3 + 2];
01821
01822 }
01823
01824
01825
01826
01827
01828
01829 return Transform(d);
01830
01831 }
01832
01833 float Icosahedral2Sym::get_az_alignment_offset() const { return 234.0; }
01834
01835 Transform Icosahedral2Sym::get_sym(const int n) const
01836 {
01837 static float matrices[60*9] = {
01838 1, 0, 0, 0, 1, 0, 0, 0, 1,
01839 0.30902, -0.80902, 0.5, 0.80902, 0.5, 0.30902, -0.5, 0.30902, 0.80902,
01840 -0.80902, -0.5, 0.30902, 0.5, -0.30902, 0.80902, -0.30902, 0.80902, 0.5,
01841 -0.80902, 0.5, -0.30902, -0.5, -0.30902, 0.80902, 0.30902, 0.80902, 0.5,
01842 0.30902, 0.80902, -0.5, -0.80902, 0.5, 0.30902, 0.5, 0.30902, 0.80902,
01843 -1, 0, 0, 0, -1, 0, 0, 0, 1,
01844 -0.30902, -0.80902, 0.5, 0.80902, -0.5, -0.30902, 0.5, 0.30902, 0.80902,
01845 0.30902, 0.80902, 0.5, -0.80902, 0.5, -0.30902, -0.5, -0.30902, 0.80902,
01846 -0.30902, 0.80902, 0.5, -0.80902, -0.5, 0.30902, 0.5, -0.30902, 0.80902,
01847 -0.5, 0.30902, 0.80902, 0.30902, -0.80902, 0.5, 0.80902, 0.5, 0.30902,
01848 0.5, -0.30902, 0.80902, -0.30902, 0.80902, 0.5, -0.80902, -0.5, 0.30902,
01849 0.80902, 0.5, 0.30902, -0.5, 0.30902, 0.80902, 0.30902, -0.80902, 0.5,
01850 0.80902, -0.5, -0.30902, 0.5, 0.30902, 0.80902, -0.30902, -0.80902, 0.5,
01851 0.5, 0.30902, -0.80902, 0.30902, 0.80902, 0.5, 0.80902, -0.5, 0.30902,
01852 -0.5, -0.30902, -0.80902, -0.30902, -0.80902, 0.5, -0.80902, 0.5, 0.30902,
01853 -0.30902, -0.80902, -0.5, 0.80902, -0.5, 0.30902, -0.5, -0.30902, 0.80902,
01854 0.30902, -0.80902, -0.5, 0.80902, 0.5, -0.30902, 0.5, -0.30902, 0.80902,
01855 -0.30902, 0.80902, -0.5, -0.80902, -0.5, -0.30902, -0.5, 0.30902, 0.80902,
01856 0.5, -0.30902, -0.80902, -0.30902, 0.80902, -0.5, 0.80902, 0.5, 0.30902,
01857 -0.5, 0.30902, -0.80902, 0.30902, -0.80902, -0.5, -0.80902, -0.5, 0.30902,
01858 -0.80902, -0.5, -0.30902, 0.5, -0.30902, -0.80902, 0.30902, -0.80902, 0.5,
01859 0.80902, 0.5, -0.30902, -0.5, 0.30902, -0.80902, -0.30902, 0.80902, 0.5,
01860 0.80902, -0.5, 0.30902, 0.5, 0.30902, -0.80902, 0.30902, 0.80902, 0.5,
01861 -0.80902, 0.5, 0.30902, -0.5, -0.30902, -0.80902, -0.30902, -0.80902, 0.5,
01862 -0.5, -0.30902, 0.80902, -0.30902, -0.80902, -0.5, 0.80902, -0.5, 0.30902,
01863 0.5, 0.30902, 0.80902, 0.30902, 0.80902, -0.5, -0.80902, 0.5, 0.30902,
01864 0, 0, 1, 1, 0, 0, 0, 1, 0,
01865 0, 1, 0, 0, 0, 1, 1, 0, 0,
01866 0, -1, 0, 0, 0, 1, -1, 0, 0,
01867 0, 0, -1, -1, 0, 0, 0, 1, 0,
01868 0, -1, 0, 0, 0, -1, 1, 0, 0,
01869 0, 1, 0, 0, 0, -1, -1, 0, 0,
01870 -0.80902, -0.5, 0.30902, -0.5, 0.30902, -0.80902, 0.30902, -0.80902, -0.5,
01871 0.80902, -0.5, -0.30902, -0.5, -0.30902, -0.80902, 0.30902, 0.80902, -0.5,
01872 0.5, 0.30902, -0.80902, -0.30902, -0.80902, -0.5, -0.80902, 0.5, -0.30902,
01873 -0.30902, -0.80902, -0.5, -0.80902, 0.5, -0.30902, 0.5, 0.30902, -0.80902,
01874 -0.80902, 0.5, -0.30902, 0.5, 0.30902, -0.80902, -0.30902, -0.80902, -0.5,
01875 -0.5, -0.30902, -0.80902, 0.30902, 0.80902, -0.5, 0.80902, -0.5, -0.30902,
01876 -0.5, 0.30902, -0.80902, -0.30902, 0.80902, 0.5, 0.80902, 0.5, -0.30902,
01877 0, 0, -1, 1, 0, 0, 0, -1, 0,
01878 -0.80902, 0.5, 0.30902, 0.5, 0.30902, 0.80902, 0.30902, 0.80902, -0.5,
01879 0.80902, 0.5, -0.30902, 0.5, -0.30902, 0.80902, 0.30902, -0.80902, -0.5,
01880 -0.30902, 0.80902, -0.5, 0.80902, 0.5, 0.30902, 0.5, -0.30902, -0.80902,
01881 0.5, -0.30902, -0.80902, 0.30902, -0.80902, 0.5, -0.80902, -0.5, -0.30902,
01882 -0.80902, -0.5, -0.30902, -0.5, 0.30902, 0.80902, -0.30902, 0.80902, -0.5,
01883 -0.30902, -0.80902, 0.5, -0.80902, 0.5, 0.30902, -0.5, -0.30902, -0.80902,
01884 -0.30902, 0.80902, 0.5, 0.80902, 0.5, -0.30902, -0.5, 0.30902, -0.80902,
01885 1, 0, 0, 0, -1, 0, 0, 0, -1,
01886 0.30902, 0.80902, -0.5, 0.80902, -0.5, -0.30902, -0.5, -0.30902, -0.80902,
01887 0.30902, -0.80902, -0.5, -0.80902, -0.5, 0.30902, -0.5, 0.30902, -0.80902,
01888 -1, 0, 0, 0, 1, 0, 0, 0, -1,
01889 0.80902, 0.5, 0.30902, 0.5, -0.30902, -0.80902, -0.30902, 0.80902, -0.5,
01890 0.30902, -0.80902, 0.5, -0.80902, -0.5, -0.30902, 0.5, -0.30902, -0.80902,
01891 -0.5, 0.30902, 0.80902, -0.30902, 0.80902, -0.5, -0.80902, -0.5, -0.30902,
01892 0, 0, 1, -1, 0, 0, 0, -1, 0,
01893 0.5, -0.30902, 0.80902, 0.30902, -0.80902, -0.5, 0.80902, 0.5, -0.30902,
01894 0.30902, 0.80902, 0.5, 0.80902, -0.5, 0.30902, 0.5, 0.30902, -0.80902,
01895 0.80902, -0.5, 0.30902, -0.5, -0.30902, 0.80902, -0.30902, -0.80902, -0.5,
01896 -0.5, -0.30902, 0.80902, 0.30902, 0.80902, 0.5, -0.80902, 0.5, -0.30902,
01897 0.5, 0.30902, 0.80902, -0.30902, -0.80902, 0.5, 0.80902, -0.5, -0.30902
01898 };
01899
01900 int idx = n % 60;
01901
01902 std::vector<float> matrix(12, 0);
01903 for (int r = 0; r < 3; ++r) {
01904 for (int c = 0; c < 3; ++c) {
01905 matrix[r*4 + c] = matrices[idx*9 + r*3 + c];
01906 }
01907 }
01908
01909 Transform t3d(matrix);
01910 return t3d;
01911 }
01912
01913 Transform OctahedralSym::get_sym(const int n) const
01914 {
01915
01916
01917 static double lvl0=0.;
01918 static double lvl1=90.;
01919 static double lvl2=180.;
01920
01921 static double OCT[72] = {
01922 0,lvl0,0, 0,lvl0,90, 0,lvl0,180, 0,lvl0,270,
01923 0,lvl1,0, 0,lvl1,90, 0,lvl1,180, 0,lvl1,270,
01924 90,lvl1,0, 90,lvl1,90, 90,lvl1,180, 90,lvl1,270,
01925 180,lvl1,0, 180,lvl1,90, 180,lvl1,180, 180,lvl1,270,
01926 270,lvl1,0, 270,lvl1,90, 270,lvl1,180, 270,lvl1,270,
01927 0,lvl2,0, 0,lvl2,90, 0,lvl2,180, 0,lvl2,270
01928 };
01929
01930 int idx = n % 24;
01931
01932
01933 Dict d("type","eman");
01934 d["az"] = (float)OCT[idx * 3 ];
01935 d["alt"] = (float)OCT[idx * 3 + 1];
01936 d["phi"] = (float)OCT[idx * 3 + 2];
01937 return Transform(d);
01938
01939 }
01940
01941 float TetrahedralSym::get_az_alignment_offset() const { return 0.0; }
01942
01943 bool TetrahedralSym::is_in_asym_unit(const float& altitude, const float& azimuth, const bool inc_mirror) const
01944 {
01945 Dict d = get_delimiters(inc_mirror);
01946 float alt_max = d["alt_max"];
01947 float az_max = d["az_max"];
01948
01949 if ( altitude >= 0 && altitude <= alt_max && azimuth <= az_max && azimuth >= 0) {
01950
01951 float tmpaz = (float)(EMConsts::deg2rad * azimuth);
01952
01953 float cap_sig = platonic_params["az_max"];
01954 float alt_max = platonic_params["alt_max"];
01955 if ( tmpaz > ( cap_sig/2.0f ) )tmpaz = cap_sig - tmpaz;
01956
01957 float lower_alt_bound = platonic_alt_lower_bound(tmpaz, alt_max );
01958
01959
01960 float tmpalt = (float)(EMConsts::deg2rad * altitude);
01961 if ( lower_alt_bound > tmpalt ) {
01962 if ( !inc_mirror ) {
01963 float upper_alt_bound = platonic_alt_lower_bound( tmpaz, alt_max/2.0f);
01964
01965 if ( upper_alt_bound < tmpalt ) return false;
01966 else return true;
01967 }
01968 else return true;
01969 }
01970 return false;
01971 }
01972 else return false;
01973 }
01974
01975
01976 Transform TetrahedralSym::get_sym(const int n) const
01977 {
01978
01979
01980 static double lvl0=0;
01981 static double lvl1=109.4712;
01982
01983 static double TET[36] = {
01984 0,lvl0,0, 0,lvl0,120, 0,lvl0,240,
01985 0,lvl1,60, 0,lvl1,180, 0,lvl1,300,
01986 120,lvl1,60, 120,lvl1,180, 120,lvl1,300,
01987 240,lvl1,60, 240,lvl1,180, 240,lvl1,300
01988 };
01989
01990 int idx = n % 12;
01991
01992
01993
01994
01995 Dict d("type","eman");
01996 d["az"] = (float)TET[idx * 3 ];
01997 d["alt"] = (float)TET[idx * 3 + 1];
01998 d["phi"] = (float)TET[idx * 3 + 2];
01999 return Transform(d);
02000
02001 }
02002
02003
02004 vector<Vec3f> TetrahedralSym::get_asym_unit_points(bool inc_mirror) const
02005 {
02006 vector<Vec3f> ret;
02007
02008 Vec3f b = Vec3f(0,0,1);
02009 ret.push_back(b);
02010 float theta_c_on_two = (float)platonic_params["theta_c_on_two"];
02011 float theta_c = 2*theta_c_on_two;
02012
02013 Vec3f c_on_two = Vec3f(0,-sin(theta_c_on_two),cos(theta_c_on_two));
02014 Vec3f c = Vec3f(0,-sin(theta_c),cos(theta_c));
02015 ret.push_back(c_on_two);
02016 float cap_sig = platonic_params["az_max"];
02017 if ( inc_mirror ) {
02018 Vec3f a = Vec3f(sin(theta_c)*sin(cap_sig),-sin(theta_c)*cos(cap_sig),cos(theta_c));
02019
02020 Vec3f f = a+b+c;
02021 f.normalize();
02022
02023 ret.push_back(f);
02024 }
02025
02026 Vec3f a_on_two = Vec3f(sin(theta_c_on_two)*sin(cap_sig),-sin(theta_c_on_two)*cos(cap_sig),cos(theta_c_on_two));
02027 ret.push_back(a_on_two);
02028
02029
02030 if ( get_az_alignment_offset() != 0 ) {
02031 Dict d("type","eman");
02032 d["az"] = get_az_alignment_offset();
02033 d["phi"] = 0.0f;
02034 d["alt"] = 0.0f;
02035 Transform t(d);
02036 for (vector<Vec3f>::iterator it = ret.begin(); it != ret.end(); ++it )
02037 {
02038 *it = (*it)*t;
02039 }
02040 }
02041
02042 return ret;
02043 }
02044
02045
02046