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

Generated on Thu May 3 10:06:28 2012 for EMAN2 by  doxygen 1.4.7