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

Generated on Tue Jun 11 12:40:26 2013 for EMAN2 by  doxygen 1.4.7