diff --git a/modules/geom/src/quat.cc b/modules/geom/src/quat.cc
index ceb178713e8641862d763a06f09023faa039004b..14aa5e60c55934c22930a433d707baea4eea6543 100644
--- a/modules/geom/src/quat.cc
+++ b/modules/geom/src/quat.cc
@@ -404,6 +404,54 @@ Quat Slerp(const Quat& q0, const Quat& q1, Real t)
return nrvo;
}
+Quat Inv(const Quat& q)
+{
+ Real l=q.w*q.w+q.x*q.x+q.y*q.y+q.z*q.z;
+ Quat ret;
+ if(l>0.0) {
+ ret.w=q.w/l;
+ ret.x=-q.x/l;
+ ret.y=-q.y/l;
+ ret.z=-q.z/l;
+ }
+ return ret;
+}
+
+Quat Exp(const Quat& q)
+{
+ Real a = sqrt(q.x*q.x+q.y*q.y+q.z*q.z);
+ Real sina = sin(a);
+ Real cosa = cos(a);
+ Quat ret;
+
+ ret.w = cosa;
+ if(a>0.0) {
+ ret.x = sina * q.x/a;
+ ret.y = sina * q.y/a;
+ ret.z = sina * q.z/a;
+ } else {
+ ret.x = ret.y = ret.z = 0.0;
+ }
+ return ret;
+}
+
+Quat Log(const Quat& q)
+{
+ Real a = acos(q.w);
+ Real sina = sin(a);
+ Quat ret;
+
+ ret.w = 0.0;
+ if(sina>0.0) {
+ ret.x = a*q.x/sina;
+ ret.y = a*q.y/sina;
+ ret.z = a*q.z/sina;
+ } else {
+ ret.x= ret.y= ret.z= 0.0;
+ }
+ return ret;
+}
+
Vec3 Quat::Rotate(const Vec3& vec) const {
Quat tmp(0.0, vec[0], vec[1], vec[2]);
Quat conj=Conjugate(*this);
diff --git a/modules/geom/src/quat.hh b/modules/geom/src/quat.hh
index f8babf8668e8763c5b98b42d23abaaf99317919f..d2d1ea3e45c2841df0d402ccba5675e6c894a49d 100644
--- a/modules/geom/src/quat.hh
+++ b/modules/geom/src/quat.hh
@@ -111,6 +111,10 @@ Real DLLEXPORT_OST_GEOM Dot(const Quat& q0, const Quat& q1);
// spherical linear interpolation, with t in range [0,1]
Quat DLLEXPORT_OST_GEOM Slerp(const Quat& q0, const Quat& q1, Real t);
+Quat DLLEXPORT_OST_GEOM Inv(const Quat& q);
+Quat DLLEXPORT_OST_GEOM Exp(const Quat& q);
+Quat DLLEXPORT_OST_GEOM Log(const Quat& q);
+
Quat DLLEXPORT_OST_GEOM Grassmann(const Quat& lhs, const Quat& rhs);
//normalize quaternion
diff --git a/modules/gfx/src/gfx_object.cc b/modules/gfx/src/gfx_object.cc
index e151ccd9e0bdaa031f3fbf4eda62e24bb6a27c30..c6eec4a3a1ff878ac66821a27c6d0395d043a453 100644
--- a/modules/gfx/src/gfx_object.cc
+++ b/modules/gfx/src/gfx_object.cc
@@ -603,6 +603,7 @@ void GfxObj::CleanColorOps(){
void GfxObj::Debug(unsigned int flags)
{
+ debug_flags_=flags;
va_.DrawNormals(flags & 0x1);
}
diff --git a/modules/gfx/src/impl/backbone_trace.cc b/modules/gfx/src/impl/backbone_trace.cc
index 2a0315380b75dfe55ef774155d60416c59f07fe9..c99508826de934526fcb3ee659a3f41bf9cae09f 100644
--- a/modules/gfx/src/impl/backbone_trace.cc
+++ b/modules/gfx/src/impl/backbone_trace.cc
@@ -185,13 +185,18 @@ void BackboneTrace::PrepList(NodeEntryList& nelist)
geom::Vec3 p12 = p2-p1;
if(p10==-p12 || p10==p12) p12+=geom::Vec3(0.001,0.001,0.001);
e1->v1=e1->normal;
+ e1->normal=geom::Normalize(geom::Cross(p10,p12));
+
// twist avoidance
- if(geom::Dot(e0->v1,e1->v1)<0.0) {
+ if(geom::Dot(geom::Normalize(geom::Cross(nref,p10)),
+ geom::Normalize(geom::Cross(p10,e1->normal)))>0.0) {
e1->v1=-e1->v1;
+ e1->normal=-e1->normal;
}
- e1->normal=geom::Normalize(geom::Cross(p10,p12));
+ nref = e1->normal;
+
float omega=0.5*acos(geom::Dot(geom::Normalize(p10),geom::Normalize(p12)));
- geom::Vec3 orth=geom::AxisRotation(e1->normal, -omega)*p12;
+ geom::Vec3 orth=geom::Normalize(geom::AxisRotation(e1->normal, -omega)*p12);
e1->direction=geom::Normalize(geom::Cross(e1->normal,orth));
// align normals to avoid twisting
diff --git a/modules/gfx/src/impl/cartoon_renderer.cc b/modules/gfx/src/impl/cartoon_renderer.cc
index b52d496f0c32a1b3008098b0c14a599cab98e502..0f9b5bb5a87327e61e05ca75e38e8f009c3b199b 100644
--- a/modules/gfx/src/impl/cartoon_renderer.cc
+++ b/modules/gfx/src/impl/cartoon_renderer.cc
@@ -108,7 +108,7 @@ void CartoonRenderer::PrepareRendering(const BackboneTrace& subset,
#if !defined(NDEBUG)
unsigned int tmp_count=0;
#endif
- for(SplineEntryListList::const_iterator sit=tmp_sll.begin();sit!=tmp_sll.end();++sit) {
+ for(SplineEntryListList::iterator sit=tmp_sll.begin();sit!=tmp_sll.end();++sit) {
if(sit->size()==2 and sit->at(0).type==6) {
// don't intpol cylinders
spline_list_list.push_back(*sit);
diff --git a/modules/gfx/src/impl/entity_detail.cc b/modules/gfx/src/impl/entity_detail.cc
index 40ba7351876fbf25ea823b383fdaca532fef41d6..d449fb0f66b5a0ddb96cd548109d28df5eaf2f5d 100644
--- a/modules/gfx/src/impl/entity_detail.cc
+++ b/modules/gfx/src/impl/entity_detail.cc
@@ -27,6 +27,7 @@
namespace ost {
using namespace mol;
+using namespace geom; // this is used _so_ often, its worth pulling it in
namespace gfx { namespace impl {
@@ -37,7 +38,7 @@ static const float default_radius=0.28;
struct BlurQuadEntry
{
float zdist;
- geom::Vec3 p1,p2,p3,p4;
+ Vec3 p1,p2,p3,p4;
gfx::Color c1,c2,c3,c4;
};
@@ -63,17 +64,17 @@ void DoRenderBlur(BondEntryList& bl, float bf1, float bf2)
if(!it->atom1 || !it->atom2) continue;
- const geom::Vec3 p0=tf.Apply(it->atom1->atom.GetPos());
- const geom::Vec3 p2=tf.Apply(it->atom2->atom.GetPos());
- geom::Vec3 p1=(p0+p2)*0.5;
+ const Vec3 p0=tf.Apply(it->atom1->atom.GetPos());
+ const Vec3 p2=tf.Apply(it->atom2->atom.GetPos());
+ Vec3 p1=(p0+p2)*0.5;
- const geom::Vec3 q0=tf.Apply(it->pp1);
- const geom::Vec3 q2=tf.Apply(it->pp2);
- geom::Vec3 q1=(q0+q2)*0.5;
+ const Vec3 q0=tf.Apply(it->pp1);
+ const Vec3 q2=tf.Apply(it->pp2);
+ Vec3 q1=(q0+q2)*0.5;
- float ll0 = geom::Length2(p0-q0);
- float ll1 = geom::Length2(p1-q1);
- float ll2 = geom::Length2(p2-q2);
+ float ll0 = Length2(p0-q0);
+ float ll1 = Length2(p1-q1);
+ float ll2 = Length2(p2-q2);
if(ll0<1e-2 && ll1<1e-2 && ll2<1e-2) continue;
@@ -240,7 +241,7 @@ static int bsplineGet(float *xa, float *ya, float *y2a, int n, float x, float *y
sublist.at((size-1)*nsub). COMPONENT = v; \
}
-SplineEntryList Spline::Generate(const SplineEntryList& entry_list, int nsub, uint color_blend_mode)
+SplineEntryList Spline::Generate(SplineEntryList& entry_list, int nsub, uint color_blend_mode)
{
if(nsub<=0) {
return entry_list;
@@ -304,44 +305,139 @@ SplineEntryList Spline::Generate(const SplineEntryList& entry_list, int nsub, ui
SPLINE_ENTRY_INTERPOLATE(rad);
LOG_DEBUG("SplineGenerate: assigning direction and normal components");
+
+#if 1
+ std::vector<Quat> quats(size+3);
+
+ Vec3 prevn;
+ for(int c=0;c<size;++c) {
+ Vec3 d1,d2;
+ if(c==0) {
+ d1=Normalize(entry_list.at(c).position-entry_list.at(c+1).position);
+ d2=Normalize(entry_list.at(c+1).position-entry_list.at(c+2).position);
+ } else if(c==size-1) {
+ d1=Normalize(entry_list.at(c-2).position-entry_list.at(c-1).position);
+ d2=Normalize(entry_list.at(c-1).position-entry_list.at(c).position);
+ } else {
+ d1=Normalize(entry_list.at(c-1).position-entry_list.at(c).position);
+ d2=Normalize(entry_list.at(c).position-entry_list.at(c+1).position);
+ }
+ Vec3 d=d2;
+ Vec3 n=Normalize(Cross(d1,d2));
+ //Vec3 n=entry_list.at(c).v1;
+ Vec3 o = Normalize(Cross(d,n));
+
+ if(c>0) {
+ if(Dot(Normalize(Cross(d,prevn)),Normalize(Cross(d,n)))<0.0) {
+ n=-n;
+ o=-o;
+ }
+ }
+ prevn=n;
+
+ //n = Normalize(Cross(o,d));
+ Quat q(Mat3(d[0],n[0],o[0],
+ d[1],n[1],o[1],
+ d[2],n[2],o[2]));
+ quats[c+1]=q;
+ }
+ quats.at(0)=quats.at(1);
+ quats.at(size)=quats.at(size-1);
+ quats.at(size+1)=quats.at(size-1);
+
+ std::vector<Quat> squats(size*nsub);
+ // now for the quaternion interpolation, using squad
+ for(int c=0;c<size;++c) {
+ Quat q0=quats.at(c+0);
+ Quat q1=quats.at(c+1);
+ Quat q2=quats.at(c+2);
+ Quat q3=quats.at(c+3);
+
+ Quat s1=q1*Exp(-0.25*Log(Inv(q1)*q2)+Log(Inv(q1)*q0));
+ Quat s2=q2*Exp(-0.25*Log(Inv(q2)*q3)+Log(Inv(q2)*q1));
+
+ for(int d=0;d<nsub;++d) {
+ float u=static_cast<float>(d)*i_nsub;
+ squats.at(c*nsub+d) = Slerp(Slerp(q1,q2,u),Slerp(s1,s2,u),2.0*u*(1.0-u));
+ }
+ }
+
+ assert(sublist.size()<=squats.size());
+
+ for(unsigned int i=0;i<sublist.size();++i) {
+ Mat3 m=squats[i].ToRotationMatrix();
+ sublist.at(i).direction=Normalize(Vec3(m(0,0),m(1,0),m(2,0)));
+ sublist.at(i).normal=Normalize(Vec3(m(0,1),m(1,1),m(2,1)));
+ sublist.at(i).v2=Normalize(Vec3(m(0,2),m(1,2),m(2,2)));
+ }
+
+ std::vector<Vec3> nlist(sublist.size());
+ for(uint i=0;i<sublist.size();++i) {
+ Vec3 p0,p1;
+ if(i<sublist.size()-1) {
+ p0 = sublist.at(i).position;
+ p1 = sublist.at(i+1).position;
+ } else {
+ p0 = sublist.at(i-1).position;
+ p1 = sublist.at(i).position;
+ }
+ sublist.at(i).direction=Normalize(p1-p0);
+ Vec3 orth = Normalize(Cross(sublist.at(i).direction,sublist.at(i).normal));
+ nlist[i] = Normalize(Cross(orth,sublist.at(i).direction));
+ sublist.at(i).v2=orth;
+ if(i>0) {
+ if(Dot(Normalize(Cross(sublist.at(i).direction,nlist[i-1])),orth)<0.0) {
+ nlist[i]=-nlist[i];
+ sublist.at(i).v2=-sublist.at(i).v2;
+ }
+ }
+ }
+
+ sublist.at(0).normal=nlist[0];
+ sublist.back().normal=nlist.back();
+ for(uint i=1;i<sublist.size()-1;++i) {
+ sublist.at(i).normal=(nlist[i-1]+nlist[i]+nlist[i+1])/3.0;
+ }
+#else
+
// assign direction and normal
// entity trace has the same algorithm
- geom::Vec3 p0 = sublist.at(0).position;
- geom::Vec3 p1 = sublist.at(1).position;
- geom::Vec3 p2 = ipsize>2 ? sublist.at(2).position : p1+(p1-p0);
+ Vec3 p0 = sublist.at(0).position;
+ Vec3 p1 = sublist.at(1).position;
+ Vec3 p2 = ipsize>2 ? sublist.at(2).position : p1+(p1-p0);
// normal of 0 is set at the end
- sublist.at(0).direction=geom::Normalize(p1-p0);
- sublist.at(0).v1=geom::Normalize(sublist.at(0).v1);
- geom::Vec3 orth = geom::Cross(sublist.at(0).direction,sublist.at(0).v1);
- sublist.at(0).v0 = geom::Normalize(geom::Cross(orth,sublist.at(0).direction));
+ sublist.at(0).direction=Normalize(p1-p0);
+ sublist.at(0).v1=Normalize(sublist.at(0).v1);
+ Vec3 orth = Cross(sublist.at(0).direction,sublist.at(0).v1);
+ sublist.at(0).v0 = Normalize(Cross(orth,sublist.at(0).direction));
// reference normal to avoid twisting
- //geom::Vec3 nref=geom::Normalize(geom::Cross(p0-p1,p2-p1));
+ //Vec3 nref=Normalize(Cross(p0-p1,p2-p1));
unsigned int i=1;
for(;i<sublist.size()-1;++i) {
- geom::Vec3 p10 = p0-p1;
- geom::Vec3 p12 = p2-p1;
+ Vec3 p10 = p0-p1;
+ Vec3 p12 = p2-p1;
// correction for perfectly aligned consecutive directions
- if(p10==-p12 || p10==p12) p12+=geom::Vec3(0.001,0.001,0.001);
- sublist.at(i).normal=geom::Normalize(geom::Cross(p10,p12));
+ if(p10==-p12 || p10==p12) p12+=Vec3(0.001,0.001,0.001);
+ sublist.at(i).normal=Normalize(Cross(p10,p12));
// paranoid error checking due to occasional roundoff troubles
- float cosw = geom::Dot(geom::Normalize(p10),geom::Normalize(p12));
+ float cosw = Dot(Normalize(p10),Normalize(p12));
cosw = std::min(float(1.0),std::max(float(-1.0),cosw));
float omega=0.5*acos(cosw);
- orth=geom::AxisRotation(sublist.at(i).normal, -omega)*p12;
- sublist.at(i).direction=geom::Normalize(geom::Cross(sublist.at(i).normal,orth));
+ orth=AxisRotation(sublist.at(i).normal, -omega)*p12;
+ sublist.at(i).direction=Normalize(Cross(sublist.at(i).normal,orth));
// twist avoidance
- sublist.at(i).v1=geom::Normalize(sublist.at(i).v1);
- orth = geom::Cross(sublist.at(i).direction,sublist.at(i).v1);
- sublist.at(i).v0 = geom::Normalize(geom::Cross(orth,sublist.at(i).direction));
- if(geom::Dot(sublist.at(i-1).v0,sublist.at(i).v0)<0.0) {
+ sublist.at(i).v1=Normalize(sublist.at(i).v1);
+ orth = Cross(sublist.at(i).direction,sublist.at(i).v1);
+ sublist.at(i).v0 = Normalize(Cross(orth,sublist.at(i).direction));
+ if(Dot(sublist.at(i-1).v0,sublist.at(i).v0)<0.0) {
sublist.at(i).v0=-sublist.at(i).v0;
//sublist.at(i).nflip = !sublist.at(i).nflip;
}
// avoid twisting
- //if(geom::Dot(sublist.at(i).normal,nref)<0.0) sublist.at(i).normal=-sublist.at(i).normal;
+ //if(Dot(sublist.at(i).normal,nref)<0.0) sublist.at(i).normal=-sublist.at(i).normal;
//nref=sublist.at(i).normal;
// skip over shift for the last iteration
if(i==sublist.size()-2) break;
@@ -352,12 +448,12 @@ SplineEntryList Spline::Generate(const SplineEntryList& entry_list, int nsub, ui
}
// assign remaining ones
sublist.at(0).normal=sublist.at(1).normal;
- sublist.at(i+1).direction=geom::Normalize(p2-p1);
+ sublist.at(i+1).direction=Normalize(p2-p1);
sublist.at(i+1).normal=sublist.at(i).normal;
- sublist.at(i+1).v1=geom::Normalize(sublist.at(i+1).v1);
- orth = geom::Cross(sublist.at(i+1).direction,sublist.at(i+1).v1);
- sublist.at(i+1).v0 = geom::Normalize(geom::Cross(orth,sublist.at(i+1).direction));
- if(geom::Dot(sublist.at(i).v0,sublist.at(i+1).v0)<0.0) {
+ sublist.at(i+1).v1=Normalize(sublist.at(i+1).v1);
+ orth = Cross(sublist.at(i+1).direction,sublist.at(i+1).v1);
+ sublist.at(i+1).v0 = Normalize(Cross(orth,sublist.at(i+1).direction));
+ if(Dot(sublist.at(i).v0,sublist.at(i+1).v0)<0.0) {
sublist.at(i+1).v0=-sublist.at(i+1).v0;
}
@@ -367,6 +463,8 @@ SplineEntryList Spline::Generate(const SplineEntryList& entry_list, int nsub, ui
sublist.at(i).normal = sublist.at(i).v0;
}
+#endif
+
LOG_DEBUG("SplineGenerate: assigning non-interpolated entry components");
// finally the non-interpolated type
// with some tweaks for proper strand rendering
@@ -414,18 +512,18 @@ SplineEntryList Spline::Generate(const SplineEntryList& entry_list, int nsub, ui
while(c<sublist.size()-1) {
int cstart=c;
if(sublist.at(c).type==1 && sublist.at(c+1).type==1) {
- geom::Vec3 n = geom::Normalize(geom::Cross(sublist.at(c).normal,
+ Vec3 n = Normalize(Cross(sublist.at(c).normal,
sublist.at(c).direction));
- geom::Vec3 p0 = sublist.at(c).position+n;
- geom::Vec3 q0 = sublist.at(c).position-n;
+ Vec3 p0 = sublist.at(c).position+n;
+ Vec3 q0 = sublist.at(c).position-n;
float psum=0.0;
float qsum=0.0;
++c;
while(c<sublist.size() && sublist.at(c).type==1) {
- n = geom::Normalize(geom::Cross(sublist.at(c).normal,
+ n = Normalize(Cross(sublist.at(c).normal,
sublist.at(c).direction));
- geom::Vec3 p1 = sublist.at(c).position+n;
- geom::Vec3 q1 = sublist.at(c).position-n;
+ Vec3 p1 = sublist.at(c).position+n;
+ Vec3 q1 = sublist.at(c).position-n;
psum+=Length(p1-p0);
qsum+=Length(q1-q0);
p0=p1;
diff --git a/modules/gfx/src/impl/entity_detail.hh b/modules/gfx/src/impl/entity_detail.hh
index 2ae0a848f14c3e7c711f67dc31e436969a4a22d0..c41a4a98cbe1592c7e78041df0f533c8cecfbba1 100644
--- a/modules/gfx/src/impl/entity_detail.hh
+++ b/modules/gfx/src/impl/entity_detail.hh
@@ -158,7 +158,7 @@ typedef std::vector<SplineEntryList> SplineEntryListList;
class DLLEXPORT_OST_GFX Spline {
public:
- static SplineEntryList Generate(const SplineEntryList& entry_list,int nsub,uint color_blend_mode=0);
+ static SplineEntryList Generate(SplineEntryList& entry_list,int nsub,uint color_blend_mode=0);
};
}}} // ns
diff --git a/modules/gfx/src/impl/sline_renderer.cc b/modules/gfx/src/impl/sline_renderer.cc
index 62ef7ba45aa756d0ce5c4f45d4e4b22c0418b78a..92cf81bd54184bb5a005c5d514a71b786f9b0ddd 100644
--- a/modules/gfx/src/impl/sline_renderer.cc
+++ b/modules/gfx/src/impl/sline_renderer.cc
@@ -81,8 +81,14 @@ void SlineRenderer::PrepareRendering(const BackboneTrace& trace_subset,
VertexID p0=va.Add(sit->position, geom::Vec3(),sit->color1);
for (++sit; sit!=sel.end(); ++sit) {
VertexID p1 = va.Add(sit->position, geom::Vec3(),sit->color1);
- va.AddLine(p0,p1);
+ //va.AddLine(p0,p1);
p0=p1;
+ VertexID p2 = va.Add(sit->position+sit->direction,geom::Vec3(),Color(0,1,0));
+ VertexID p3 = va.Add(sit->position+sit->normal,geom::Vec3(),Color(1,0,0));
+ VertexID p4 = va.Add(sit->position+sit->v2,geom::Vec3(),Color(1,0,1));
+ va.AddLine(p0,p2);
+ va.AddLine(p0,p3);
+ va.AddLine(p0,p4);
}
}
}
diff --git a/modules/img/alg/pymod/export_transcendentals.cc b/modules/img/alg/pymod/export_transcendentals.cc
index 2a9f42cc72add271b1f37ad800470b5cb6eaacd6..e2f22e8c13a20af0b3342f9caf3ee2e737201e0b 100644
--- a/modules/img/alg/pymod/export_transcendentals.cc
+++ b/modules/img/alg/pymod/export_transcendentals.cc
@@ -28,16 +28,16 @@ using namespace boost::python;
#include <ost/img/alg/transcendentals.hh>
using namespace ost::img;
-using namespace ost::img::alg;
+//using namespace ost::img::alg;
void export_Transcendentals()
{
- class_<Cos, bases<ConstModIPAlgorithm> >("Cos", init<>());
- class_<Exp, bases<ConstModIPAlgorithm> >("Exp", init<>());
- class_<Log, bases<ConstModIPAlgorithm> >("Log", init<>());
- class_<Log10, bases<ConstModIPAlgorithm> >("Log10", init<>());
- class_<Sin, bases<ConstModIPAlgorithm> >("Sin", init<>());
- class_<Sqrt, bases<ConstModIPAlgorithm> >("Sqrt", init<>());
- class_<Tan, bases<ConstModIPAlgorithm> >("Tan", init<>());
- class_<Pow, bases<ConstModIPAlgorithm> >("Pow", init<Real>());
+ class_<alg::Cos, bases<ConstModIPAlgorithm> >("Cos", init<>());
+ class_<alg::Exp, bases<ConstModIPAlgorithm> >("Exp", init<>());
+ class_<alg::Log, bases<ConstModIPAlgorithm> >("Log", init<>());
+ class_<alg::Log10, bases<ConstModIPAlgorithm> >("Log10", init<>());
+ class_<alg::Sin, bases<ConstModIPAlgorithm> >("Sin", init<>());
+ class_<alg::Sqrt, bases<ConstModIPAlgorithm> >("Sqrt", init<>());
+ class_<alg::Tan, bases<ConstModIPAlgorithm> >("Tan", init<>());
+ class_<alg::Pow, bases<ConstModIPAlgorithm> >("Pow", init<Real>());
}
diff --git a/modules/img/base/src/vecmat.hh b/modules/img/base/src/vecmat.hh
index 0fe9d299dbd292e8365aa8e41a45292fcc72d8e4..d2c82ba9055e44dbba5c934621fe6c9b2cbe25e3 100644
--- a/modules/img/base/src/vecmat.hh
+++ b/modules/img/base/src/vecmat.hh
@@ -37,6 +37,7 @@
namespace ost { namespace img {
// this pulls the geom namespace into the img one
+// which is a bad idea
using namespace ::geom;