||
计算分子表面的时候,需要计算直线和球面的交点
https://en.wikipedia.org/wiki/Line%E2%80%93sphere_intersection
直线和球面的交点:
///begin with a segment from r0 to r1
///
///The points p on the line are described as p = r0 + u*(r1-r0), or, more specifically,
///
///p[0] = r0[0] + u*(r1[0]-r0[0])
///p[1] = r0[1] + u*(r1[1]-r0[1])
///p[2] = r0[2] + u*(r1[2]-r0[2])
///
///The sphere centered at c with radius r is described as
///
/// (x-c[0])^2 + (y-c[1])^2 + (z - c[2])^2 = r^2
///
///If x,y,z is the point of intersection with the line, then x=p[0], y=p[1], z=p[2], and
///we can substitute the expressions above into the sphere equation. If we evaluate the
///squared quantities ( (x-c[0])^2, etc ) and separate out the terms into terms with
///u^2, terms with u, and constant terms, we get a quadratic equation of the form
///
/// a*(u^2) + b*u + c = 0
///
///where
///
/// a = (r1[0]-r0[0])^2 + (r1[1]-r0[1])^2 + (r1[2]-r0[2])^2
/// b = 2*( ((r1[0]-r0[0])*(r[0]-c[0])) + ((r1[1]-r0[1])*(r[1]-c[1])) + ((r1[2]-r0[2])*(r[2]-c[2])) )
/// c = c[0]^2 + c[1]^2 + c[2]^2 + r0[0]^2 + r0[1]^2 + r0[2]^2 + 2( (c[0]*r0[0]) + (c[1]*r0[1]) + (c[2]*r0[2]) ) - r^2
///
/// Since this is quadratic, its solutions are
///
/// -b +- sqrt( b^2 - 4*a*c )
/// u = -----------------------------
/// 2a
///
/// Which has a descriminant (b^2 - 4*a*c).
/// CASE 1: If the discriminant is less than zero, the segment does not intersect the sphere
/// CASE 2: If the discriminant is equal to zero, then the segment intersects the sphere at u = (-b)/(2*a)
/// CASE 3: If the discriminant is greater than zero, then the line intersects the sphere at two places.
///begin by computing the components of the quadratic:
--------------------------------------
Vector3 intersect_Seg_Sphere(Vector3 pt1,Vector3 pt2,Vector3 origin, float r)
{
Vector3 pt;
float dist1=pt1.distance(origin);
float dist2=pt2.distance(origin);
if( ((dist1 >r)&& (dist2 >r)) || ((dist1 <r) &&(dist2<r) ))
{
return pt;
}
//https://en.wikipedia.org/wiki/Line%E2%80%93sphere_intersection
float quad_a,quad_b,quad_c;
quad_a=(pt2.distance(pt1))*(pt2.distance(pt1));
quad_b=2*(
((pt2.xvalue()-pt1.xvalue())*(pt1.xvalue()-origin.xvalue()))+
((pt2.yvalue()-pt1.yvalue())*((pt1.yvalue()-origin.yvalue())))+
((pt2.zvalue()-pt1.zvalue())*(pt1.zvalue()-origin.zvalue())));
float disc=(quad_b*quad_b)-4*quad_a*quad_c;
if(disc <-.00001){return pt;};
if((disc<0.00001) &&(disc > -.00001))
{
float u=(-quad_b)/(2*quad_a);
pt.setx(pt1.xvalue()+u*(pt2.xvalue()-pt1.xvalue()));
pt.sety(pt1.yvalue()+u*(pt2.yvalue()-pt1.yvalue()));
pt.setz(pt1.zvalue()+u*(pt2.zvalue()-pt1.zvalue()));
return pt;
}
if(disc > 0.00001)
{
float u1 = (-quad_b - sqrt( (quad_b*quad_b) - (4*quad_a*quad_c) ) )/(2*quad_a);
float u2 = (-quad_b + sqrt( (quad_b*quad_b) - (4*quad_a*quad_c) ) )/(2*quad_a);
//intersect point
//pt=pt1+u(pt2-pt1)
//intersect point should local between the pt1 and pt2
double u;
if( u1<0 || u1>1 ){ u = u2; }
if( u2<0 || u2>1 ){ u = u1; }
pt.setx(pt1.xvalue()+u*(pt2.xvalue()-pt1.xvalue()));
pt.sety(pt1.yvalue()+u*(pt2.yvalue()-pt1.yvalue()));
pt.setz(pt1.zvalue()+u*(pt2.zvalue()-pt1.zvalue()));
return pt;
}
}
--------------------------------------
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-10-20 01:54
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社