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