@@ -235,6 +235,11 @@ void GlyphGeom::addEllipsoid(const Point& p, Tensor& t, double scale, int resolu
235235 generateEllipsoid (p, t, scale, resolution, color, false , normalize);
236236}
237237
238+ void GlyphGeom::addSuperEllipsoid (const Point& p, Tensor& t, double scale, int resolution, const ColorRGB& color, bool normalize, double emphasis)
239+ {
240+ generateSuperEllipsoid (p, t, scale, resolution, color, normalize, emphasis);
241+ }
242+
238243void GlyphGeom::addCylinder (const Point& p1, const Point& p2, double radius, int resolution,
239244 const ColorRGB& color1, const ColorRGB& color2,
240245 bool renderBase1, bool renderBase2)
@@ -980,6 +985,199 @@ void GlyphGeom::generateEllipsoid(const Point& center, Tensor& t, double scale,
980985 }
981986}
982987
988+ inline double spow (double e, double x)
989+ {
990+ // This for round off of very small numbers.
991+ if ( abs ( e ) < 1.0e-6 )
992+ e = 0.0 ;
993+
994+ if (e < 0.0 )
995+ {
996+ return (double )(pow (abs (e), x) * -1.0 );
997+ }
998+ else
999+ {
1000+ return (double )(pow (e, x));
1001+ }
1002+ }
1003+
1004+ void GlyphGeom::generateSuperEllipsoid (const Point& center, Tensor& t, double scale, int resolution, const ColorRGB& color, bool normalize, double emphasis)
1005+ {
1006+ double zeroThreshold = 0.000001 ;
1007+ std::vector<Vector> eigvectors (3 );
1008+ t.get_eigenvectors (eigvectors[0 ], eigvectors[1 ], eigvectors[2 ]);
1009+
1010+ double eigval1, eigval2, eigval3;
1011+ t.get_eigenvalues (eigval1, eigval2, eigval3);
1012+ Vector eigvals = Vector (fabs (eigval1), fabs (eigval2), fabs (eigval3));
1013+ if (normalize)
1014+ eigvals.normalize ();
1015+ eigvals *= scale;
1016+
1017+ // Checks to see if eigenvalues are close to 0
1018+ bool eig_x_0 = eigvals.x () <= zeroThreshold;
1019+ bool eig_y_0 = eigvals.y () <= zeroThreshold;
1020+ bool eig_z_0 = eigvals.z () <= zeroThreshold;
1021+
1022+ // Set to 0 if below threshold
1023+ eigvals[0 ] = (!eig_x_0) * eigvals[0 ];
1024+ eigvals[1 ] = (!eig_y_0) * eigvals[1 ];
1025+ eigvals[2 ] = (!eig_z_0) * eigvals[2 ];
1026+
1027+ bool flatTensor = (eig_x_0 + eig_y_0 + eig_z_0) >= 1 ;
1028+ Vector zero_norm;
1029+
1030+ if (flatTensor)
1031+ {
1032+ reorderTensor (eigvectors, eigvals);
1033+
1034+ eig_x_0 = eigvals.x () <= zeroThreshold;
1035+ eig_y_0 = eigvals.y () <= zeroThreshold;
1036+ eig_z_0 = eigvals.z () <= zeroThreshold;
1037+ // Check for zero eigenvectors
1038+ if (eig_x_0)
1039+ {
1040+ zero_norm = Cross (eigvectors[1 ], eigvectors[2 ]);
1041+ eigvectors[0 ] = zero_norm;
1042+ }
1043+ else if (eig_y_0)
1044+ {
1045+ zero_norm = Cross (eigvectors[0 ], eigvectors[2 ]);
1046+ eigvectors[1 ] = zero_norm;
1047+ }
1048+ else if (eig_z_0)
1049+ {
1050+ zero_norm = Cross (eigvectors[0 ], eigvectors[1 ]);
1051+ eigvectors[2 ] = zero_norm;
1052+ }
1053+ }
1054+
1055+ Transform rotate (Point (0.0 , 0.0 , 0.0 ), eigvectors[0 ], eigvectors[1 ], eigvectors[2 ]);
1056+ Transform trans = rotate;
1057+ trans.pre_translate ((Vector) center);
1058+
1059+ trans.post_scale (Vector (1.0 ,1.0 ,1.0 ) * eigvals);
1060+ rotate.post_scale (Vector (1.0 ,1.0 ,1.0 ) / eigvals);
1061+
1062+ int nu = resolution + 1 ;
1063+ int nv = resolution;
1064+
1065+ SinCosTable tab1 (nu, 0 , 2 * M_PI);
1066+ SinCosTable tab2 (nv, 0 , M_PI);
1067+
1068+ double cl = (eigvals[0 ] - eigvals[1 ]) / (eigvals[0 ] + eigvals[1 ] + eigvals[2 ]);
1069+ double cp = 2.0 * (eigvals[1 ] - eigvals[2 ]) / (eigvals[0 ] + eigvals[1 ] + eigvals[2 ]);
1070+ double A = spow ((1.0 - cl), emphasis);
1071+ double B = spow ((1.0 - cp), emphasis);
1072+
1073+ double nr[2 ];
1074+ double nz[2 ];
1075+
1076+ for (int v=0 ; v < nv-1 ; v++)
1077+ {
1078+ nr[0 ] = tab2.sin (v+1 );
1079+ nr[1 ] = tab2.sin (v);
1080+
1081+ nz[0 ] = tab2.cos (v+1 );
1082+ nz[1 ] = tab2.cos (v);
1083+
1084+ for (int u=0 ; u<nu; u++)
1085+ {
1086+ double nx = tab1.sin (u);
1087+ double ny = tab1.cos (u);
1088+
1089+ uint32_t offset = static_cast <uint32_t >(numVBOElements_);
1090+ for ( unsigned int i=0 ; i<2 ; i++ )
1091+ {
1092+ // Transorm points and add to points list
1093+ const double x = spow (nr[i], B) * spow (nx, A);
1094+ const double y = spow (nr[i], B) * spow (ny, A);
1095+ const double z = spow (nz[i], B);
1096+ Vector point = Vector (trans * Point ( x, y, z ));
1097+ points_.push_back (point);
1098+
1099+ // Add normals
1100+ const float nnx = spow (nr[i], 2.0 -B) * spow (nx, 2.0 -A);
1101+ const float nny = spow (nr[i], 2.0 -B) * spow (ny, 2.0 -A);
1102+ const float nnz = spow (nz[i], 2.0 -B);
1103+ Vector normal = rotate * Vector ( nnx, nny, nnz );
1104+ normal.safe_normalize ();
1105+ normals_.push_back (normal);
1106+
1107+ // Add color vectors from parameters
1108+ colors_.push_back (color);
1109+
1110+ numVBOElements_ ++;
1111+ }
1112+
1113+ indices_.push_back (0 + offset);
1114+ indices_.push_back (1 + offset);
1115+ indices_.push_back (2 + offset);
1116+ indices_.push_back (2 + offset);
1117+ indices_.push_back (1 + offset);
1118+ indices_.push_back (3 + offset);
1119+ }
1120+ }
1121+ for (int jj = 0 ; jj < 6 ; jj++) indices_.pop_back ();
1122+ }
1123+
1124+ // template <class T>
1125+ // void GeomGlyph::gen_torus(const Point& center, const T& t,
1126+ // double major_radius, double minor_radius,
1127+ // int nu, int nv,
1128+ // std::vector< QuadStrip >& quadstrips )
1129+ // {
1130+ // nu++; //Bring nu to expected value for shape.
1131+
1132+ // SinCosTable tab1(nu, 0, 2*M_PI);
1133+ // SinCosTable tab2(nv, 0, 2*M_PI, minor_radius);
1134+
1135+ // Transform trans;
1136+ // Transform rotate;
1137+ // gen_transforms( center, t, trans, rotate );
1138+
1139+ // Draw the torus
1140+ // for (int v=0; v<nv-1; v++)
1141+ // {
1142+ // double z1 = tab2.cos(v+1);
1143+ // double z2 = tab2.cos(v);
1144+
1145+ // double nr1 = tab2.sin(v+1);
1146+ // double nr2 = tab2.sin(v);
1147+
1148+ // double r1 = major_radius + nr1;
1149+ // double r2 = major_radius + nr2;
1150+
1151+ // QuadStrip quadstrip;
1152+
1153+ // for (int u=0; u<nu; u++)
1154+ // {
1155+ // double nx = tab1.sin(u);
1156+ // double ny = tab1.cos(u);
1157+
1158+ // double x1 = r1 * nx;
1159+ // double y1 = r1 * ny;
1160+
1161+ // double x2 = r2 * nx;
1162+ // double y2 = r2 * ny;
1163+
1164+ // Point p1 = trans * Point(x1, y1, z1);
1165+ // Point p2 = trans * Point(x2, y2, z2);
1166+
1167+ // Vector v1 = rotate * Vector(nr1*nx, nr1*ny, z1);
1168+ // Vector v2 = rotate * Vector(nr2*nx, nr2*ny, z2);
1169+
1170+ // v1.safe_normalize();
1171+ // v2.safe_normalize();
1172+
1173+ // quadstrip.push_back( std::make_pair(p1, v1) );
1174+ // quadstrip.push_back( std::make_pair(p2, v2) );
1175+ // }
1176+
1177+ // quadstrips.push_back( quadstrip );
1178+ // }
1179+ // }
1180+
9831181void GlyphGeom::generateLine (const Point& p1, const Point& p2, const ColorRGB& color1, const ColorRGB& color2)
9841182{
9851183 points_.push_back (Vector (p1));
0 commit comments