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