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