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