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

symmetry.cpp

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

Generated on Tue Jul 12 13:49:29 2011 for EMAN2 by  doxygen 1.3.9.1