diff --git a/modules/geom/src/vecmat3_op.cc b/modules/geom/src/vecmat3_op.cc index 4194cac02880ea0d4fb7d3cd1fd528c247dacec9..3066deef9a85a45d42fc6a77641d673667a7e89f 100644 --- a/modules/geom/src/vecmat3_op.cc +++ b/modules/geom/src/vecmat3_op.cc @@ -101,6 +101,14 @@ Real Angle(const Vec3& v1, const Vec3& v2) return std::acos(dot_product); } +Real SignedAngle(const Vec3& v1, const Vec3& v2, const Vec3& ref_normal ) +{ + Vec3 a=Normalize(v1); + Vec3 b=Normalize(v2); + Vec3 c=Normalize(ref_normal); + return std::atan2(Dot(c,Cross(a,b)), Dot(a,b)); +} + Mat3 EulerTransformation(Real theta, Real phi, Real xi) { Real costheta=cos(theta); diff --git a/modules/geom/src/vecmat3_op.hh b/modules/geom/src/vecmat3_op.hh index 76bd719d282ffc4261fdcc15d3104c9f2e4ff079..94c3c77b32d2fe7cb45a34d0a66c5bc001a9a3fe 100644 --- a/modules/geom/src/vecmat3_op.hh +++ b/modules/geom/src/vecmat3_op.hh @@ -151,9 +151,12 @@ Real DLLEXPORT_OST_GEOM Minor(const Mat3& m, unsigned int i, unsigned int j); // determinant Real DLLEXPORT_OST_GEOM Det(const Mat3& m); -// angle between zwo vectors +// angle between two vectors Real DLLEXPORT_OST_GEOM Angle(const Vec3& v1, const Vec3& v2); +// signed angle between two vectors, based on a reference normal +Real DLLEXPORT_OST_GEOM SignedAngle(const Vec3& v1, const Vec3& v2, const Vec3& ref); + Mat3 DLLEXPORT_OST_GEOM EulerTransformation(Real theta, Real phi, Real xi); Mat3 DLLEXPORT_OST_GEOM AxisRotation(const Vec3& axis, Real angle); diff --git a/modules/gfx/src/impl/backbone_trace.cc b/modules/gfx/src/impl/backbone_trace.cc index e3a14b7a7106414d3fb98cdfbab962b8a0a34c95..235ad226b232465e037342333493de62063343b9 100644 --- a/modules/gfx/src/impl/backbone_trace.cc +++ b/modules/gfx/src/impl/backbone_trace.cc @@ -83,7 +83,7 @@ public: NodeEntry entry={ca, GfxObj::Ele2Color(ca.GetElement()), GfxObj::Ele2Color(ca.GetElement()), geom::Vec3(), // this will be set by the gfx trace obj - res.GetCentralNormal(), + geom::Normalize(res.GetCentralNormal()), 1.0, geom::Vec3(),geom::Vec3(),geom::Vec3(), // for later use in NA rendering false,id_counter_++}; @@ -159,46 +159,36 @@ void BackboneTrace::AddNodeEntryList(const NodeEntryList& l) void BackboneTrace::PrepList(NodeEntryList& nelist) { - // assign direction and normal vectors for each entry - // they are composed of the i-1 -> i and i->i+1 directions - // - // this same algorithm is used in the spline generation, so - // perhaps all of this here is not necessary ?! - // + // orthogonalize the residue normals with + // twist detection; important for later + // spline generation + if(nelist.size()<3) return; NodeEntry* e0=&nelist[0]; NodeEntry* e1=&nelist[1]; NodeEntry* e2=&nelist[2]; geom::Vec3 p0 = e0->atom.GetPos(); geom::Vec3 p1 = e1->atom.GetPos(); geom::Vec3 p2 = e2->atom.GetPos(); - + + geom::Vec3 dir=geom::Normalize(p1-p0); e0->direction=geom::Normalize(p1-p0); - // e0->normal is set afterwards to normal of second one - // backup residue normal - e0->v1 = e0->normal; - - //reference normal to avoid twisting - geom::Vec3 nref=geom::Normalize(geom::Cross(p0-p1,p2-p1)); - - // start loop with the second + geom::Vec3 orth=geom::Normalize(geom::Cross(e0->direction,e0->normal)); + geom::Vec3 norm=geom::Normalize(geom::Cross(orth,dir)); + e0->normal=norm; + // start loop with the second residue unsigned int i=1; for(;i<nelist.size()-1;++i) { - geom::Vec3 p10 = p0-p1; - geom::Vec3 p12 = p2-p1; - if(p10==-p12 || p10==p12) p12+=geom::Vec3(0.001,0.001,0.001); - e1->v1=e1->normal; - // twist avoidance - if(geom::Dot(e0->v1,e1->v1)<0.0) { - e1->v1=-e1->v1; + geom::Vec3 p10=p0-p1; + geom::Vec3 p12=p2-p1; + dir=geom::Normalize(p2-p0); + e1->direction = dir; + orth=geom::Normalize(geom::Cross(dir,e1->normal)); + norm=geom::Normalize(geom::Cross(orth,dir)); + // twist check + if(geom::Dot(geom::Cross(e0->normal,dir),geom::Cross(norm,dir))<0.0) { + norm=-norm; } - e1->normal=geom::Normalize(geom::Cross(p10,p12)); - float omega=0.5*acos(geom::Dot(geom::Normalize(p10),geom::Normalize(p12))); - geom::Vec3 orth=geom::AxisRotation(e1->normal, -omega)*p12; - e1->direction=geom::Normalize(geom::Cross(e1->normal,orth)); - - // align normals to avoid twisting - //if(geom::Dot(e1->normal,nref)<0.0) e1->normal=-e1->normal; - //nref=e1->normal; + e1->normal=norm; // skip over shift for the last iteration if(i==nelist.size()-2) break; // shift to i+1 for next iteration @@ -209,16 +199,14 @@ void BackboneTrace::PrepList(NodeEntryList& nelist) p1 = e1->atom.GetPos(); p2 = e2->atom.GetPos(); } - // finish remaining values - // i is at size-2 due to break statement above - nelist[0].normal=nelist[1].normal; - nelist[i+1].direction=geom::Normalize(p2-p1); - nelist[i+1].v1=nelist[i+1].normal; - if (geom::Dot(nelist[i].v1, - nelist[i+1].v1)<0.0) { - nelist[i+1].v1=-nelist[i+1].v1; + dir=geom::Normalize(p2-p1); + e2->direction=dir; + orth=geom::Normalize(geom::Cross(dir,e2->normal)); + norm=geom::Normalize(geom::Cross(orth,dir)); + if(geom::Dot(geom::Cross(e1->normal,dir),geom::Cross(norm,dir))<0.0) { + norm=-norm; } - nelist[i+1].normal=nelist[i].normal; + e2->normal=norm; } BackboneTrace BackboneTrace::CreateSubset(const mol::EntityView& subview) diff --git a/modules/gfx/src/impl/cartoon_renderer.cc b/modules/gfx/src/impl/cartoon_renderer.cc index 1c5bc481bc2d574f62b35aa18aab155a3b02b01b..b1ec4dad4f0f0cbe5092a16a9f3e654c348b14fb 100644 --- a/modules/gfx/src/impl/cartoon_renderer.cc +++ b/modules/gfx/src/impl/cartoon_renderer.cc @@ -389,33 +389,33 @@ void CartoonRenderer::RebuildSplineObj(IndexedVertexArray& va, std::vector<TraceProfile> profiles; float factor=is_sel ? 0.2 : 0.0; profiles.push_back(GetCircProfile(detail, - options_->GetTubeRadius()*options_->GetTubeRatio()+factor, - options_->GetTubeRadius()+factor, - options_->GetTubeProfileType(), - 1.0)); // profile 0 = tube + options_->GetTubeRadius()*options_->GetTubeRatio()+factor, + options_->GetTubeRadius()+factor, + options_->GetTubeProfileType(), + 1.0)); // profile 0 = tube if (!force_tube_) { profiles.push_back(GetCircProfile(detail, - options_->GetHelixWidth()+factor, - options_->GetHelixThickness()+factor, - options_->GetHelixProfileType(), - options_->GetHelixEcc())); // profile 1 = helix + options_->GetHelixWidth()+factor, + options_->GetHelixThickness()+factor, + options_->GetHelixProfileType(), + options_->GetHelixEcc())); // profile 1 = helix profiles.push_back(GetCircProfile(detail, - options_->GetStrandWidth()+factor, - options_->GetStrandThickness()+factor, - options_->GetStrandProfileType(), - options_->GetStrandEcc())); // profile 2 = strand + options_->GetStrandWidth()+factor, + options_->GetStrandThickness()+factor, + options_->GetStrandProfileType(), + options_->GetStrandEcc())); // profile 2 = strand profiles.push_back(profiles.back()); // profile 3==2, strand - + profiles.push_back(GetCircProfile(detail, - 1.7*options_->GetStrandWidth()+factor, - 1.1*options_->GetStrandThickness()+factor, - options_->GetStrandProfileType(), - options_->GetStrandEcc())); // profile 4 = arrow start + 1.7*options_->GetStrandWidth()+factor, + 1.1*options_->GetStrandThickness()+factor, + options_->GetStrandProfileType(), + options_->GetStrandEcc())); // profile 4 = arrow start profiles.push_back(GetCircProfile(detail, - 0.01*options_->GetStrandWidth()+factor, - 1.1*options_->GetStrandThickness()+factor, - options_->GetStrandProfileType(), - options_->GetStrandEcc())); // profile 5 = arrow end + 0.01*options_->GetStrandWidth()+factor, + 1.1*options_->GetStrandThickness()+factor, + options_->GetStrandProfileType(), + options_->GetStrandEcc())); // profile 5 = arrow end } @@ -447,14 +447,23 @@ void CartoonRenderer::RebuildSplineObj(IndexedVertexArray& va, TraceProfile tprof1=TransformAndAddProfile(profiles,slist[0],va); CapProfile(tprof1,slist[0],true,va); TraceProfile tprof2; - unsigned int sc=1; + size_t sc=1; + size_t psize=tprof1.size()-1; + double psized=static_cast<double>(psize); + float of=static_cast<float>(psize)/(2.0*M_PI); for (; sc<slist.size(); ++sc) { + //double dd=geom::Dot(slist.at(sc-1).normal,slist.at(sc).normal); + geom::Vec3 dir=slist.at(sc).position-slist.at(sc-1).position; + double ang=-geom::SignedAngle(geom::Cross(slist.at(sc-1).normal,dir), + geom::Cross(slist.at(sc).normal,dir), + dir); + size_t offset=psize+static_cast<size_t>(round(ang*of)+psized); if(slist.at(sc).type==3) { if(slist.at(sc-1).type!=3) { // boundary to arrow SplineEntry se(slist[sc]); tprof2=TransformAndAddProfile(profiles,se, va); - AssembleProfile(tprof1,tprof2,va); + AssembleProfile(tprof1,tprof2,va,offset); tprof1=tprof2; se.type=2; se.type1=4; @@ -469,7 +478,7 @@ void CartoonRenderer::RebuildSplineObj(IndexedVertexArray& va, } else { tprof2=TransformAndAddProfile(profiles,slist.at(sc), va); } - AssembleProfile(tprof1,tprof2,va); + AssembleProfile(tprof1,tprof2,va,offset); tprof1=tprof2; } CapProfile(tprof1,slist.at(sc-1),false,va); @@ -551,49 +560,22 @@ namespace { void CartoonRenderer::AssembleProfile(const TraceProfile& prof1, const TraceProfile& prof2, - IndexedVertexArray& va) + IndexedVertexArray& va, + size_t offset) { /* the wrap around algorithm used here needs to take into account that the TraceProfile has a duplicate entry for prof*[0] and prof*[N-1], which fall onto the same point except and have the same normal but a different texture coordinate. Hence - all the mods are done with size()-1, but in the assembly - routine prof*[0] is turned into prof*[N-1] if necessary + size()-1 */ size_t size=prof1.size()-1; - - // first get the best correction offset - float accum[]={0.0,0.0,0.0,0.0,0.0}; - for(size_t i=0;i<size;++i) { - int i1=(i+0)%(size); - int i2=(i+1)%(size); - geom::Vec3 v1=va.GetVert(prof1[i1].id); - geom::Vec3 v2=va.GetVert(prof1[i2].id); - for(int k=-2;k<=2;++k) { - int i3=(i+k+0+size)%(size); - int i4=(i+k+1+size)%(size); - geom::Vec3 v3=va.GetVert(prof2[i3].id); - geom::Vec3 v4=va.GetVert(prof2[i4].id); - accum[k+2]+=spread(v1,v2,v3,v4); - } - } - - float best_spread=accum[0]; - int best_off=-2; - for(int k=-1;k<=2;++k) { - if(accum[k+2]<best_spread) { - best_spread=accum[k+2]; - best_off=k; - } - } - best_off=(best_off+(size))%(size); - // now assemble the triangles for(unsigned int i1=0;i1<size;++i1) { - unsigned int i2=(i1+1)%(size); - unsigned int i3=(i1+best_off)%(size); - unsigned int i4=(i1+best_off+1)%(size); + size_t i2=(i1+1)%(size); + size_t i3=(i1+offset)%(size); + size_t i4=(i1+offset+1)%(size); // wrap around correction for proper texture coordinates i2 = (i2==0) ? size : i2; diff --git a/modules/gfx/src/impl/cartoon_renderer.hh b/modules/gfx/src/impl/cartoon_renderer.hh index 95a7b3844802e14c9f3fd087b4f67b3021c47ab4..475c1db2b110055f73e0275bea73dcb5373114f9 100644 --- a/modules/gfx/src/impl/cartoon_renderer.hh +++ b/modules/gfx/src/impl/cartoon_renderer.hh @@ -66,7 +66,8 @@ private: void AssembleProfile(const TraceProfile& prof1, const TraceProfile& prof2, - IndexedVertexArray& va); + IndexedVertexArray& va, + size_t offset); TraceProfile TransformAndAddProfile(const std::vector<TraceProfile>& profiles, const SplineEntry& se, diff --git a/modules/gfx/src/impl/entity_detail.cc b/modules/gfx/src/impl/entity_detail.cc index 40ba7351876fbf25ea823b383fdaca532fef41d6..55fab9f53d8b7553dbb1d8bd4fa20ed3ea6327a2 100644 --- a/modules/gfx/src/impl/entity_detail.cc +++ b/modules/gfx/src/impl/entity_detail.cc @@ -296,76 +296,37 @@ SplineEntryList Spline::Generate(const SplineEntryList& entry_list, int nsub, ui // interpolate position and helper vectors for(int k=0;k<3;++k) { SPLINE_ENTRY_INTERPOLATE(position[k]); - //SPLINE_ENTRY_INTERPOLATE(v0[k]); - SPLINE_ENTRY_INTERPOLATE(v1[k]); - //SPLINE_ENTRY_INTERPOLATE(v2[k]); + SPLINE_ENTRY_INTERPOLATE(normal[k]); } SPLINE_ENTRY_INTERPOLATE(rad); LOG_DEBUG("SplineGenerate: assigning direction and normal components"); - // assign direction and normal - // entity trace has the same algorithm + // assign direction and then re-assign normal 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); - // normal of 0 is set at the end + geom::Vec3 p2 = sublist.size()>2 ? sublist.at(2).position : p1+p1-p0; 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)); - - // reference normal to avoid twisting - //geom::Vec3 nref=geom::Normalize(geom::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; - // 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)); - // paranoid error checking due to occasional roundoff troubles - float cosw = geom::Dot(geom::Normalize(p10),geom::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)); - // 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).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; - //nref=sublist.at(i).normal; - // skip over shift for the last iteration + geom::Vec3 dir=geom::Normalize(p2-p0); + sublist.at(i).direction=dir; + geom::Vec3 orth=geom::Normalize(geom::Cross(dir,sublist[i].normal)); + geom::Vec3 norm=geom::Normalize(geom::Cross(orth,dir)); + sublist[i].normal=norm; if(i==sublist.size()-2) break; // shift to i+1 for next iteration p0 = sublist.at(i).position; p1 = sublist.at(i+1).position; p2 = sublist.at(i+2).position; } - // assign remaining ones - sublist.at(0).normal=sublist.at(1).normal; - sublist.at(i+1).direction=geom::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).v0=-sublist.at(i+1).v0; - } + sublist[0].normal=sublist[0].normal; + sublist[i+1].direction=geom::Normalize(p2-p1); + sublist[i+1].normal=sublist[i].normal; - // hack - // TODO: merge this with above routine - for(unsigned int i=0;i<sublist.size();++i) { - sublist.at(i).normal = sublist.at(i).v0; - } LOG_DEBUG("SplineGenerate: assigning non-interpolated entry components"); // finally the non-interpolated type diff --git a/modules/gfx/src/impl/line_trace_renderer.cc b/modules/gfx/src/impl/line_trace_renderer.cc index f183c6b5b0439696f8fa357baaeb24951c65dc49..abf91fcad5c3bfafd27f5a093ea88951eccb9136 100644 --- a/modules/gfx/src/impl/line_trace_renderer.cc +++ b/modules/gfx/src/impl/line_trace_renderer.cc @@ -88,43 +88,51 @@ void LineTraceRenderer::PrepareRendering(const BackboneTrace& trace_subset, va.SetAALines(options_->GetAALines()); if(is_sel) { for (int node_list=0; node_list<trace_subset.GetListCount(); ++node_list) { - const NodeEntryList& nl=trace_subset.GetList(node_list); - if(nl.size()<1) continue; - for(unsigned int i=0;i<nl.size();++i) { - const NodeEntry& entry=nl[i]; - if(sel_.FindAtom(entry.atom).IsValid()) { - geom::Vec3 apos = entry.atom.GetPos(); - VertexID p0=va.Add(apos, geom::Vec3(),sel_clr); - if(i>0) { - VertexID p1 =va.Add(apos+0.5*(nl[i-1].atom.GetPos()-apos), geom::Vec3(), sel_clr); - va.AddLine(p0, p1); - } - if(i<nl.size()-1) { - VertexID p1 =va.Add(apos+0.5*(nl[i+1].atom.GetPos()-apos), geom::Vec3(), sel_clr); - va.AddLine(p0, p1); - } - } - } + const NodeEntryList& nl=trace_subset.GetList(node_list); + if(nl.size()<1) continue; + for(unsigned int i=0;i<nl.size();++i) { + const NodeEntry& entry=nl[i]; + if(sel_.FindAtom(entry.atom).IsValid()) { + geom::Vec3 apos = entry.atom.GetPos(); + VertexID p0=va.Add(apos, geom::Vec3(),sel_clr); + if(i>0) { + VertexID p1 =va.Add(apos+0.5*(nl[i-1].atom.GetPos()-apos), geom::Vec3(), sel_clr); + va.AddLine(p0, p1); + } + if(i<nl.size()-1) { + VertexID p1 =va.Add(apos+0.5*(nl[i+1].atom.GetPos()-apos), geom::Vec3(), sel_clr); + va.AddLine(p0, p1); + } + } + } } } else { for (int node_list=0; node_list<trace_subset.GetListCount(); ++node_list) { - const NodeEntryList& nl=trace_subset.GetList(node_list); - - if (nl.size()<2) continue; - - VertexID p0=va.Add(nl[0].atom.GetPos(), geom::Vec3(), - nl[0].color1); - for (unsigned int i=1; i<nl.size()-1;++i) { - const NodeEntry& entry=nl[i]; - VertexID p1 =va.Add(entry.atom.GetPos(), geom::Vec3(), - entry.color1); - va.AddLine(p0, p1); - p0=p1; - } - const NodeEntry& entry=nl.back(); - VertexID p1 =va.Add(entry.atom.GetPos(), geom::Vec3(), - entry.color1); - va.AddLine(p0, p1); + const NodeEntryList& nl=trace_subset.GetList(node_list); + + if (nl.size()<2) continue; + + VertexID p0=va.Add(nl[0].atom.GetPos(), geom::Vec3(), + nl[0].color1); + for (unsigned int i=1; i<nl.size()-1;++i) { + const NodeEntry& entry=nl[i]; + VertexID p1 =va.Add(entry.atom.GetPos(), geom::Vec3(), + entry.color1); + va.AddLine(p0, p1); + p0=p1; +#if 1 + VertexID p2 =va.Add(entry.atom.GetPos()+entry.direction, geom::Vec3(), + Color(1,0,0)); + VertexID p3 =va.Add(entry.atom.GetPos()+entry.normal, geom::Vec3(), + Color(0,1,1)); + va.AddLine(p0,p2); + va.AddLine(p0,p3); +#endif + } + const NodeEntry& entry=nl.back(); + VertexID p1 =va.Add(entry.atom.GetPos(), geom::Vec3(), + entry.color1); + va.AddLine(p0, p1); } } } diff --git a/modules/gfx/src/impl/sline_renderer.cc b/modules/gfx/src/impl/sline_renderer.cc index 2b54046dd634acc1da4413aa9d69b7f4960889a9..92aa0f4836b4b8d603d0279fd655c37206789652 100644 --- a/modules/gfx/src/impl/sline_renderer.cc +++ b/modules/gfx/src/impl/sline_renderer.cc @@ -72,7 +72,7 @@ void SlineRenderer::PrepareRendering(const BackboneTrace& trace_subset, entry.normal, entry.rad, is_sel ? sel_clr : entry.color1, is_sel ? sel_clr : entry.color2, - 0, entry.id); + 0, entry.id); ee.v1 = entry.v1; spl.push_back(ee); } @@ -83,13 +83,11 @@ void SlineRenderer::PrepareRendering(const BackboneTrace& trace_subset, VertexID p1 = va.Add(sit->position, geom::Vec3(),sit->color1); va.AddLine(p0,p1); p0=p1; -#if 0 - 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)); +#if 1 + VertexID p2 = va.Add(sit->position+sit->direction,geom::Vec3(),Color(1,0,0)); + VertexID p3 = va.Add(sit->position+sit->normal,geom::Vec3(),Color(0,1,1)); va.AddLine(p0,p2); va.AddLine(p0,p3); - va.AddLine(p0,p4); #endif } } diff --git a/modules/mol/base/src/impl/residue_impl.cc b/modules/mol/base/src/impl/residue_impl.cc index 75590b1fea096169c76226df15442543fafd697e..6d0856b111c071316e4604b7c6d9a18d1c953085 100644 --- a/modules/mol/base/src/impl/residue_impl.cc +++ b/modules/mol/base/src/impl/residue_impl.cc @@ -212,24 +212,28 @@ ChainImplPtr ResidueImpl::GetChain() const geom::Vec3 ResidueImpl::GetCentralNormal() const { + geom::Vec3 nrvo(1,0,0); if (chem_class_.IsPeptideLinking()) { AtomImplPtr a1 = FindAtom("C"); AtomImplPtr a2 = FindAtom("O"); - geom::Vec3 nrvo; if(a1 && a2) { - nrvo = geom::Normalize(a1->GetPos()-a2->GetPos()); + nrvo = geom::Normalize(a2->GetPos()-a1->GetPos()); } else { - geom::Vec3 v0=GetCentralAtom()->GetPos(); - nrvo=geom::Cross(geom::Normalize(v0), - geom::Normalize(geom::Vec3(-v0[2],v0[0],v0[1]))); - LOG_VERBOSE("warning: could not find atoms for proper central normal calculation"); + a1 = FindAtom("CB"); + a2 = FindAtom("CA"); + if(a1 && a2) { + nrvo = geom::Normalize(a2->GetPos()-a1->GetPos()); + } else { + geom::Vec3 v0=GetCentralAtom()->GetPos(); + nrvo=geom::Cross(geom::Normalize(v0), + geom::Normalize(geom::Vec3(-v0[2],v0[0],v0[1]))); + LOG_VERBOSE("warning: could not find atoms for proper central normal calculation"); + } } - return nrvo; } else if (chem_class_.IsNucleotideLinking()) { AtomImplPtr a1 = FindAtom("P"); AtomImplPtr a2 = FindAtom("OP1"); AtomImplPtr a3 = FindAtom("OP2"); - geom::Vec3 nrvo; if(a1 && a2 && a3) { nrvo = geom::Normalize(a1->GetPos()-(a2->GetPos()+a3->GetPos())*.5); } else { @@ -238,9 +242,8 @@ geom::Vec3 ResidueImpl::GetCentralNormal() const geom::Normalize(geom::Vec3(-v0[2],v0[0],v0[1]))); LOG_VERBOSE("warning: could not find atoms for proper central normal calculation"); } - return nrvo; } - return geom::Vec3(); + return nrvo; }