Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

symmetry.cpp

Go to the documentation of this file.
00001 /*
00002  * Author: David Woolford, 09/23/2008 (woolford@bcm.edu)
00003  * Copyright (c) 2000-2006 Baylor College of Medicine
00004  *
00005  * This software is issued under a joint BSD/GNU license. You may use the
00006  * source code in this file under either license. However, note that the
00007  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00008  * so you are responsible for compliance with the licenses of these packages
00009  * if you opt to use BSD licensing. The warranty disclaimer below holds
00010  * in either instance.
00011  *
00012  * This complete copyright notice must be included in any revised version of the
00013  * source code. Additional authorship citations may be added, but existing
00014  * author citations must be preserved.
00015  *
00016  * This program is free software; you can redistribute it and/or modify
00017  * it under the terms of the GNU General Public License as published by
00018  * the Free Software Foundation; either version 2 of the License, or
00019  * (at your option) any later version.
00020  *
00021  * This program is distributed in the hope that it will be useful,
00022  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00023  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00024  * GNU General Public License for more details.
00025  *
00026  * You should have received a copy of the GNU General Public License
00027  * along with this program; if not, write to the Free Software
00028  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00029  *
00030  * */
00031 
00032 #include "symmetry.h"
00033 #include "transform.h"
00034 #include "vec3.h"
00035 #include "exception.h"
00036 #include "util.h"
00037 
00038 using namespace EMAN;
00039 
00040 const string CSym::NAME = "c";
00041 const string DSym::NAME = "d";
00042 const string HSym::NAME = "h";
00043 const string TetrahedralSym::NAME = "tet";
00044 const string OctahedralSym::NAME = "oct";
00045 const string IcosahedralSym::NAME = "icos";
00046 const string EmanOrientationGenerator::NAME = "eman";
00047 const string SaffOrientationGenerator::NAME = "saff";
00048 const string EvenOrientationGenerator::NAME = "even";
00049 const string RandomOrientationGenerator::NAME = "rand";
00050 const string OptimumOrientationGenerator::NAME = "opt";
00051 
00052 template <> Factory < Symmetry3D >::Factory()
00053 {
00054         force_add<CSym>();
00055         force_add<DSym>();
00056         force_add<HSym>();
00057         force_add<TetrahedralSym>();
00058         force_add<OctahedralSym>();
00059         force_add<IcosahedralSym>();
00060 }
00061 
00062 void EMAN::dump_symmetries()
00063 {
00064         dump_factory < Symmetry3D > ();
00065 }
00066 
00067 map<string, vector<string> > EMAN::dump_symmetries_list()
00068 {
00069         return dump_factory_list < Symmetry3D > ();
00070 }
00071 
00072 template <>
00073 Symmetry3D* Factory < Symmetry3D >::get(const string & instancename_)
00074 {
00075         init();
00076 
00077         const string instancename = Util::str_to_lower(instancename_);
00078 
00079         unsigned int n = instancename.size();
00080         if ( n == 0 ) throw NotExistingObjectException(instancename, "Empty instance name!");
00081 
00082         char leadingchar = instancename[0];
00083         if (leadingchar == 'c' || leadingchar == 'd' || leadingchar == 'h' ) {
00084                 Dict parms;
00085                 if (n > 1) {
00086                         int nsym = atoi(instancename.c_str() + 1);
00087                         parms["nsym"] = nsym;
00088                 }
00089 
00090                 if (leadingchar == 'c') {
00091                         return get("c",parms);
00092                 }
00093                 if (leadingchar == 'd') {
00094                         return get("d",parms);
00095                 }
00096                 if (leadingchar == 'h') {
00097                         return get("h",parms);
00098                 }
00099 
00100 //              delete lc;
00101         }
00102         else if ( instancename == "icos" || instancename == "oct" || instancename == "tet" )
00103         {
00104                 map < string, InstanceType >::iterator fi =
00105                                 my_instance->my_dict.find(instancename);
00106                 if (fi != my_instance->my_dict.end()) {
00107                         return my_instance->my_dict[instancename] ();
00108                 }
00109 
00110                 string lower = instancename;
00111                 for (unsigned int i=0; i<lower.length(); i++) lower[i]=tolower(lower[i]);
00112 
00113                 fi = my_instance->my_dict.find(lower);
00114                 if (fi != my_instance->my_dict.end()) {
00115                         return my_instance->my_dict[lower] ();
00116                 }
00117 
00118                 throw NotExistingObjectException(instancename, "No such an instance existing");
00119         }
00120         else throw NotExistingObjectException(instancename, "No such an instance existing");
00121 
00122         throw NotExistingObjectException(instancename, "No such an instance existing");
00123 }
00124 
00125 template <> Factory < OrientationGenerator >::Factory()
00126 {
00127         force_add<EmanOrientationGenerator>();
00128         force_add<RandomOrientationGenerator>();
00129         force_add<EvenOrientationGenerator>();
00130         force_add<SaffOrientationGenerator>();
00131         force_add<OptimumOrientationGenerator>();
00132 }
00133 
00134 
00135 
00136 void EMAN::dump_orientgens()
00137 {
00138         dump_factory < OrientationGenerator > ();
00139 }
00140 
00141 map<string, vector<string> > EMAN::dump_orientgens_list()
00142 {
00143         return dump_factory_list < OrientationGenerator > ();
00144 }
00145 
00146 vector<Transform> Symmetry3D::gen_orientations(const string& generatorname, const Dict& parms)
00147 {
00148         ENTERFUNC;
00149         vector<Transform> ret;
00150         OrientationGenerator *g = Factory < OrientationGenerator >::get(Util::str_to_lower(generatorname), parms);
00151         if (g) {
00152                 ret = g->gen_orientations(this);
00153                 if( g )
00154                 {
00155                         delete g;
00156                         g = 0;
00157                 }
00158         }
00159         else throw;
00160 
00161         EXITFUNC;
00162 
00163         return ret;
00164 }
00165 
00166 void OrientationGenerator::get_az_max(const Symmetry3D* const sym, const float& altmax, const bool inc_mirror, const float& alt_iterator,const float& h,bool& d_odd_mirror_flag, float& azmax_adjusted) const
00167 {
00168 
00169         if ( sym->is_d_sym() && alt_iterator == altmax && ( (sym->get_nsym())/2 % 2 == 1 )) {
00170                 if (inc_mirror) {
00171                         azmax_adjusted /= 4.0;
00172                         d_odd_mirror_flag = true;
00173                 }
00174                 else azmax_adjusted /= 2.0;
00175         }
00176         else if (sym->is_d_sym() && alt_iterator == altmax && ( (sym->get_nsym())/2 % 2 == 0 ) && inc_mirror) {
00177                 azmax_adjusted /= 2.0;
00178         }
00179         // if this is odd c symmetry, and we're at the equator, and we're excluding the mirror then
00180         // half the equator is redundant (it is the mirror of the other half)
00181         else if (sym->is_c_sym() && !inc_mirror && alt_iterator == altmax && (sym->get_nsym() % 2 == 1 ) ){
00182                 azmax_adjusted /= 2.0;
00183         }
00184         // at the azimuthal boundary in c symmetry and tetrahedral symmetry we have come
00185         // full circle, we must not include it
00186         else if (sym->is_c_sym() || sym->is_tet_sym() ) {
00187                 azmax_adjusted -=  h/4.0f;
00188         }
00189         // If we're including the mirror then in d and icos and oct symmetry the azimuthal
00190         // boundary represents coming full circle, so must be careful to exclude it
00191         else if (inc_mirror && ( sym->is_d_sym() || sym->is_platonic_sym() ) )  {
00192                 azmax_adjusted -=  h/4.0f;
00193         }
00194         // else do nothing - this means that we're including the great arc traversing
00195         // the full range of permissable altitude angles at azmax.
00196         // This happens in d symmetry, and in the icos and oct symmetries, when the mirror
00197         // portion of the asymmetric unit is being excluded
00198 
00199 }
00200 
00201 
00202 //vector<Vec3f> OrientationGenerator::get_graph_opt(const Symmetry3D* const sym) const {
00203 //      bool inc_mirror = params.set_default("inc_mirror",false);
00204 //      vector<Vec3f> seeds = get_asym_unit_points(inc_mirror);
00205 //      vector<Vec3i> connections;
00206 //      if (seeds.size() == 3) {
00207 //              connections.push_back(Vec2i(0,1,2));
00208 //      } else {
00209 //              throw;
00210 //
00211 //      }
00212 //}
00213 
00214 
00215 float OrientationGenerator::get_optimal_delta(const Symmetry3D* const sym, const int& n) const
00216 {
00217 
00218 //      float delta_soln = 360.0f/sym->get_max_csym();
00219         float delta_soln = 180.0;
00220         float delta_upper_bound = delta_soln;
00221         float delta_lower_bound = 0.0;
00222 
00223         int prev_tally = -1;
00224         // This is an example of a divide and conquer approach, the possible values of delta are searched
00225         // like a binary tree
00226 
00227         bool soln_found = false;
00228 
00229         while ( soln_found == false ) {
00230                 int tally = get_orientations_tally(sym,delta_soln);
00231                 if ( tally == n ) soln_found = true;
00232                 else if ( (delta_upper_bound - delta_lower_bound) < 0.0001 ) {
00233                         // If this is the case, the requested number of projections is practically infeasible
00234                         // in which case just return the nearest guess
00235                         soln_found = true;
00236                         delta_soln = (delta_upper_bound+delta_lower_bound)/2.0f;
00237                 }
00238                 else if (tally < n) {
00239                         delta_upper_bound = delta_soln;
00240                         delta_soln = delta_soln - (delta_soln-delta_lower_bound)/2.0f;
00241                 }
00242                 else  /* tally > n*/{
00243                         delta_lower_bound = delta_soln;
00244                         delta_soln = delta_soln  + (delta_upper_bound-delta_soln)/2.0f;
00245                 }
00246                 prev_tally = tally;
00247         }
00248 
00249         return delta_soln;
00250 }
00251 
00252 bool OrientationGenerator::add_orientation(vector<Transform>& v, const float& az, const float& alt) const
00253 {
00254         bool randphi = params.set_default("random_phi",false);
00255         float phi = 0.0;
00256         if (randphi) phi = Util::get_frand(0.0f,359.99999f);
00257         float phitoo = params.set_default("phitoo",0.0f);
00258         if ( phitoo < 0 ) throw InvalidValueException(phitoo, "Error, if you specify phitoo is must be positive");
00259         Dict d;
00260         d["type"] = "eman";
00261         d["az"] = az;
00262         d["alt"] = alt;
00263         d["phi"] = phi;
00264         Transform t(d);
00265         v.push_back(t);
00266         if ( phitoo != 0 ) {
00267                 if (phitoo < 0) return false;
00268                 else {
00269                         for ( float p = phitoo; p <= 360.0-phitoo; p+= phitoo )
00270                         {
00271                                 d["phi"] = fmod(phi+p,360);
00272                                 Transform t(d);
00273                                 v.push_back(t);
00274                         }
00275                 }
00276         }
00277         return true;
00278 }
00279 
00280 float EmanOrientationGenerator::get_az_delta(const float& delta,const float& altitude, const int maxcsym) const
00281 {
00282         // convert altitude into radians
00283         float tmp = (float)(EMConsts::deg2rad * altitude);
00284 
00285         // This is taken from EMAN1 project3d.C
00286         // This wasn't working like it was supposed to. Rather than
00287         // figuring it out, I'm just replacing it --steve
00288 /*      float h=floor(360.0f/(delta*1.1547f));  // the 1.1547 makes the overall distribution more like a hexagonal mesh
00289         h=(int)floor(h*sin(tmp)+.5f);
00290         if (h==0) h=1;
00291         h=abs(maxcsym)*floor(h/(float)abs(maxcsym)+.5f);
00292         if ( h == 0 ) h = (float)maxcsym;
00293         h=2.0f*M_PI/h;
00294 
00295         return (float)(EMConsts::rad2deg*h);*/
00296 
00297         return altitude==0?360.0f:delta/sin(tmp);
00298 
00299 }
00300 
00301 
00302 int EmanOrientationGenerator::get_orientations_tally(const Symmetry3D* const sym, const float& delta) const
00303 {
00304         //FIXME THIS IS SO SIMILAR TO THE gen_orientations function that they should be probably use
00305         // a common routine - SAME ISSUE FOR OTHER ORIENTATION GENERATORS
00306         bool inc_mirror = params.set_default("inc_mirror",false);
00307         Dict delimiters = sym->get_delimiters(inc_mirror);
00308         float altmax = delimiters["alt_max"];
00309         float azmax = delimiters["az_max"];
00310 
00311         float alt_iterator = 0.0f;
00312 
00313         // #If it's a h symmetry then the alt iterator starts at very close
00314         // #to the altmax... the object is a h symmetry then it knows its alt_min...
00315         if (sym->is_h_sym()) alt_iterator = delimiters["alt_min"];
00316 
00317         int tally = 0;
00318         while ( alt_iterator <= altmax ) {
00319                 float h = get_az_delta(delta,alt_iterator, sym->get_max_csym() );
00320 
00321                 // not sure what this does code taken from EMAN1 - FIXME original author add comments
00322                 if ( (alt_iterator > 0) && ( (azmax/h) < 2.8) ) h = azmax / 2.1f;
00323                 else if (alt_iterator == 0) h = azmax;
00324 
00325                 float az_iterator = 0.0;
00326 
00327                 float azmax_adjusted = azmax;
00328                 bool d_odd_mirror_flag = false;
00329                 get_az_max(sym,altmax, inc_mirror,alt_iterator, h,d_odd_mirror_flag, azmax_adjusted);
00330 
00331                 while ( az_iterator <= azmax_adjusted ) {
00332                         // FIXME: add an intelligent comment - this was copied from old code
00333 //                      if ( az_iterator > 180.0 && alt_iterator > 180.0/(2.0-0.001) && alt_iterator < 180.0/(2.0+0.001) ) {
00334 //                              az_iterator +=  h;
00335 //                              continue;
00336 //                      }
00337 
00338                         if (sym->is_platonic_sym()) {
00339                                 if ( sym->is_in_asym_unit(alt_iterator, az_iterator,inc_mirror) == false ) {
00340                                         az_iterator += h;
00341                                         continue;
00342                                 }
00343                         }
00344 
00345                         tally++;
00346                         if ( sym->is_h_sym() && inc_mirror && alt_iterator != (float) delimiters["alt_min"] ) {
00347                                 tally++;
00348                         }
00349                         az_iterator += h;
00350                         if ( (az_iterator > azmax_adjusted) && d_odd_mirror_flag) {
00351                                 azmax_adjusted = azmax;
00352                                 az_iterator += azmax/2.0f;
00353                         }
00354                 }
00355                 alt_iterator += delta;
00356         }
00357         return tally;
00358 }
00359 
00360 vector<Transform> EmanOrientationGenerator::gen_orientations(const Symmetry3D* const sym) const
00361 {
00362         float delta = params.set_default("delta", 0.0f);
00363         int n = params.set_default("n", 0);
00364 
00365         if ( delta <= 0 && n <= 0 ) throw InvalidParameterException("Error, you must specify a positive non-zero delta or n");
00366         if ( delta > 0 && n > 0 ) throw InvalidParameterException("Error, the delta and the n arguments are mutually exclusive");
00367 
00368         if ( n > 0 ) {
00369                 delta = get_optimal_delta(sym,n);
00370         }
00371 
00372         bool inc_mirror = params.set_default("inc_mirror",false);
00373         Dict delimiters = sym->get_delimiters(inc_mirror);
00374         float altmax = delimiters["alt_max"];
00375         float azmax = delimiters["az_max"];
00376 
00377         bool perturb = params.set_default("perturb",false);
00378 
00379         float alt_iterator = 0.0;
00380 
00381         // #If it's a h symmetry then the alt iterator starts at very close
00382         // #to the altmax... the object is a h symmetry then it knows its alt_min...
00383         if (sym->is_h_sym()) alt_iterator = delimiters["alt_min"];
00384 
00385         vector<Transform> ret;
00386         while ( alt_iterator <= altmax ) {
00387                 float h = get_az_delta(delta,alt_iterator, sym->get_max_csym() );
00388 
00389                 // not sure what this does code taken from EMAN1 - FIXME original author add comments
00390                 if ( (alt_iterator > 0) && ( (azmax/h) < 2.8) ) h = azmax / 2.1f;
00391                 else if (alt_iterator == 0) h = azmax;
00392 
00393                 float az_iterator = 0.0;
00394 
00395                 float azmax_adjusted = azmax;
00396 
00397                 bool d_odd_mirror_flag = false;
00398                 get_az_max(sym,altmax, inc_mirror,alt_iterator, h,d_odd_mirror_flag, azmax_adjusted);
00399 
00400 
00401                 while ( az_iterator <= azmax_adjusted ) {
00402                         // FIXME: add an intelligent comment - this was copied from old code
00403 //                      if ( az_iterator > 180.0 && alt_iterator > 180.0/(2.0-0.001) && alt_iterator < 180.0/(2.0+0.001) ) {
00404 //                              az_iterator +=  h;
00405 //                              continue;
00406 //                      }
00407 //                      // Now that I am handling the boundaries very specifically, I don't think we need
00408                         // the above if statement. But I am leaving it there in case I need to reconsider.
00409 
00410                         if (alt_iterator == 0 && az_iterator > 0){
00411                                 az_iterator += h;
00412                                 continue; // We only ever need to generate on orientation at alt=0
00413                         }
00414 
00415 
00416                         float alt_soln = alt_iterator;
00417                         float az_soln = az_iterator;
00418 
00419                         if (sym->is_platonic_sym()) {
00420                                 if ( sym->is_in_asym_unit(alt_soln, az_soln,inc_mirror) == false ) {
00421                                         az_iterator += h;
00422                                         continue;
00423                                 }
00424                                 // Some objects have alignment offsets (icos and tet particularly)
00425                                 az_soln += sym->get_az_alignment_offset();
00426                         }
00427 
00428                         if ( perturb &&  alt_soln != 0 ) {
00429                                 alt_soln += Util::get_gauss_rand(0.0f,.25f*delta);
00430                                 az_soln += Util::get_gauss_rand(0.0f,h/4.0f);
00431                         }
00432 
00433                         add_orientation(ret,az_soln,alt_soln);
00434 
00435                         // Add helical symmetry orientations on the other side of the equator (if we're including
00436                         // mirror orientations)
00437                         if ( sym->is_h_sym() && inc_mirror && alt_iterator != (float) delimiters["alt_min"] ) {
00438                                 add_orientation(ret, az_soln,2.0f*(float)delimiters["alt_min"]-alt_soln);
00439                         }
00440 
00441                         az_iterator += h;
00442                         if ( (az_iterator > azmax_adjusted) && d_odd_mirror_flag) {
00443                                 azmax_adjusted = azmax;
00444                                 az_iterator += azmax/2.0f;
00445                         }
00446                 }
00447                 alt_iterator += delta;
00448         }
00449 
00450         return ret;
00451 }
00452 
00453 vector<Transform> RandomOrientationGenerator::gen_orientations(const Symmetry3D* const sym) const
00454 {
00455         int n = params.set_default("n", 0);
00456 
00457         if ( n <= 0 ) throw InvalidParameterException("You must specify a positive, non zero n for the Random Orientation Generator");
00458 
00459         bool phitoo = params.set_default("phitoo", false);
00460         bool inc_mirror = params.set_default("inc_mirror", false);
00461 
00462         vector<Transform> ret;
00463 
00464         int i = 0;
00465         Dict d("type","eman");
00466         while ( i < n ){
00467                 float u1 =  Util::get_frand(-1.0f,1.0f);
00468                 float u2 =  Util::get_frand(-1.0f,1.0f);
00469                 float s = u1*u1 + u2*u2;
00470                 if ( s > 1.0 ) continue;
00471                 float alpha = 2.0f*sqrtf(1.0f-s);
00472                 float x = alpha * u1;
00473                 float y = alpha * u2;
00474                 float z = 2.0f*s-1.0f;
00475 
00476                 float altitude = (float)EMConsts::rad2deg*acos(z);
00477                 float azimuth = (float)EMConsts::rad2deg*atan2(y,x);
00478 
00479                 float phi = 0.0;
00480                 if ( phitoo ) phi = Util::get_frand(0.0f,359.9999f);
00481 
00482                 d["az"] = azimuth; d["phi"] = phi; d["alt"] = altitude;
00483                 Transform t(d);
00484 
00485                 if ( !(sym->is_c_sym() && sym->get_nsym() == 1)) t = sym->reduce(t); //reduce doesn't make sense for C1 symmetry
00486 
00487                 if ( !sym->is_in_asym_unit(altitude,azimuth,inc_mirror) ){
00488                         // is_in_asym_unit has returned the wrong value!
00489                         // FIXME
00490 //                      cout << "warning, there is an unresolved issue - email D Woolford" << endl;
00491                 }
00492                 ret.push_back(t);
00493                 i++;
00494         }
00495         return ret;
00496 }
00497 
00498 int EvenOrientationGenerator::get_orientations_tally(const Symmetry3D* const sym, const float& delta) const
00499 {
00500         bool inc_mirror = params.set_default("inc_mirror",false);
00501         Dict delimiters = sym->get_delimiters(inc_mirror);
00502         float altmax = delimiters["alt_max"];
00503         float azmax = delimiters["az_max"];
00504 
00505         float altmin = 0.0;
00506         // #If it's a h symmetry then the alt iterator starts at very close
00507         // #to the altmax... the object is a h symmetry then it knows its alt_min...
00508         if (sym->is_h_sym()) altmin = delimiters["alt_min"];
00509 
00510         int tally = 0;
00511 
00512         for (float alt = altmin; alt <= altmax; alt += delta) {
00513                 float detaz;
00514                 int lt;
00515                 if ((0.0 == alt)||(180.0 == alt)) {
00516                         detaz = 360.0f;
00517                         lt = 1;
00518                 } else {
00519                         detaz = delta/(float)sin(alt*EMConsts::deg2rad);
00520                         lt = int(azmax/detaz)-1;
00521                         if (lt < 1) lt = 1;
00522                         detaz = azmax/(float)lt;
00523                 }
00524 //              bool d_odd_mirror_flag = false;
00525 //              get_az_max(sym,altmax, inc_mirror,alt, lt,d_odd_mirror_flag, detaz);
00526                 for (int i = 0; i < lt; i++) {
00527                         float az = (float)i*detaz;
00528                         if (sym->is_platonic_sym()) {
00529                                 if ( sym->is_in_asym_unit(alt, az,inc_mirror) == false ) continue;
00530                         }
00531                         tally++;
00532                         if ( sym->is_h_sym() && inc_mirror && alt != altmin ) {
00533                                 tally++;
00534                         }
00535                 }
00536         }
00537 
00538         return tally;
00539 }
00540 
00541 vector<Transform> EvenOrientationGenerator::gen_orientations(const Symmetry3D* const sym) const
00542 {
00543         float delta = params.set_default("delta", 0.0f);
00544         int n = params.set_default("n", 0);
00545 
00546         if ( delta <= 0 && n <= 0 ) throw InvalidParameterException("Error, you must specify a positive non-zero delta or n");
00547         if ( delta > 0 && n > 0 ) throw InvalidParameterException("Error, the delta and the n arguments are mutually exzclusive");
00548 
00549         if ( n > 0 ) {
00550                 delta = get_optimal_delta(sym,n);
00551         }
00552 
00553         bool inc_mirror = params.set_default("inc_mirror",false);
00554         Dict delimiters = sym->get_delimiters(inc_mirror);
00555         float altmax = delimiters["alt_max"];
00556         float azmax = delimiters["az_max"];
00557 
00558         float altmin = 0.0;
00559         // If it's a h symmetry then the alt iterator starts at very close
00560         // to the altmax... the object is a h symmetry then it knows its alt_min...
00561         if (sym->is_h_sym()) altmin = delimiters["alt_min"];
00562 
00563         vector<Transform> ret;
00564 
00565         for (float alt = altmin; alt <= altmax; alt += delta) {
00566                 float detaz;
00567                 int lt;
00568                 if ((0.0 == alt)||(180.0 == alt)) {
00569                         detaz = 360.0f;
00570                         lt = 1;
00571                 } else {
00572                         detaz = delta/(float)sin(alt*EMConsts::deg2rad);
00573                         lt = int(azmax/detaz)-1;
00574                         if (lt < 1) lt = 1;
00575                         detaz = azmax/(float)lt;
00576                 }
00577 //              bool d_odd_mirror_flag = false;
00578 //              get_az_max(sym,altmax, inc_mirror,alt, lt,d_odd_mirror_flag, detaz);
00579 
00580                 for (int i = 0; i < lt; i++) {
00581                         float az = (float)i*detaz;
00582                         if (sym->is_platonic_sym()) {
00583                                 if ( sym->is_in_asym_unit(alt, az,inc_mirror) == false ) continue;
00584                                 else {
00585                                         az += sym->get_az_alignment_offset(); // Align to the symmetry axes
00586                                 }
00587                         }
00588                         add_orientation(ret,az,alt);
00589                         if ( sym->is_h_sym() && inc_mirror && alt != altmin ) {
00590                                 add_orientation(ret,az,2.0f*altmin-alt);
00591                         }
00592                 }
00593         }
00594 
00595         return ret;
00596 }
00597 
00598 int SaffOrientationGenerator::get_orientations_tally(const Symmetry3D* const sym, const float& delta) const
00599 {
00600         bool inc_mirror = params.set_default("inc_mirror",false);
00601         Dict delimiters = sym->get_delimiters(inc_mirror);
00602         float altmax = delimiters["alt_max"];
00603         float azmax = delimiters["az_max"];
00604 
00605         float altmin = 0.0;
00606         // #If it's a h symmetry then the alt iterator starts at very close
00607         // #to the altmax... the object is a h symmetry then it knows its alt_min...
00608         if (sym->is_h_sym()){
00609                 altmin = delimiters["alt_min"];
00610                 if (inc_mirror) {
00611                         altmin -= (float) sym->get_params()["equator_range"];
00612                 }
00613         }
00614 
00615         float Deltaz = (float)(cos(altmax*EMConsts::deg2rad)-cos(altmin*EMConsts::deg2rad));
00616         float s = delta*M_PI/180.0f;
00617         float NFactor = 3.6f/s;
00618         float wedgeFactor = fabs( Deltaz*(azmax)/720.0f) ;
00619         int NumPoints   =  static_cast<int> (NFactor*NFactor*wedgeFactor);
00620 
00621         int tally = 0;
00622         if (!sym->is_h_sym()) ++tally;
00623         float az = 0.0f;
00624         float dz = (float)cos(altmin*EMConsts::deg2rad);
00625         for(int i = 1; i < NumPoints; ++i ){
00626                 float z = dz + Deltaz* (float)i/ float(NumPoints-1);
00627                 float r= sqrt(1.0f-z*z);
00628                 az = fmod(az + delta/r,azmax);
00629                 float alt = (float)(acos(z)*EMConsts::rad2deg);
00630                 if (sym->is_platonic_sym()) {
00631                         if ( sym->is_in_asym_unit(alt,az,inc_mirror) == false ) continue;
00632                 }
00633                 tally++;
00634         }
00635 
00636         return tally;
00637 }
00638 
00639 vector<Transform> SaffOrientationGenerator::gen_orientations(const Symmetry3D* const sym) const
00640 {
00641         float delta = params.set_default("delta", 0.0f);
00642         int n = params.set_default("n", 0);
00643 
00644         if ( delta <= 0 && n <= 0 ) throw InvalidParameterException("Error, you must specify a positive non-zero delta or n");
00645         if ( delta > 0 && n > 0 ) throw InvalidParameterException("Error, the delta and the n arguments are mutually exclusive");
00646 
00647         if ( n > 0 ) {
00648                 delta = get_optimal_delta(sym,n);
00649         }
00650 
00651 //      if ( sym->is_platonic_sym() ) return gen_platonic_orientations(sym, delta);
00652 
00653         bool inc_mirror = params.set_default("inc_mirror",false);
00654         Dict delimiters = sym->get_delimiters(inc_mirror);
00655         float altmax = delimiters["alt_max"];
00656         float azmax = delimiters["az_max"];
00657 
00658         float altmin = 0.0;
00659         // #If it's a h symmetry then the alt iterator starts at very close
00660         // #to the altmax... the object is a h symmetry then it knows its alt_min...
00661         if (sym->is_h_sym()){
00662                 altmin = delimiters["alt_min"];
00663                 if (inc_mirror) {
00664                         altmin -= (float) sym->get_params()["equator_range"];
00665                 }
00666         }
00667 
00668         float Deltaz = (float)(cos(altmax*EMConsts::deg2rad)-cos(altmin*EMConsts::deg2rad));
00669         float s = delta*M_PI/180.0f;
00670         float NFactor = 3.6f/s;
00671         float wedgeFactor = fabs( Deltaz*(azmax)/720.0f) ;
00672         int NumPoints   =  static_cast<int> (NFactor*NFactor*wedgeFactor);
00673 
00674         vector<Transform> ret;
00675 
00676         if (!sym->is_h_sym()) add_orientation(ret,0,0);
00677         float az = 0.0;
00678         float dz = (float)cos(altmin*EMConsts::deg2rad);
00679         for(int i = 1; i < NumPoints; ++i ){
00680                 float z = dz + Deltaz* (float)i/ float(NumPoints-1);
00681                 float r= sqrt(1.0f-z*z);
00682                 az = fmod(az + delta/r,azmax);
00683                 float alt = (float)(acos(z)*EMConsts::rad2deg);
00684                 if (sym->is_platonic_sym()) {
00685                         if ( sym->is_in_asym_unit(alt,az,inc_mirror) == false ) continue;
00686                         else {
00687                                 az += sym->get_az_alignment_offset(); // Align to the symmetry axes
00688                         }
00689                 }
00690                 add_orientation(ret,az,alt);
00691         }
00692 
00693         return ret;
00694 }
00695 
00696 int OptimumOrientationGenerator::get_orientations_tally(const Symmetry3D* const sym, const float& delta) const
00697 {
00698         string deltaoptname = params.set_default("use","saff");
00699         Dict a;
00700         a["inc_mirror"] = (bool)params.set_default("inc_mirror",false);
00701         OrientationGenerator *g = Factory < OrientationGenerator >::get(deltaoptname,a);
00702         if (g) {
00703                 int tally = g->get_orientations_tally(sym,delta);
00704                 delete g;
00705                 g = 0;
00706                 return tally;
00707         }
00708         else throw;
00709 }
00710 
00711 vector<Transform> OptimumOrientationGenerator::gen_orientations(const Symmetry3D* const sym) const
00712 {
00713         float delta = params.set_default("delta", 0.0f);
00714         int n = params.set_default("n", 0);
00715 
00716         bool inc_mirror = params.set_default("inc_mirror",false);
00717 
00718         if ( delta <= 0 && n <= 0 ) throw InvalidParameterException("Error, you must specify a positive non-zero delta or n");
00719         if ( delta > 0 && n > 0 ) throw InvalidParameterException("Error, the delta and the n arguments are mutually exclusive");
00720 
00721         string generatorname = params.set_default("use","saff");
00722 
00723         if ( n > 0 && generatorname != RandomOrientationGenerator::NAME ) {
00724                 params["delta"] = get_optimal_delta(sym,n);
00725                 params["n"] = (int)0;
00726         }
00727 
00728         // Force the orientation generator to include the mirror - this is because
00729         // We will enventually use it to generate orientations over the intire sphere
00730         // which is C1 symmetry, with the inc_mirror flag set to true
00731         params["inc_mirror"] = true;
00732         OrientationGenerator* g = Factory < OrientationGenerator >::get(generatorname);
00733         g->set_params(copy_relevant_params(g));
00734 
00735 
00736         // get the starting orientation distribution
00737         CSym* unit_sphere = new CSym();
00738         Dict nsym; nsym["nsym"] = 1; unit_sphere->set_params(nsym);
00739 
00740         vector<Transform> unitsphereorientations = g->gen_orientations(unit_sphere);
00741         delete g; g = 0;
00742         delete unit_sphere; unit_sphere = 0;
00743 
00744         vector<Vec3f> angles = optimize_distances(unitsphereorientations);
00745 
00746         vector<Transform> ret;
00747         for (vector<Vec3f>::const_iterator it = angles.begin(); it != angles.end(); ++it ) {
00748                 if ( sym->is_in_asym_unit((*it)[1],(*it)[0],inc_mirror) ) {
00749                         add_orientation(ret,(*it)[0],(*it)[1]);
00750                 }
00751         }
00752 
00753         // reset the params to what they were before they were acted upon by this class
00754         params["inc_mirror"] = inc_mirror;
00755         params["delta"] = delta;
00756         params["n"] = n;
00757 
00758         return ret;
00759 }
00760 
00761 vector<Vec3f> OptimumOrientationGenerator::optimize_distances(const vector<Transform>& v) const
00762 {
00763         vector<Vec3f> points;
00764 
00765         for (vector<Transform>::const_iterator it = v.begin(); it != v.end(); ++it ) {
00766                 points.push_back(Vec3f(0,0,1)*(*it));
00767         }
00768 
00769         if ( points.size() >= 2 ) {
00770                 int max_it = 100;
00771                 float percentage = 0.01f;
00772 
00773                 for ( int i = 0; i < max_it; ++i ){
00774                         unsigned int p1 = 0;
00775                         unsigned int p2 = 1;
00776 
00777                         float distsquared = (points[p1]-points[p2]).squared_length();
00778 
00779                         // Find the nearest points
00780                         for(unsigned int j = 0; j < points.size(); ++j) {
00781                                 for(unsigned int k = j+1; k < points.size(); ++k) {
00782                                         float d = (points[j]-points[k]).squared_length();
00783                                         if ( d < distsquared ) {
00784                                                 distsquared = d;
00785                                                 p1 = j;
00786                                                 p2 = k;
00787                                         }
00788                                 }
00789                         }
00790 
00791                         // Move them apart by a small fraction
00792                         Vec3f delta = percentage*(points[p2]-points[p1]);
00793 
00794                         points[p2] += delta;
00795                         points[p2].normalize();
00796                         points[p1] -= delta;
00797                         points[p1].normalize();
00798                 }
00799         }
00800 
00801         vector<Vec3f> ret;
00802         for (vector<Vec3f>::const_iterator it = points.begin(); it != points.end(); ++it ) {
00803                 float altitude = (float)(EMConsts::rad2deg*acos((*it)[2]));
00804                 float azimuth = (float)(EMConsts::rad2deg*atan2((*it)[1],(*it)[0]));
00805                 ret.push_back(Vec3f(90.0f+azimuth,altitude,0));
00806         }
00807 
00808         return ret;
00809 }
00810 // THIS IS DWOOLFORDS FIRST SHOT AT EXTRACTING PHIL'S PLATONIC STUFF FROM SPARX
00811 // It didn't work.
00812 // vector<Transform3D> SaffOrientationGenerator::gen_platonic_orientations(const Symmetry3D* const sym, const float& delta) const
00813 // {
00814 //      float scrunch = 0.9; //closeness factor to eliminate oversampling corners
00815 //
00816 //      float fudge; //# fudge is a factor used to adjust phi steps
00817 //      float m = static_cast<float>(sym->get_max_csym());
00818 //      if ( sym->get_name() == TetrahedralSym::NAME ) fudge=0.9;
00819 //      else if ( sym->get_name() == OctahedralSym::NAME ) fudge=0.8;
00820 //      else if ( sym->get_name() == IcosahedralSym::NAME) fudge=0.95;
00821 //      else throw; // this should not happen
00822 //
00823 //      float n=3.0;
00824 //      float OmegaR = 2.0*M_PI/m;
00825 //      float cosOmega= cos(OmegaR);
00826 //      int Edges  = static_cast<int>(2.0*m*n/(2.0*(m+n)-m*n));
00827 //      int Faces  = static_cast<int>(2*Edges/n);
00828 //      float Area   = 4*M_PI/Faces/3.0; // also equals  2*pi/3 + Omega
00829 //      float costhetac = cosOmega/(1-cosOmega);
00830 //      float deltaRad= delta*M_PI/180.0;
00831 //
00832 //      int     NumPoints = static_cast<int>(Area/(deltaRad*deltaRad));
00833 //      float fheight = 1.0f/sqrt(3.0f)/(tan(OmegaR/2.0f));
00834 //      float z0      = costhetac; // initialize loop
00835 //      float z       = z0;
00836 //      float phi     = 0;
00837 //      float Deltaz  = (1-costhetac);
00838 //
00839 //      vector<Transform3D> ret;
00840 //      ret.push_back(Transform3D(phi, acos(z)*EMConsts::rad2deg, 0 ));
00841 //
00842 //      vector<Vec3f> points;
00843 //      points.push_back(Vec3f(sin(acos(z))*cos(phi*EMConsts::deg2rad) ,  sin(acos(z))*sin(phi*EMConsts::deg2rad) , z) );
00844 //      //nLast=  [ sin(acos(z))*cos(phi*piOver) ,  sin(acos(z))*sin(phi*piOver) , z]
00845 //      //nVec.append(nLast)
00846 //
00847 //      for(int k = 0; k < NumPoints-1; ++k ) {
00848 //              z = z0 + Deltaz*(float)k/(float)(NumPoints-1);
00849 //              float r= sqrt(1-z*z);
00850 //              float phiMax =180.0*OmegaR/M_PI/2.0;
00851 //              // Is it higher than fhat or lower
00852 //              if (z<= fheight && false) {
00853 //                      float thetaR   = acos(z);
00854 //                      float cosStuff = (cos(thetaR)/sin(thetaR))*sqrt(1. - 2 *cosOmega);
00855 //                      phiMax   =  180.0*( OmegaR - acos(cosStuff))/M_PI;
00856 //              }
00857 //              float angleJump = fudge* delta/r;
00858 //              phi = fmod(phi + angleJump,phiMax);
00859 // //           anglesNew = [phi,180.0*acos(z)/pi,0.];
00860 // //           Vec3f nNew( sin(acos(z))*cos(phi*EMConsts::deg2rad) ,  sin(acos(z))*sin(phi*EMConsts::deg2rad) , z);
00861 // //           float mindiff = acos(nNew.dot(points[0]));
00862 // //           for(unsigned int l = 0; l < points.size(); ++l ) {
00863 // //                   float dx = acos(nNew.dot(points[l]));
00864 // //                   if (dx < mindiff ) mindiff = dx;
00865 // //           }
00866 // //           if (mindiff > (angleJump*EMConsts::deg2rad *scrunch) ){
00867 // //                   points.push_back(nNew);
00868 //                      cout << "phi " << phi << " alt " << acos(z)*EMConsts::rad2deg << " " << z << endl;
00869 //                      ret.push_back(Transform3D(phi, acos(z)*EMConsts::rad2deg, 0 ));
00870 // //           }
00871 //      }
00872 //      ret.push_back(Transform3D(0,0,0 ));
00873 //
00874 //      return ret;
00875 // }
00876 
00877 
00878 
00879 Symmetry3D::Symmetry3D() : cached_au_planes(0),cache_size(0),num_triangles(0),au_sym_triangles() {}
00880 Symmetry3D::~Symmetry3D() {
00881         if (cached_au_planes != 0 ) {
00882                 delete_au_planes();
00883         }
00884 }
00885 
00886 void verify(const Vec3f& tmp, float * plane, const string& message )
00887 {
00888         cout << message << " residual " << plane[0]*tmp[0]+plane[1]*tmp[1]+plane[2]*tmp[2] + plane[3]  << endl;
00889 }
00890 
00891 Transform Symmetry3D::reduce(const Transform& t, int n) const
00892 {
00893         // Determine which asym unit the given asym unit is in
00894         int soln = in_which_asym_unit(t);
00895 
00896         // This should never happen
00897         if ( soln == -1 ) {
00898                 cout << "error, no solution found!" << endl;
00899                 throw;
00900         }
00901 
00902         // Get the symmetry operation corresponding to the intersection asymmetric unit
00903         Transform nt = get_sym(soln);
00904         // Transpose it (invert it)
00905         nt.invert();
00906         // Now we can transform the argument orientation into the default asymmetric unit
00907         nt  = t*nt;
00908         // Now that we're at the default asymmetric unit, we can map into the requested asymmunit by doing this
00909         if ( n != 0 ) {
00910                 nt = nt*get_sym(n);
00911         }
00912         // Done!
00913         return nt;
00914 
00915 }
00916 
00917 int Symmetry3D::in_which_asym_unit(const Transform& t3d) const
00918 {
00919         // Here it is assumed that final destination of the orientation (as encapsulated in the t3d object) is
00920         // in the z direction, so in essence we will start in the direction z and 'undo' the orientation to get the real
00921         // direction
00922         Vec3f p(0,0,1);
00923 
00924         Transform o(t3d);
00925         // Orientations are alway transposed when dealing with asymmetric units, projections,etc
00926         // We take the transpose to 'undo' the transform and get the true direction of the point.
00927         o.invert();
00928         // Figure out where the point would end up. No we could just as easily not transpose and do
00929         // left multiplation (as in what occurs in the FourierReconstructor during slice insertion)
00930         p = o*p;
00931 
00932         return point_in_which_asym_unit(p);
00933 }
00934 
00935 
00936 void Symmetry3D::cache_au_planes() const {
00937         if (cached_au_planes == 0 ) {
00938                 vector< vector<Vec3f> > au_triangles = get_asym_unit_triangles(true);
00939                 num_triangles = au_triangles.size();
00940                 cache_size = get_nsym()*au_triangles.size();
00941 
00942                 cached_au_planes = new float*[cache_size];
00943                 float** fit = cached_au_planes;
00944                 for(int i =0; i < cache_size; ++i,++fit) {
00945                         float *t = new float[4];
00946                         *fit = t;
00947                 }
00948 
00949 
00950                 int k = 0;
00951                 for(int i = 0; i < get_nsym(); ++i) {
00952 
00953                         for( ncit it = au_triangles.begin(); it != au_triangles.end(); ++it, ++k)
00954                         {
00955                                 // For each given triangle
00956                                 vector<Vec3f> points = *it;
00957                                 if ( i != 0 ) {
00958                                         for (vector<Vec3f>::iterator iit = points.begin(); iit != points.end(); ++iit ) {
00959                                                 // Rotate the points in the triangle so that the triangle occupies the
00960                                                 // space of the current asymmetric unit
00961                                                 *iit = (*iit)*get_sym(i);
00962                                         }
00963                                 }
00964 
00965                                 au_sym_triangles.push_back(points);
00966 
00967                                 // Determine the equation of the plane for the points, store it in plane
00968                                 Util::equation_of_plane(points[0],points[2],points[1],cached_au_planes[k]);
00969                         }
00970                 }
00971         }
00972         else throw UnexpectedBehaviorException("Attempt to generate a cache when cache exists");
00973 }
00974 
00975 void Symmetry3D::delete_au_planes() {
00976         if (cached_au_planes == 0 ) throw UnexpectedBehaviorException("Attempt to delete a cache that does not exist");
00977         float** fit = cached_au_planes;
00978         for(int i =0; i < cache_size; ++i,++fit) {
00979                 if (*fit == 0) throw UnexpectedBehaviorException("Attempt to delete a cache that does not exist");
00980                 delete [] *fit;
00981                 *fit = 0;
00982         }
00983 
00984         delete [] cached_au_planes;
00985         cached_au_planes = 0;
00986 }
00987 
00988 int Symmetry3D::point_in_which_asym_unit(const Vec3f& p) const
00989 {
00990         if (cached_au_planes == 0) {
00991                 cache_au_planes();
00992         }
00993 
00994         int k = 0;
00995         for(int i = 0; i < get_nsym(); ++i) {
00996                 for( int j = 0; j < num_triangles; ++j,++k) {
00997                         vector<Vec3f> points = au_sym_triangles[k];
00998 
00999                         float* plane = cached_au_planes[k];
01000                         Vec3f tmp = p;
01001 
01002                         // Determine the intersection of p with the plane - do this by finding out how much p should be scaled by
01003                         float scale = plane[0]*tmp[0]+plane[1]*tmp[1]+plane[2]*tmp[2];
01004                         if ( scale != 0 )
01005                                 scale = -plane[3]/scale;
01006                         else {
01007                                 // parralel!
01008                                 continue;
01009                         }
01010 
01011                         // If the scale factor is less than zero, then p is definitely not in this asymmetric unit
01012                         if (scale <= 0) continue;
01013 
01014                         // This is the intersection point
01015                         Vec3f pp = tmp*scale;
01016 
01017                         // Now we have to see if the point p is inside the region bounded by the points, or if it is outside
01018                         // If it is inside the region then p is in this asymmetric unit.
01019 
01020                         // This formula take from FIXME fill in once I get to work
01021                         Vec3f v = points[2]-points[0];
01022                         Vec3f u = points[1]-points[0];
01023                         Vec3f w = pp - points[0];
01024 
01025                         float udotu = u.dot(u);
01026                         float udotv = u.dot(v);
01027                         float udotw = u.dot(w);
01028                         float vdotv = v.dot(v);
01029                         float vdotw = v.dot(w);
01030 
01031                         float d = 1.0f/(udotv*udotv - udotu*vdotv);
01032                         float s = udotv*vdotw - vdotv*udotw;
01033                         s *= d;
01034 
01035                         float t = udotv*udotw - udotu*vdotw;
01036                         t *= d;
01037 
01038                         // We've done a few multiplications, so detect when there are tiny residuals that may throw off the final
01039                         // decision
01040                         if (fabs(s) < Transform::ERR_LIMIT ) s = 0;
01041                         if (fabs(t) < Transform::ERR_LIMIT ) t = 0;
01042 
01043                         if ( fabs((fabs(s)-1.0)) < Transform::ERR_LIMIT ) s = 1;
01044                         if ( fabs((fabs(t)-1.0)) < Transform::ERR_LIMIT ) t = 1;
01045 
01046                         // The final decision, if this is true then we've hit the jackpot
01047                         if ( s >= 0 && t >= 0 && (s+t) <= 1 ) {
01048                                 return i;
01049                         }
01050                 }
01051         }
01052 
01053         return -1;
01054 }
01055 vector<Transform> Symmetry3D::get_touching_au_transforms(bool inc_mirror) const
01056 {
01057         vector<Transform>  ret;
01058         vector<int> hit_cache;
01059 
01060         vector<Vec3f> points = get_asym_unit_points(inc_mirror);
01061         // Warning, this is a gross hack because it is assuming that the asym_unit_points
01062         // returned by DSym are in a particular orientation with respect to symmetric axes
01063         // if the internals of DSym change it could change what we should do here...
01064         // but for the time being it will do
01065         if (inc_mirror && is_d_sym() && (get_nsym()/2 % 2 == 0)) {
01066                 Dict delim = get_delimiters(false);
01067                 float angle = (float)(EMConsts::deg2rad*float(delim["az_max"]));
01068                 float y = -cos(angle);
01069                 float x = sin(angle);
01070                 points.push_back(Vec3f(x,y,0));
01071         }
01072         else if ( is_d_sym() && (get_nsym()/2 % 2 == 1)) {
01073                 Dict delim = get_delimiters(false);
01074                 float angle = float(delim["az_max"])/2.0f;
01075 //              cout << "Odd dsym using " << angle << endl;
01076                 angle *= (float)EMConsts::deg2rad;
01077                 float y = -cos(angle);
01078                 float x = sin(angle);
01079                 points.push_back(Vec3f(x,y,0));
01080 
01081                 if ( inc_mirror ) {
01082                         angle = 3.0f*(float(delim["az_max"]))/2.0f;
01083                         angle *= (float)EMConsts::deg2rad;
01084                         float y = -cos(angle);
01085                         float x = sin(angle);
01086                         points.push_back(Vec3f(x,y,0));
01087                 }
01088         }
01089 
01090         typedef vector<Vec3f>::const_iterator const_point_it;
01091         for(const_point_it point = points.begin(); point != points.end(); ++point ) {
01092 
01093                 for(int i = 1; i < get_nsym(); ++i) {
01094 
01095                         if ( find(hit_cache.begin(),hit_cache.end(),i) != hit_cache.end() ) continue;
01096                         Transform t = get_sym(i);
01097                         Vec3f result = (*point)*t;
01098 
01099                         if (is_platonic_sym()) {
01100                                 for(const_point_it tmp = points.begin(); tmp != points.end(); ++tmp ) {
01101                                         Vec3f tt = result-(*tmp);
01102                                         if (tt.squared_length() < 0.01) {
01103                                                 hit_cache.push_back(i);
01104                                                 ret.push_back(t);
01105                                         }
01106 
01107                                 }
01108                         }
01109                         else {
01110                                 result -= *point;
01111                                 if (result.squared_length() < 0.05) {
01112                                         hit_cache.push_back(i);
01113                                         ret.push_back(t);
01114                                 }
01115                         }
01116                 }
01117 
01118         }
01119 
01120         return ret;
01121 }
01122 
01123 
01124 vector<Transform> Symmetry3D::get_syms() const
01125 {
01126 
01127         vector<Transform> ret;
01128 //      if (t.is_identity()) {
01129                 for(int i = 0; i < get_nsym(); ++i ) {
01130                         ret.push_back(get_sym(i));
01131                 }
01132 //      } else {
01133 //              for(int i = 0; i < get_nsym(); ++i ) {
01134 //                      ret.push_back(get_sym(i)*t);
01135 //              }
01136 //      }
01137         return ret;
01138 }
01139 
01140 vector<Transform> Symmetry3D::get_symmetries(const string& symmetry)
01141 {
01142         Symmetry3D* sym = Factory<Symmetry3D>::get(Util::str_to_lower(symmetry));
01143         vector<Transform> ret = sym->get_syms();
01144         delete sym;
01145         return ret;
01146 }
01147 
01148 // C Symmetry stuff
01149 Dict CSym::get_delimiters(const bool inc_mirror) const {
01150         Dict returnDict;
01151         // Get the parameters of interest
01152         int nsym = params.set_default("nsym",0);
01153         if ( nsym <= 0 ) throw InvalidValueException(nsym,"Error, you must specify a positive non zero nsym");
01154 
01155         if ( inc_mirror ) returnDict["alt_max"] = 180.0;
01156         else  returnDict["alt_max"] = 90.0;
01157 
01158         returnDict["az_max"] = 360.0/(float)nsym;
01159 
01160         return returnDict;
01161 }
01162 
01163 bool CSym::is_in_asym_unit(const float& altitude, const float& azimuth, const bool inc_mirror = false) const
01164 {
01165         Dict d = get_delimiters(inc_mirror);
01166         float alt_max = d["alt_max"];
01167         float az_max = d["az_max"];
01168 
01169         int nsym = params.set_default("nsym",0);
01170         if ( nsym != 1 && azimuth < 0) return false;
01171         if ( altitude <= alt_max && azimuth <= az_max ) return true;
01172         return false;
01173 }
01174 
01175 vector<vector<Vec3f> > CSym::get_asym_unit_triangles(bool inc_mirror ) const{
01176         vector<Vec3f> v = get_asym_unit_points(inc_mirror);
01177         int nsym = params.set_default("nsym",0);
01178 
01179         vector<vector<Vec3f> > ret;
01180         if (v.size() == 0) return ret; // nsym == 1 and inc_mirror == true, this is the entire sphere!
01181         if (nsym == 1 && !inc_mirror) {
01182                 Vec3f z(0,0,1);
01183                 vector<Vec3f> tmp;
01184                 tmp.push_back(z);
01185                 tmp.push_back(v[1]);
01186                 tmp.push_back(v[0]);
01187                 ret.push_back(tmp);
01188 
01189                 vector<Vec3f> tmp2;
01190                 tmp2.push_back(z);
01191                 tmp2.push_back(v[2]);
01192                 tmp2.push_back(v[1]);
01193                 ret.push_back(tmp2);
01194 
01195                 vector<Vec3f> tmp3;
01196                 tmp3.push_back(z);
01197                 tmp3.push_back(v[3]);
01198                 tmp3.push_back(v[2]);
01199                 ret.push_back(tmp3);
01200 
01201                 vector<Vec3f> tmp4;
01202                 tmp4.push_back(z);
01203                 tmp4.push_back(v[0]);
01204                 tmp4.push_back(v[3]);
01205                 ret.push_back(tmp4);
01206         }
01207         else if (nsym == 2 && inc_mirror) {
01208                 Vec3f x(1,0,0);
01209                 vector<Vec3f> tmp;
01210                 tmp.push_back(v[1]);
01211                 tmp.push_back(v[0]);
01212                 tmp.push_back(x);
01213                 ret.push_back(tmp);
01214 
01215                 vector<Vec3f> tmp2;
01216                 tmp2.push_back(v[2]);
01217                 tmp2.push_back(v[1]);
01218                 tmp2.push_back(x);
01219                 ret.push_back(tmp2);
01220 
01221                 vector<Vec3f> tmp3;
01222                 tmp3.push_back(v[3]);
01223                 tmp3.push_back(v[2]);
01224                 tmp3.push_back(x);
01225                 ret.push_back(tmp3);
01226 
01227                 vector<Vec3f> tmp4;
01228                 tmp4.push_back(v[0]);
01229                 tmp4.push_back(v[3]);
01230                 tmp4.push_back(x);
01231                 ret.push_back(tmp4);
01232         }
01233         else if (nsym == 2 && !inc_mirror) {
01234                 vector<Vec3f> tmp;
01235                 tmp.push_back(v[0]);
01236                 tmp.push_back(v[2]);
01237                 tmp.push_back(v[1]);
01238                 ret.push_back(tmp);
01239 
01240                 vector<Vec3f> tmp2;
01241                 tmp2.push_back(v[2]);
01242                 tmp2.push_back(v[0]);
01243                 tmp2.push_back(v[3]);
01244                 ret.push_back(tmp2);
01245         }
01246         else if (v.size() == 3) {
01247                 vector<Vec3f> tmp;
01248                 tmp.push_back(v[0]);
01249                 tmp.push_back(v[2]);
01250                 tmp.push_back(v[1]);
01251                 ret.push_back(tmp);
01252         }
01253         else if (v.size() == 4) {
01254                 vector<Vec3f> tmp;
01255                 tmp.push_back(v[0]);
01256                 tmp.push_back(v[3]);
01257                 tmp.push_back(v[1]);
01258                 ret.push_back(tmp);
01259 
01260                 vector<Vec3f> tmp2;
01261                 tmp2.push_back(v[1]);
01262                 tmp2.push_back(v[3]);
01263                 tmp2.push_back(v[2]);
01264                 ret.push_back(tmp2);
01265         }
01266 
01267         return ret;
01268 }
01269 
01270 vector<Vec3f> CSym::get_asym_unit_points(bool inc_mirror) const
01271 {
01272         Dict delim = get_delimiters(inc_mirror);
01273         int nsym = params.set_default("nsym",0);
01274         vector<Vec3f> ret;
01275 
01276         if ( nsym == 1 ) {
01277                 if (inc_mirror == false ) {
01278                         ret.push_back(Vec3f(0,-1,0));
01279                         ret.push_back(Vec3f(1,0,0));
01280                         ret.push_back(Vec3f(0,1,0));
01281                         ret.push_back(Vec3f(-1,0,0));
01282                 }
01283                 // else return ret; // an empty vector! this is fine
01284         }
01285         else if (nsym == 2 && !inc_mirror) {
01286                 ret.push_back(Vec3f(0,0,1));
01287                 ret.push_back(Vec3f(0,-1,0));
01288                 ret.push_back(Vec3f(1,0,0));
01289                 ret.push_back(Vec3f(0,1,0));
01290         }
01291         else {
01292                 ret.push_back(Vec3f(0,0,1));
01293                 ret.push_back(Vec3f(0,-1,0));
01294                 if (inc_mirror == true) {
01295                         ret.push_back(Vec3f(0,0,-1));
01296                 }
01297                 float angle = (float)(EMConsts::deg2rad*float(delim["az_max"]));
01298                 float y = -cos(angle);
01299                 float x = sin(angle);
01300                 ret.push_back(Vec3f(x,y,0));
01301         }
01302 
01303         return ret;
01304 
01305 }
01306 
01307 Transform CSym::get_sym(const int n) const {
01308         int nsym = params.set_default("nsym",0);
01309         if ( nsym <= 0 ) throw InvalidValueException(n,"Error, you must specify a positive non zero nsym");
01310 
01311         Dict d("type","eman");
01312         // courtesy of Phil Baldwin
01313         d["az"] = (n%nsym) * 360.0f / nsym;
01314         d["alt"] = 0.0f;
01315         d["phi"] = 0.0f;
01316         return Transform(d);
01317 }
01318 
01319 // D symmetry stuff
01320 Dict DSym::get_delimiters(const bool inc_mirror) const {
01321         Dict returnDict;
01322 
01323         // Get the parameters of interest
01324         int nsym = params.set_default("nsym",0);
01325         if ( nsym <= 0 ) throw InvalidValueException(nsym,"Error, you must specify a positive non zero nsym");
01326 
01327         returnDict["alt_max"] = 90.0;
01328 
01329         if ( inc_mirror )  returnDict["az_max"] = 360.0/(float)nsym;
01330         else returnDict["az_max"] = 180.0/(float)nsym;
01331 
01332         return returnDict;
01333 }
01334 
01335 bool DSym::is_in_asym_unit(const float& altitude, const float& azimuth, const bool inc_mirror = false) const
01336 {
01337         Dict d = get_delimiters(inc_mirror);
01338         float alt_max = d["alt_max"];
01339         float az_max = d["az_max"];
01340 
01341         int nsym = params.set_default("nsym",0);
01342 
01343         if ( nsym == 1 && inc_mirror ) {
01344                 if (altitude >= 0 && altitude <= alt_max && azimuth <= az_max ) return true;
01345         }
01346         else {
01347                 if ( altitude >= 0 && altitude <= alt_max && azimuth <= az_max && azimuth >= 0 ) return true;
01348         }
01349         return false;
01350 }
01351 
01352 Transform DSym::get_sym(const int n) const
01353 {
01354         int nsym = 2*params.set_default("nsym",0);
01355         if ( nsym <= 0 ) throw InvalidValueException(n,"Error, you must specify a positive non zero nsym");
01356 
01357         Dict d("type","eman");
01358         // courtesy of Phil Baldwin
01359         if (n >= nsym / 2) {
01360                 d["az"] = ( (n%nsym) - nsym/2) * 360.0f / (nsym / 2);
01361                 d["alt"] = 180.0f;
01362                 d["phi"] = 0.0f;
01363         }
01364         else {
01365                 d["az"] = (n%nsym) * 360.0f / (nsym / 2);
01366                 d["alt"] = 0.0f;
01367                 d["phi"] = 0.0f;
01368         }
01369         return Transform(d);
01370 }
01371 
01372 vector<vector<Vec3f> > DSym::get_asym_unit_triangles(bool inc_mirror) const{
01373         vector<Vec3f> v = get_asym_unit_points(inc_mirror);
01374         int nsym = params.set_default("nsym",0);
01375         vector<vector<Vec3f> > ret;
01376         if ( (nsym == 1 && inc_mirror == false) || (nsym == 2 && inc_mirror)) {
01377                 vector<Vec3f> tmp;
01378                 tmp.push_back(v[0]);
01379                 tmp.push_back(v[2]);
01380                 tmp.push_back(v[1]);
01381                 ret.push_back(tmp);
01382 
01383                 vector<Vec3f> tmp2;
01384                 tmp2.push_back(v[2]);
01385                 tmp2.push_back(v[0]);
01386                 tmp2.push_back(v[3]);
01387                 ret.push_back(tmp2);
01388         }
01389         else if (nsym == 1) {
01390                 Vec3f z(0,0,1);
01391                 vector<Vec3f> tmp;
01392                 tmp.push_back(z);
01393                 tmp.push_back(v[1]);
01394                 tmp.push_back(v[0]);
01395                 ret.push_back(tmp);
01396 
01397                 vector<Vec3f> tmp2;
01398                 tmp2.push_back(z);
01399                 tmp2.push_back(v[2]);
01400                 tmp2.push_back(v[1]);
01401                 ret.push_back(tmp2);
01402 
01403                 vector<Vec3f> tmp3;
01404                 tmp3.push_back(z);
01405                 tmp3.push_back(v[3]);
01406                 tmp3.push_back(v[2]);
01407                 ret.push_back(tmp3);
01408 
01409                 vector<Vec3f> tmp4;
01410                 tmp4.push_back(z);
01411                 tmp4.push_back(v[0]);
01412                 tmp4.push_back(v[3]);
01413                 ret.push_back(tmp4);
01414         }
01415         else {
01416 //              if v.size() == 3
01417                 vector<Vec3f> tmp;
01418                 tmp.push_back(v[0]);
01419                 tmp.push_back(v[2]);
01420                 tmp.push_back(v[1]);
01421                 ret.push_back(tmp);
01422         }
01423 
01424         return ret;
01425 }
01426 
01427 vector<Vec3f> DSym::get_asym_unit_points(bool inc_mirror) const
01428 {
01429         Dict delim = get_delimiters(inc_mirror);
01430 
01431         vector<Vec3f> ret;
01432         int nsym = params.set_default("nsym",0);
01433         if ( nsym == 1 ) {
01434                 if (inc_mirror == false ) {
01435                         ret.push_back(Vec3f(0,0,1));
01436                         ret.push_back(Vec3f(0,-1,0));
01437                         ret.push_back(Vec3f(1,0,0));
01438                         ret.push_back(Vec3f(0,1,0));
01439                 }
01440                 else {
01441                         ret.push_back(Vec3f(0,-1,0));
01442                         ret.push_back(Vec3f(1,0,0));
01443                         ret.push_back(Vec3f(0,1,0));
01444                         ret.push_back(Vec3f(-1,0,0));
01445                 }
01446         }
01447         else if ( nsym == 2 && inc_mirror ) {
01448                 ret.push_back(Vec3f(0,0,1));
01449                 ret.push_back(Vec3f(0,-1,0));
01450                 ret.push_back(Vec3f(1,0,0));
01451                 ret.push_back(Vec3f(0,1,0));
01452         }
01453         else {
01454                 float angle = (float)(EMConsts::deg2rad*float(delim["az_max"]));
01455                 ret.push_back(Vec3f(0,0,1));
01456                 ret.push_back(Vec3f(0,-1,0));
01457                 float y = -cos(angle);
01458                 float x = sin(angle);
01459                 ret.push_back(Vec3f(x,y,0));
01460         }
01461 
01462         return ret;
01463 
01464 }
01465 
01466 // H symmetry stuff
01467 Dict HSym::get_delimiters(const bool) const {
01468         Dict returnDict;
01469 
01470         // Get the parameters of interest
01471         int nsym = params.set_default("nsym",0);
01472         if ( nsym <= 0 ) throw InvalidValueException(nsym,"Error, you must specify a positive non zero nsym");
01473 
01474         float equator_range = params.set_default("equator_range",5.0f);
01475 
01476         returnDict["alt_max"] = 90.0 + equator_range;
01477 
01478         returnDict["alt_min"] = 90.0;
01479 
01480         returnDict["az_max"] = 360.0/(float)nsym;
01481 
01482         return returnDict;
01483 }
01484 
01485 bool HSym::is_in_asym_unit(const float& altitude, const float& azimuth, const bool inc_mirror = false) const
01486 {
01487         Dict d = get_delimiters(inc_mirror);
01488         float alt_max = d["alt_max"];
01489         float alt_min = d["alt_min"];
01490 
01491         if (inc_mirror) {
01492                 float e = params.set_default("equator_range",5.0f);
01493                 alt_min -= e;
01494         }
01495 
01496         float az_max = d["az_max"];
01497 
01498         if ( altitude >=alt_min && altitude <= alt_max && azimuth <= az_max && azimuth >= 0 ) return true;
01499         return false;
01500 }
01501 
01502 vector<vector<Vec3f> > HSym::get_asym_unit_triangles(bool inc_mirror) const{
01503 
01504         vector<vector<Vec3f> > ret;
01505         return ret;
01506 }
01507 
01508 vector<Vec3f> HSym::get_asym_unit_points(bool inc_mirror) const
01509 {
01510         vector<Vec3f> ret;
01511 
01512         Dict delim = get_delimiters(inc_mirror);
01513         int nsym = params.set_default("nsym",0);
01514         float az = -(float)delim["az_max"];
01515 
01516 
01517         bool tracing_arcs = false;
01518 
01519 
01520         if ( !tracing_arcs) {
01521                 Vec3f a(0,-1,0);
01522                 ret.push_back(a);
01523 
01524                 if ( nsym > 2 ) {
01525                         Dict d("type","eman");
01526                         d["phi"] = 0.0f;
01527                         d["alt"] = 0.0f;
01528                         d["az"] = az;
01529                         Vec3f b = Transform(d)*a;
01530                         ret.push_back(b);
01531                 }
01532                 else
01533                 {
01534                         ret.push_back(Vec3f(1,0,0));
01535 
01536                         ret.push_back(Vec3f(0,1,0));
01537 
01538                         if ( nsym == 1 ) {
01539                                 ret.push_back(Vec3f(-1,0,0));
01540                                 ret.push_back(a);
01541                         }
01542                 }
01543         }
01544         return ret;
01545 
01546 }
01547 
01548 Transform HSym::get_sym(const int n) const
01549 {
01550         float daz = params.set_default("daz",0.0f);
01551         float apix = params.set_default("apix",1.0f);
01552         float dz = params.set_default("tz", 0)/apix;
01553 
01554         Dict d("type","eman");
01555         // courtesy of Phil Baldwin
01556         d["az"] = n * daz;
01557         d["alt"] = 0.0f;
01558         d["phi"] = 0.0f;
01559         Transform ret(d);
01560         ret.set_trans(0,0,n*dz);
01561         return ret;
01562 }
01563 
01564 // Generic platonic symmetry stuff
01565 void PlatonicSym::init()
01566 {
01567         //See the manuscript "The Transform Class in Sparx and EMAN2", Baldwin & Penczek 2007. J. Struct. Biol. 157 (250-261)
01568         //In particular see pages 257-259
01569         //cap_sig is capital sigma in the Baldwin paper
01570         float cap_sig =  2.0f*M_PI/ get_max_csym();
01571         //In EMAN2 projection cap_sig is really az_max
01572         platonic_params["az_max"] = cap_sig;
01573 
01574         // Alpha is the angle between (immediately) neighborhing 3 fold axes of symmetry
01575         // This follows the conventions in the Baldwin paper
01576         float alpha = acos(1.0f/(sqrtf(3.0f)*tan(cap_sig/2.0f)));
01577         // In EMAN2 projection alpha is really al_maz
01578         platonic_params["alt_max"] = alpha;
01579 
01580         // This is half of "theta_c" as in the conventions of the Balwin paper. See also http://blake.bcm.edu/emanwiki/EMAN2/Symmetry.
01581         platonic_params["theta_c_on_two"] = 1.0/2.0*acos( cos(cap_sig)/(1.0-cos(cap_sig)));
01582 
01583 }
01584 
01585 
01586 Dict PlatonicSym::get_delimiters(const bool inc_mirror) const
01587 {
01588         Dict ret;
01589         ret["az_max"] = EMConsts::rad2deg * (float) platonic_params["az_max"];
01590         // For icos and oct symmetries, excluding the mirror means halving az_maz
01591         if ( inc_mirror == false )
01592                 if ( get_name() ==  IcosahedralSym::NAME || get_name() == OctahedralSym::NAME )
01593                         ret["az_max"] = 0.5*EMConsts::rad2deg * (float) platonic_params["az_max"];
01594                 //else
01595                 //the alt_max variable should probably be altered if the symmetry is tet, but
01596                 //this is taken care of in TetSym::is_in_asym_unit
01597 
01598         ret["alt_max"] = EMConsts::rad2deg * (float) platonic_params["alt_max"];
01599         return ret;
01600 }
01601 
01602 //.Warning, this function only returns valid answers for octahedral and icosahedral symmetries.
01603 bool PlatonicSym::is_in_asym_unit(const float& altitude, const float& azimuth, const bool inc_mirror) const
01604 {
01605         Dict d = get_delimiters(inc_mirror);
01606         float alt_max = d["alt_max"];
01607         float az_max = d["az_max"];
01608 
01609         if ( altitude >= 0 &&  altitude <= alt_max && azimuth <= az_max && azimuth >= 0) {
01610 
01611                 // Convert azimuth to radians
01612                 float tmpaz = (float)(EMConsts::deg2rad * azimuth);
01613 
01614                 float cap_sig = platonic_params["az_max"];
01615                 float alt_max = platonic_params["alt_max"];
01616                 if ( tmpaz > ( cap_sig/2.0 ) )tmpaz = cap_sig - tmpaz;
01617 
01618                 float lower_alt_bound = platonic_alt_lower_bound(tmpaz, alt_max );
01619 
01620                 // convert altitude to radians
01621                 float tmpalt = (float)(EMConsts::deg2rad * altitude);
01622                 if ( lower_alt_bound > tmpalt ) {
01623                         if ( inc_mirror == false )
01624                         {
01625                                 if ( cap_sig/2.0 < tmpaz ) return false;
01626                                 else return true;
01627                         }
01628                         else return true;
01629                 }
01630                 return false;
01631         }
01632         return false;
01633 }
01634 
01635 float PlatonicSym::platonic_alt_lower_bound(const float& azimuth, const float& alpha) const
01636 {
01637         float cap_sig = platonic_params["az_max"];
01638         float theta_c_on_two = platonic_params["theta_c_on_two"];
01639 
01640         float baldwin_lower_alt_bound = sin(cap_sig/2.0f-azimuth)/tan(theta_c_on_two);
01641         baldwin_lower_alt_bound += sin(azimuth)/tan(alpha);
01642         baldwin_lower_alt_bound *= 1/sin(cap_sig/2.0f);
01643         baldwin_lower_alt_bound = atan(1/baldwin_lower_alt_bound);
01644 
01645         return baldwin_lower_alt_bound;
01646 }
01647 
01648 vector<vector<Vec3f> > PlatonicSym::get_asym_unit_triangles(bool inc_mirror) const{
01649         vector<Vec3f> v = get_asym_unit_points(inc_mirror);
01650         vector<vector<Vec3f> > ret;
01651         if (v.size() == 3) {
01652                 vector<Vec3f> tmp;
01653                 tmp.push_back(v[0]);
01654                 tmp.push_back(v[2]);
01655                 tmp.push_back(v[1]);
01656                 ret.push_back(tmp);
01657         }
01658         else /* v.size() == 4*/ {
01659                 vector<Vec3f> tmp;
01660                 tmp.push_back(v[0]);
01661                 tmp.push_back(v[2]);
01662                 tmp.push_back(v[1]);
01663                 ret.push_back(tmp);
01664 
01665                 vector<Vec3f> tmp2;
01666                 tmp2.push_back(v[0]);
01667                 tmp2.push_back(v[3]);
01668                 tmp2.push_back(v[2]);
01669                 ret.push_back(tmp2);
01670         }
01671 
01672         return ret;
01673 }
01674 
01675 vector<Vec3f> PlatonicSym::get_asym_unit_points(bool inc_mirror) const
01676 {
01677         vector<Vec3f> ret;
01678 
01679         Vec3f b = Vec3f(0,0,1);
01680         ret.push_back(b);
01681         float theta_c_on_two = (float)platonic_params["theta_c_on_two"]; // already in radians
01682         float theta_c = 2*theta_c_on_two;
01683 
01684         Vec3f c_on_two = Vec3f(0,-sin(theta_c_on_two),cos(theta_c_on_two));
01685         Vec3f c = Vec3f(0,-sin(theta_c),cos(theta_c));
01686         ret.push_back(c_on_two);
01687 
01688         float cap_sig = platonic_params["az_max"];
01689         Vec3f a = Vec3f(sin(theta_c)*sin(cap_sig),-sin(theta_c)*cos(cap_sig),cos(theta_c));
01690 
01691         Vec3f f = a+b+c;
01692         f.normalize();
01693 
01694         ret.push_back(f);
01695 
01696         if ( inc_mirror ) {
01697                 Vec3f a_on_two = Vec3f(sin(theta_c_on_two)*sin(cap_sig),-sin(theta_c_on_two)*cos(cap_sig),cos(theta_c_on_two));
01698                 ret.push_back(a_on_two);
01699         }
01700 
01701         if ( get_az_alignment_offset() != 0 ) {
01702                 Dict d("type","eman");
01703                 d["az"] = get_az_alignment_offset();
01704                 d["phi"] = 0.0f;
01705                 d["alt"] = 0.0f;
01706                 Transform t(d);
01707                 for (vector<Vec3f>::iterator it = ret.begin(); it != ret.end(); ++it )
01708                 {
01709                         *it = (*it)*t;
01710                 }
01711         }
01712         //
01713         return ret;
01714 
01715 }
01716 
01717 float IcosahedralSym::get_az_alignment_offset() const { return 234.0; } // This offset positions a 3 fold axis on the positive x axis
01718 
01719 Transform IcosahedralSym::get_sym(const int n) const
01720 {
01721         // These rotations courtesy of Phil Baldwin
01722         static double  lvl0=0.; //  there is one pentagon on top; five-fold along z
01723         static double  lvl1= 63.4349; // that is atan(2)  // there are 5 pentagons with centers at this height (angle)
01724         static double  lvl2=116.5651; //that is 180-lvl1  // there are 5 pentagons with centers at this height (angle)
01725         static double lvl3=180.0;
01726 
01727         static double ICOS[180] = { // This is with a pentagon normal to z
01728                 0,lvl0,0,    0,lvl0,288,   0,lvl0,216,   0,lvl0,144,  0,lvl0,72,
01729   0,lvl1,36,   0,lvl1,324,   0,lvl1,252,   0,lvl1,180,  0,lvl1,108,
01730   72,lvl1,36,  72,lvl1,324,  72,lvl1,252,  72,lvl1,180,  72,lvl1,108,
01731   144,lvl1,36, 144,lvl1,324, 144,lvl1,252, 144,lvl1,180, 144,lvl1,108,
01732   216,lvl1,36, 216,lvl1,324, 216,lvl1,252, 216,lvl1,180, 216,lvl1,108,
01733   288,lvl1,36, 288,lvl1,324, 288,lvl1,252, 288,lvl1,180, 288,lvl1,108,
01734   36,lvl2,0,   36,lvl2,288,  36,lvl2,216,  36,lvl2,144,  36,lvl2,72,
01735   108,lvl2,0,  108,lvl2,288, 108,lvl2,216, 108,lvl2,144, 108,lvl2,72,
01736   180,lvl2,0,  180,lvl2,288, 180,lvl2,216, 180,lvl2,144, 180,lvl2,72,
01737   252,lvl2,0,  252,lvl2,288, 252,lvl2,216, 252,lvl2,144, 252,lvl2,72,
01738   324,lvl2,0,  324,lvl2,288, 324,lvl2,216, 324,lvl2,144, 324,lvl2,72,
01739   0,lvl3,0,    0,lvl3,288,   0,lvl3,216,   0,lvl3,144,   0,lvl3,72
01740         };
01741 
01742         int idx = n % 60;
01743         Dict d("type","eman");
01744 //      Transform3D ret;
01745         if (get_az_alignment_offset() == 234.0) {
01746                 d["az"] =(float)ICOS[idx * 3 ]+90;
01747                 d["alt"] = (float)ICOS[idx * 3 + 1];
01748                 d["phi"] = (float)ICOS[idx * 3 + 2]-90;
01749 //              ret.set_rotation((float)ICOS[idx * 3 ]+90,(float)ICOS[idx * 3 + 1], (float)ICOS[idx * 3 + 2]-90);
01750         }
01751         else {
01752                 d["az"] =(float)(float)ICOS[idx * 3 ];
01753                 d["alt"] = (float)ICOS[idx * 3 + 1];
01754                 d["phi"] = (float)ICOS[idx * 3 + 2];
01755 //              ret.set_rotation((float)ICOS[idx * 3 ],(float)ICOS[idx * 3 + 1], (float)ICOS[idx * 3 + 2]);
01756         }
01757 
01758 //      ret.set_rotation((float)ICOS[idx * 3 ],(float)ICOS[idx * 3 + 1], (float)ICOS[idx * 3 + 2]);
01759 //      if ( get_az_alignment_offset() != 0 ) {
01760 //              Transform3D t(get_az_alignment_offset(),0,0);
01761 //              ret = t*ret;
01762 //      }
01763         return Transform(d);
01764 
01765 }
01766 
01767 Transform OctahedralSym::get_sym(const int n) const
01768 {
01769         // These rotations courtesy of Phil Baldwin
01770         // We have placed the OCT symmetry group with a face along the z-axis
01771         static double lvl0=0.;
01772         static double lvl1=90.;
01773         static double lvl2=180.;
01774 
01775         static double OCT[72] = {// This is with a face of a cube along z
01776                 0,lvl0,0,   0,lvl0,90,    0,lvl0,180,    0,lvl0,270,
01777   0,lvl1,0,   0,lvl1,90,    0,lvl1,180,    0,lvl1,270,
01778   90,lvl1,0,  90,lvl1,90,   90,lvl1,180,   90,lvl1,270,
01779   180,lvl1,0, 180,lvl1,90,  180,lvl1,180,  180,lvl1,270,
01780   270,lvl1,0, 270,lvl1,90,  270,lvl1,180,  270,lvl1,270,
01781   0,lvl2,0,   0,lvl2,90,    0,lvl2,180,    0,lvl2,270
01782         };
01783 
01784         int idx = n % 24;
01785 //      Transform3D ret;
01786 //      ret.set_rotation((float)OCT[idx * 3 ],(float)OCT[idx * 3 + 1], (float)OCT[idx * 3 + 2] );
01787         Dict d("type","eman");
01788         d["az"] = (float)OCT[idx * 3 ];
01789         d["alt"] = (float)OCT[idx * 3 + 1];
01790         d["phi"] = (float)OCT[idx * 3 + 2];
01791         return Transform(d);
01792 
01793 }
01794 
01795 float TetrahedralSym::get_az_alignment_offset() const { return  0.0; }
01796 
01797 bool TetrahedralSym::is_in_asym_unit(const float& altitude, const float& azimuth, const bool inc_mirror) const
01798 {
01799         Dict d = get_delimiters(inc_mirror);
01800         float alt_max = d["alt_max"];
01801         float az_max = d["az_max"];
01802 
01803         if ( altitude >= 0 &&  altitude <= alt_max && azimuth <= az_max && azimuth >= 0) {
01804                 // convert azimuth to radians
01805                 float tmpaz = (float)(EMConsts::deg2rad * azimuth);
01806 
01807                 float cap_sig = platonic_params["az_max"];
01808                 float alt_max = platonic_params["alt_max"];
01809                 if ( tmpaz > ( cap_sig/2.0 ) )tmpaz = cap_sig - tmpaz;
01810 
01811                 float lower_alt_bound = platonic_alt_lower_bound(tmpaz, alt_max );
01812 
01813                 // convert altitude to radians
01814                 float tmpalt = (float)(EMConsts::deg2rad * altitude);
01815                 if ( lower_alt_bound > tmpalt ) {
01816                         if ( !inc_mirror ) {
01817                                 float upper_alt_bound = platonic_alt_lower_bound( tmpaz, alt_max/2.0f);
01818                                 // you could change the "<" to a ">" here to get the other mirror part of the asym unit
01819                                 if ( upper_alt_bound < tmpalt ) return false;
01820                                 else return true;
01821                         }
01822                         else return true;
01823                 }
01824                 return false;
01825         }
01826         else return false;
01827 }
01828 
01829 
01830 Transform TetrahedralSym::get_sym(const int n) const
01831 {
01832         // These rotations courtesy of Phil Baldwin
01833          // It has n=m=3; F=4, E=6=nF/2, V=4=nF/m
01834         static double lvl0=0;         // There is a face along z
01835         static double lvl1=109.4712;  //  that is acos(-1/3)  // There  are 3 faces at this angle
01836 
01837         static double TET[36] = {// This is with the face along z
01838                 0,lvl0,0,   0,lvl0,120,    0,lvl0,240,
01839   0,lvl1,60,   0,lvl1,180,    0,lvl1,300,
01840   120,lvl1,60, 120,lvl1,180,  120,lvl1,300,
01841   240,lvl1,60, 240,lvl1,180,  240,lvl1,300
01842         };
01843         //
01844         int idx = n % 12;
01845 //      Transform3D ret;
01846 //      ret.set_rotation((float)TET[idx * 3 ],(float)TET[idx * 3 + 1], (float)TET[idx * 3 + 2] );
01847 //      return ret;
01848 
01849         Dict d("type","eman");
01850         d["az"] = (float)TET[idx * 3 ];
01851         d["alt"] = (float)TET[idx * 3 + 1];
01852         d["phi"] = (float)TET[idx * 3 + 2];
01853         return Transform(d);
01854 
01855 }
01856 
01857 
01858 vector<Vec3f> TetrahedralSym::get_asym_unit_points(bool inc_mirror) const
01859 {
01860         vector<Vec3f> ret;
01861 
01862         Vec3f b = Vec3f(0,0,1);
01863         ret.push_back(b);
01864         float theta_c_on_two = (float)platonic_params["theta_c_on_two"]; // already in radians
01865         float theta_c = 2*theta_c_on_two;
01866 
01867         Vec3f c_on_two = Vec3f(0,-sin(theta_c_on_two),cos(theta_c_on_two));
01868         Vec3f c = Vec3f(0,-sin(theta_c),cos(theta_c));
01869         ret.push_back(c_on_two);
01870         float cap_sig = platonic_params["az_max"];
01871         if ( inc_mirror ) {
01872                 Vec3f a = Vec3f(sin(theta_c)*sin(cap_sig),-sin(theta_c)*cos(cap_sig),cos(theta_c));
01873 
01874                 Vec3f f = a+b+c;
01875                 f.normalize();
01876 
01877                 ret.push_back(f);
01878         }
01879 
01880         Vec3f a_on_two = Vec3f(sin(theta_c_on_two)*sin(cap_sig),-sin(theta_c_on_two)*cos(cap_sig),cos(theta_c_on_two));
01881         ret.push_back(a_on_two);
01882 
01883 
01884         if ( get_az_alignment_offset() != 0 ) {
01885                 Dict d("type","eman");
01886                 d["az"] = get_az_alignment_offset();
01887                 d["phi"] = 0.0f;
01888                 d["alt"] = 0.0f;
01889                 Transform t(d);
01890                 for (vector<Vec3f>::iterator it = ret.begin(); it != ret.end(); ++it )
01891                 {
01892                         *it = (*it)*t;
01893                 }
01894         }
01895 
01896         return ret;
01897 }
01898 
01899 
01900 

Generated on Fri Apr 30 15:38:57 2010 for EMAN2 by  doxygen 1.3.9.1