Geometrical calculations : points, lines, planes : intersections, distances, angles

[Gibbz]

nice stuff :)

takes me back to highschool math days....

[prettyPixel]

line and line:

Line-Line Intersection : intersection of two lines AB and CD:
fn lineLineIntersect pA pB pC pD = (
local a=pB-pA
local b=pD-pC
local c=pC-pA
local cross1 = cross a b
local cross2 = cross c b
pA + ( a*( (dot cross2 cross1)/((length cross1)^2) ) )
)
or by using vectors:
fn lineLineIntersect pA a pB b = (
local c=pB-pA
local cross1 = cross a b
local cross2 = cross c b
pA + ( a*( (dot cross2 cross1)/((length cross1)^2) ) )
)
see also : Eric W. Weisstein. From MathWorld--A Wolfram Web Resource.
http://mathworld.wolfram.com/Line-LineIntersection.html




Line-Line Distance : distance between two skew lines
fn lineLineDist pA a pB b = (
local c=pB-pA
local crossAB=cross a b
abs (dot c crossAB) / (length crossAB)
)

see also : Eric W. Weisstein. From MathWorld--A Wolfram Web Resource.
http://mathworld.wolfram.com/Line-LineDistance.html




Line-Line // Distance : distance between two parallel lines
use "Point-Line Distance" : pA-pB is the first line and pC is a point of the second line




Line-Line Angle : angle between two lines:
fn getVectorsAngle vAB vCD = (
acos (dot (normalize vAB) (normalize vCD))
)
or
fn lineLineAngle pA pB pC pD = (
local vAB=pB-pA
local vCD=pD-pC
local angle=acos (dot (normalize vAB) (normalize vCD))
if angle<90.0 then angle else (180.0-angle)
)
see also : Eric W. Weisstein. From MathWorld--A Wolfram Web Resource.
http://mathworld.wolfram.com/Line-LineAngle.html

[prettyPixel]

plane and point:

Point-Plane Projection : find the point on the plane ABC which is the projection of the point D
fn pointPlaneProj pA pB pC pD = (
local nABC=normalize (cross (pB-pA) (pC-pA))
pD+((dot (pA-pD) nABC)*nABC)
)



Point-Plane Distance : find the distance between a plane and a point
fn pointPlaneDist pA pB pC pD = (
local nABC=normalize (cross (pB-pA) (pC-pA))
length ((dot (pA-pD) nABC)*nABC)
)
see also : Eric W. Weisstein. From MathWorld--A Wolfram Web Resource.
http://mathworld.wolfram.com/Point-PlaneDistance.html

[prettyPixel]

plane and line:

Line-Plane Intersection : intersection between a line and a plane
fn planeLineIntersect planePoint planeNormal linePoint lineVector = (
local lineVector=normalize lineVector
local d1=dot (planePoint-linePoint) planeNormal
local d2=dot lineVector planeNormal
if abs(d2)<.0000000754 then ( if abs(d1)>.0000000754 then 0 else -1 )
else ( linePoint + ( (d1/d2)*lineVector ) )
)
This script is based on an original script of Joshua Newman.
return 0 for parrallel (no intersections) results
and -1 for co-incident (infinate) results
http://forums.cgsociety.org/showthread.php?t=290255
see also : Eric W. Weisstein. From MathWorld--A Wolfram Web Resource.
http://mathworld.wolfram.com/Line-PlaneIntersection.html



Line-Plane Distance : distance between a plane and a line //
use "Point-Plane Distance" : where 'point' is a point of the line



Vector perpendicular to a line into a Plane : the perpendicular of AB into the plane ABC
vectorPerp=cross vectorAB (cross vectorAB vectorAC)
-- or
vectorPerp=cross (cross vectorAB vectorAC) vectorAB
see also : Eric W. Weisstein. From MathWorld--A Wolfram Web Resource.
http://mathworld.wolfram.com/Perpendicular.html

[j-man]

plane and line:

Line-Plane Intersection : intersection between a line and a plane
fn planeLineIntersect planePoint planeNormal linePoint lineVector = (
local lineVector=normalize lineVector
local d1=dot (planePoint-linePoint) planeNormal
local d2=dot lineVector planeNormal
if abs(d2)<.0000000754 then ( if abs(d1)>.0000000754 then 0 else -1 )
else ( linePoint + ( (d1/d2)*lineVector ) )
)
This script is based on an original script of Joshua Newman.
return 0 for parrallel (no intersections) results
and -1 for co-incident (infinate) results
http://forums.cgsociety.org/showthread.php?t=290255
see also : Eric W. Weisstein. From MathWorld--A Wolfram Web Resource.
http://mathworld.wolfram.com/Line-PlaneIntersection.html



Line-Plane Distance : distance between a plane and a line //
use "Point-Plane Distance" : where 'point' is a point of the line



Vector perpendicular to a line into a Plane : the perpendicular of AB into the plane ABC
vectorPerp=cross vectorAB (cross vectorAB vectorAC)
-- or
vectorPerp=cross (cross vectorAB vectorAC) vectorAB
see also : Eric W. Weisstein. From MathWorld--A Wolfram Web Resource.
http://mathworld.wolfram.com/Perpendicular.html


Thanks for the aknowelgement PP.

Still can't get around using the tolerence then eh?

Great resource BTW, keep up the good work!

J¬

[JHaywood]

This is great, thanks.

[erilaz]

This should be a sticky! :)

[prettyPixel]

Thanks guys :)

jman: The tolerance seems to work for me, but I use rather the function in cases where the line is not parallel to the plan. Can you give an example where the tolerance is an issue?

Well, the continuation:

[prettyPixel]

plane and plane:

Plane-Plane Intersection : find the line which is the intersection of 2 planes
p1 : a point of the plane 1
n1 : the normal of the plane 1
fn planePlaneIntersect p1 n1 p2 n2 = (
-- n1, n2 are normalized
local lineVector = cross n1 n2
local proj1=(dot n1 p1)*n1
local proj2=(dot n2 p2)*n2
local perp1=cross n1 (normalize lineVector)
local perp2=cross n2 (normalize lineVector)
local cr = cross (proj2-proj1) perp2
local intersectionPoint = proj1 + (perp1*( (dot cr lineVector) / ((length lineVector)^2)) )
ray intersectionPoint lineVector
)
I think that it is not the easiest solution. If you know a better means...
see also : Eric W. Weisstein. From MathWorld--A Wolfram Web Resource.
http://mathworld.wolfram.com/Plane-PlaneIntersection.html

Plane perpendicular to a Plane and a Line : find the plane which is perpendicular to the plane n and to the line AB
It's the plane ABC where C=A+n
see also : Eric W. Weisstein. From MathWorld--A Wolfram Web Resource.
http://mathworld.wolfram.com/ParallelPlanes.html

Plane-Plane Angle : find the angle between two planes n1 and n2 (normals)
fn getVectorsAngle n1 n2 = ( acos (dot (normalize n1) (normalize n2)) )

[j-man]

Can you give an example where the tolerance is an issue?

Yes, if the results are parrallel or co-incident the rounding error tells the function from telling which!

J¬

[JHaywood]

Hey PP, just wondering if you're done adding to this thread. I want to put together a word doc containing all these functions, but wanted to wait until you were fiinished.

[JHaywood]

I've got a math question, and this seems like a good place to ask. I tried digging through Mathworld, but it's way over my head.

I'd like to be able to find a point on an arc defined by the end point of an object rotated around it's pivot, along a specific axis. For example, let's say there's a bone object with a end bone as it's child. If I rotate the parent bone around it's Y axis 10 degrees, the end bone is now at a new position. What's the formula for finding that position?

This page at Mathworld describes what I'm talking about, except in 2D: http://mathworld.wolfram.com/Arc.html

[erilaz]

I've got a math question, and this seems like a good place to ask. I tried digging through Mathworld, but it's way over my head.

I'd like to be able to find a point on an arc defined by the end point of an object rotated around it's pivot, along a specific axis. For example, let's say there's a bone object with a end bone as it's child. If I rotate the parent bone around it's Y axis 10 degrees, the end bone is now at a new position. What's the formula for finding that position?

This page at Mathworld describes what I'm talking about, except in 2D: http://mathworld.wolfram.com/Arc.html

That should be a straight trig calculation where (skeleton code):

r = $bone01.length
theta = (rotation of your bone)
phi = 90 - theta
y = sin(theta) * r
x = cos(theta) * r
z = sin(theta) * r * cos(phi)


I'm a little hazy on the Z phi stuff, so somebody correct me if i'm wrong!

[prettyPixel]

Hey PP, just wondering if you're done adding to this thread. I want to put together a word doc containing all these functions, but wanted to wait until you were fiinished.

I think that the essential formulae are there but anyone can add some more. The possibilities are illimited.
Here is an example: this script discovers if a line intersects a sphere.
It returns true or false.

fn lineSphereIntersection sphereCenter sphereRadius linePoint lineVector = (
local nLineVector= normalize lineVector
local projPoint=linePoint+((dot (sphereCenter-linePoint) nLineVector)*nLineVector)
local dist=distance sphereCenter projPoint
if dist>sphereRadius then false else true
)
I use the projection of a point on a plane which passes by the sphere center. If the distance between this point and the center is lower or equal to the sphere radius then the line intersects the sphere. It's based on the function 'Point-Plane Projection' .

Another way of going farther is to write the same functions but on objects like triangular faces or segments. For example find the position of a line that intersects a mesh face (not a plane). The case is more restricted.

[JHaywood]

Thanks Martin. That got me going in the right direction, along with this site: http://www.kirupa.com/developer/actionscript/trigonometry.htm

Here's what I came up with to illustrate what I wanted to do:

global polylinePosArr = #()

fn drawCircle = (
gw.setColor #line green
gw.setTransform(Matrix3 1)
gw.Polyline polylinePosArr true
gw.enlargeUpdateRect #whole
gw.updateScreen()
)

fn drawObjectCirle obj radius step:36 axis:1 = (
polylinePosArr = #()
-- get the object's position
local p = obj.pos
for i = 0 to (360 - step) by step do (
-- make a point helper
-- define the current angle around the circle
local angle = i
-- calculate the X and Y coordinates of the point on the circle using the current angle
local x = cos(angle) * radius
local y = sin(angle) * radius
-- 1 = X axis
-- 2 = Y axis
-- 3 = Z axis
global p3 = case axis of (
1:[0,x,y]
2:[x,0,y]
3:[x,y,0]
)
-- using the coordinate system of the object, define the coordinates of the current point
newP = ( (matrix3 [1,0,0] [0,1,0] [0,0,1] p3) * (obj.transform) ).translationPart
append polylinePosArr newP
)
)

-- create a box
b = box length:1 width:1 height:10 pos:(random [-10,-10,-10] [10,10,10])
b.rotation = quat (random -1 1) (random -1 1) (random -1 1) 1
-- get the length of the box, which will be used as the radius for the circle calculations
l = b.height

-- be careful with the "steps" parameter, too low (around 10) and Max will crash
drawObjectCirle b l step:18 axis:2

-- draw the circle
registerRedrawViewsCallback drawCircle
redrawViews()

/* stop drawing the circle
unregisterRedrawViewsCallback drawCircle
*/

[prettyPixel]

The quaternion is easy to use to rotate a position in 3D space. You know:
1- the axis of rotation (it must be normalized)
2- the center of the rotation
3- the angle
4- the position of the point to be moved

rotAxis = [0,0,2]
rotCenter = [0,2,0]
rotAngle = 90
thePoint = [1,2,3]

q = quat rotAngle (normalize rotAxis)
pointAfterRotation = (((thePoint - rotCenter ) * q) + rotCenter)
--> [0,1,3]

The direction of the rotation can be modified by changing the direction of the axis.
if you use rotAxis = -[0,0,2], the rotation will be reversed.

Here is a variation of your function:

fn drawObjectCirle obj radius step:36 axis:1 = (
polylinePosArr = #()
local rotCenter = obj.pos
local rotAxis = obj.dir
local p
case axis of (
1: p=[radius,0,0]*obj.transform
2: p=[0,radius,0]*obj.transform
3: p=[0,0,radius]*obj.transform
)
q = quat step rotAxis
for i = 0 to (360 - step) by step do (
p = ((p-rotCenter)*q) + rotCenter
append polylinePosArr p
)
)

[proteus2002]

hi all

i am trying to make a sort of "one axis lookat" scriptcontroller.
object_a is flying around in 3d space and object_b is located at position [0,0,0]
object_b has a script_contoller at the z-rotation so that it looks at object_a.

this is what i have so far as scriptcontroller for object_b z-rotation:
-----------------------------------------

dependson $object_a
v1t = at time currenttime $object_a.position
a = normalize [v1t.x,v1t.y,0]-- set z to 0 because it is not needed
upv = [0,1,0] --defining a 2d "up vector"

cosine = a.x * upv.x + a.y * upv.y
if cosine > 1 then cosine1 = 1
if cosine < -1 then cosine1 = -1
if (a.x * upv.y - a.y * upv.x) < 0 then
(
result = -cosine
)
else
(
result = cosine
)
result


-------------------------------------------------------------------

it does not work
any ideas???

thanks

[proteus2002]

this seems to work:



dependson $object_a
v1t = at time currenttime $object_a.position
a = normalize [v1t.x,v1t.y,0]
ndir = normalize a

upvector = [0,0,1]
rightvector = normalize (cross upvector ndir )
thematrix = matrix3 ndir rightvector upvector ndir

pitch = ((thematrix.rotation as eulerangles).z +180)*pi/180


But is there a way not to use the matrix3 calculations?
any ideas???

[hblan]

that is great topic . espacially the math website , thanks .

[ajohnson]

I concur... this is really great, props to prettyPixel and everyone else that contributed!

[sleepite]

...But is there a way not to use the matrix3 calculations?
any ideas???

The code I use for doing a look-at type operation is:
fn PointAxisAt a which_axis b =
(
local c = normalize (b - a)
local angle = acos (dot which_axis c)
local axis = normalize (cross which_axis c)

local out = (angleaxis angle axis) as quat
out
)

a = $Pyramid01.position
b = $Point01.position
q = PointAxisAt a z_axis b
r = q * (inverse $Pyramid01.rotation)
rotate $Pyramid01 r

[dutch_delight]

THought this might be the best topic to post it:

I have this camera data:
posz=-12.896441
posy=10.177296
posx=-248.952271

And this for the angle:
anglez=1.000000
angley=0.000000
anglex=0.000000


this is how they calculated the angle:
angle = camera_target - camera
angle = normalize angle
anglex = angle.x
angley = angle.y
anglez = angle.z


how would i go about transforming the angle from this data?
how can i create a target or free camera from this?

[arketip]

I prefer to rename "angle" in "direction" because the angle can remind the angle of the camera (fov)

thePos = [posx,posy,posz]
theDirection = [directionx,directiony,directionz]
theLength = 100.0
theTargetPos = thePos + ((normalize theDirection)*theLength)

myCam = targetCamera pos:thePos target:(targetObject pos:theTargetPos)


Regrettably your data are not complete : the camera's FOV is missing.

[dutch_delight]

Thanks for your reply but i already solved it. I did it slightly different though:

camera_target_position = point3 (posx + directionx*10) (posy + directiony*10) (posz + directionz*10)

Which returns the correct 'direction'.


Oh and you dont need a FOV to create a camera, it uses default setting when you leave it out.

[xissburg]

I just wanna give you two link for resources about geometry (intersections , distance computations etc):

http://www.softsurfer.com (algorithms section)
http://www.geometrictools.com

Hope this helps

[davestewart]

Great thread! Seeing basic concepts as code snippets is fantastic!

I do think that the posting of questions (as opposed to answers) clutters the "repository" feel, though. For those asking (new and unrelated) questions wouldn't it be more helpful to start and finish a new thread, THEN post (a concise answer) here?

:D

[erilaz]

I do think that the posting of questions (as opposed to answers) clutters the "repository" feel, though. For those asking (new and unrelated) questions wouldn't it be more helpful to start and finish a new thread, THEN post (a concise answer) here?

:D

This is being transferred to the CGWiki bit by bit, so the information can be kept clean:
http://wiki.cgsociety.org/index.php/Geometric_Calculations_%283dsmax%29

[yuenlw]

Hy all,

I am new at maxscript and I want to calculate the shortest distance between a vertex in a skin and a bone. Can someone help me?


-edit-

ok, i modified pretty pixel's code and got this:



fn pointSegmentDist pA pB pC =
(
local vAB=pB-pA
local vAC=pC-pA
local vBC=pC-PB

local c1=dot vAB vAC
local c2=dot vAB vAB

if c1 < 0 do
return length vAC

if c2 < c1 do
return length vBC

(length (cross vAB vAC))/(length vAB)
)


not sure if it is 100% correct, i have other issues but will start a new thread for that.

[rdg]

The Algorithmic Beauty of Plants can be downloaded here as PDF:
http://algorithmicbotany.org/papers/

This is the 'definitve book' about L-Systems (see also Challenge #16 (http://forums.cgsociety.org/showthread.php?f=98&t=479116)).
Sometimes you can find a printed version in 'antique shops' - get it!
I ordered mine some weeks ago, hope it arrives soon ...

Georg

[rdg]

this is a maxscript version of Robert Penners Easing Equotations for Flash:

struct rpEFn (
-- ---
-- Easing Equotations by Robert Penner
-- http://www.robertpenner.com/
-- The BSD License: http://www.opensource.org/licenses/bsd-license.php
-- ---
-- maxScript port by Georg Duemlein
-- 2007-04-18
-- http://www.preset.de
-- ---
-- previews: http://www.gizma.com/easing/
-- ---
-- t: current time
-- b: start value
-- c: end value
-- d: duration
-- ---
-- simple linear tweening - no easing, no acceleration
fn linearTween t b c d = (
c*t/d + b
),
-- quadratic easing in - accelerating from zero velocity
fn easeInQuad t b c d = (
t /= d
c*t*t + b
),
-- quadratic easing out - decelerating to zero velocity
fn easeOutQuad t b c d = (
t /= d
-c * t * (t-2) + b
),
-- quadratic easing in/out - acceleration until halfway, then deceleration
fn easeInOutQuad t b c d = (
t /= d/2
if (t < 1) then (
c / 2 * t * t + b
)else(
t-=1
-c / 2 * (t* (t - 2) - 1) + b
)
),
-- cubic easing in - accelerating from zero velocity
fn easeInCubic t b c d = (
t /= d
c*t*t*t + b
),
-- cubic easing out - decelerating to zero velocity
fn easeOutCubic t b c d = (
t /= d
t-=1
c*(t*t*t + 1) + b
),
-- cubic easing in/out - acceleration until halfway, then deceleration
fn easeInOutCubic t b c d = (
t /= d/2
if (t < 1) then (
c/2*t*t*t + b
)else(
t -= 2
c/2*(t*t*t + 2) + b
)
),
-- quartic easing in - accelerating from zero velocity
fn easeInQuart t b c d = (
t /= d
c*t*t*t*t + b
),
-- quartic easing out - decelerating to zero velocity
fn easeOutQuart t b c d = (
t /= d
t -= 1
-c * (t*t*t*t - 1) + b
),
-- quartic easing in/out - acceleration until halfway, then deceleration
fn easeInOutQuart t b c d = (
t /= d/2
if (t < 1) then(
c/2*t*t*t*t + b
)else(
t -= 2
-c/2 * (t*t*t*t - 2) + b
)
),
-- quintic easing in - accelerating from zero velocity
fn easeInQuint t b c d = (
t /= d
c*t*t*t*t*t + b
),
-- quintic easing out - decelerating to zero velocity
fn easeOutQuint t b c d = (
t /= d
t -= 1
c*(t*t*t*t*t + 1) + b
),
-- quintic easing in/out - acceleration until halfway, then deceleration
fn easeInOutQuint t b c d = (
t /= d/2
if (t < 1) then (
c/2*t*t*t*t*t + b
)else(
t -= 2
c/2*(t*t*t*t*t + 2) + b
)
),
-- sinusoidal easing in - accelerating from zero velocity
fn easeInSine t b c d = (
-c * cos (t/d * (180 / 2)) + c + b
),
-- sinusoidal easing out - decelerating to zero velocity
fn easeOutSine t b c d = (
c * sin (t/d * 90) + b
),
-- sinusoidal easing in/out - accelerating until halfway, then decelerating
fn easeInOutSine t b c d = (
-c/2 * (cos(180 * t/d) - 1) + b
),
-- exponential easing in - accelerating from zero velocity
fn easeInExpo t b c d = (
c * pow 2 (10 * (t/d - 1)) + b
),
-- exponential easing out - decelerating to zero velocity
fn easeOutExpo t b c d = (
c * (-pow 2 (-10 * t/d) + 1 ) + b
),
-- exponential easing in/out - accelerating until halfway, then decelerating
fn easeInOutExpo t b c d = (
t /= d/2
if (t < 1) then (
c / 2 * pow 2 (10. * (t - 1)) + b
)else(
t -= 1
c / 2 * (-pow 2 (-10. * t) + 2 ) + b
)
),
-- circular easing in - accelerating from zero velocity
fn easeInCirc t b c d = (
t /= d
-c * (sqrt(1 - t*t) - 1) + b
),
-- circular easing out - decelerating to zero velocity
fn easeOutCirc t b c d = (
t /= d
t -= 1
c * sqrt(1 - t*t) + b
),
-- circular easing in/out - acceleration until halfway, then deceleration
fn easeInOutCirc t b c d = (
t /= d/2
if (t < 1) then (
-c/2. * (sqrt(1. - t*t) - 1.) + b
)else(
t -= 2
c/2. * (sqrt(1. - t*t) + 1.) + b
)
)
)

-- sample
d = 50.
t = 0.

while t < d do (
v1 = v2 = v3 = 0
v1 = rpEFn.linearTween t 0. 100. d
v2 = rpEFn.easeInOutQuad t 1. 100. d
v3 = rpEFn.easeInOutQuad t 0. 1. (d/2)
sp = sphere pos:[v1, 0, v2] radius:(.5 + v3)
sp.wirecolor = [255,255,255] * v3
-- ---
t += .5
max zoomext sel all
)


maybe this is usefull for somebody

http://www.robertpenner.com/
http://www.gizma.com/easing/

2007-05-08
I found this API documentation:
http://livedocs.adobe.com/flash/9.0/ActionScriptLangRefV3/fl/transitions/easing/Elastic.html
and c is the "total change in the animation property".

I also recomend the " Tweening chapter of my book (http://www.robertpenner.com/easing/penner_chapter7_tweening.pdf)" at Robert Penners site.
I just started reading it and it explains a lot.

Georg

[pix3l]

Hey people, lovely math here :) heheh i'v forgotten most of it heheh.
I'm trying now to get the angle between 2 points.

So i looked in the reference. there it says the following :
to calculate the angle use the following formula :
theAngle = acos(dot (normalize v1) (normalize v2))

So i tried this in max. if i select 0.0.0 and 150.661,150.661,0 i get 45. till now all okey.
But if i select 10,10,0 and 150.661,150.661,0 i get 0.0197823

Hmmmmm why?
for both two the angle should be 45. but i get different results.
Any1 who can shed some light on this ?

Many thx in advance!

[xissburg]

Very strange to get 0.0197823, you should get 0 since the vectors are pointing in the same direction(the angle between them is 0 :) ). Maybe its a float point accuracy problem at the normalization(1 inversion, 3 multiplies and one square root) and the acos. But anyway, the formula is correct.

[magicm]

So i tried this in max. if i select 0.0.0 and 150.661,150.661,0 i get 45. till now all okey.
But if i select 10,10,0 and 150.661,150.661,0 i get 0.0197823

Hmmmmm why?
The normalize function returns a unit vector, meaning that it's length equals 1 unit. Normalizing [0,0,0] will return [1,0,0], so this explains why the result of your first example is 45. Your second example should return 0, but returns 0.0197823 instead because of round-off errors.

1. [0,0,0] [150.661,150.661,0]
== [1,0,0] [0.707107,0.707107,0] (normalized)
== 45

2. [10,10,0] [150.661,150.661,0]
== [0.707107,0.707107,0] [0.707107,0.707107,0] (normalized)
== 0
Hope this helps,
Martijn

[Guimas]

Nice thread!

Any idea in how to test if a point is inside a n-dimensional polyhedron just by coordinates? If one has all the vertexes coordinates of a 3d solid and a point3, how to check if the point3 is inside?

cheers

[prettyPixel]

I don't have the complete solution but some ideas:

First try a fast test : find if the point is completely outside the volume by using the distance.
For example: you find the center of the bounding box of the volume and the point which is farthest. If your point is farther then it is not inside the volume. It's obvious but that avoid to apply a complex function each time...

After... I'm not sure that is 100% secure but if you use the normal of the faces I think we can have a good approach to do that.
For a convex volume, a point is inside a volume if this point is behind of each face of the volume.
For each face:
+ find the vector of the point-plane projection of the face
+ calculate if the vector have the same direction than the normal of the face
+ if all vectors have the same direction then your point is inside the convex volume

The problem appears when you have a non-convex volume.
For example imagine a form like a "U" (in 3D)
If your point is at the center of the U, ok that does work because the point is really outside the volume in itself.
... but if the point is inside of the volume then there are several faces where the normal is in opposite direction of the vector of the point-plane projection.

Probably the solution is to divide the concave volume in several convex volumes. After you test for each convex volume if the point is inside or not and you stop as soon as you find a point inside...

[prettyPixel]

I had another idea.
If the point is inside the volume, you can trace a ray to any direction and you intersect a face.
You compare the 2 vectors (your ray and the face normal) and if the sign of the dot product is positive then your point is inside the volume. And that could work for a concave volume :)

[nenalata]

hello,
nice one with the geometry, and thanks for starting such a thread,does anyone know how to find the intersection of a plane(or a cube, or polyhedra) and a spline?

[arketip]

Hello,

I had sent a question rather similar about splines but it seems that it is not possible to find a mathematical solution for a spline and the only solution is to convert your spline into on a certain number of lines or walking the spline with iterations.

Here is the thread:
http://forums.cgsociety.org/showthread.php?f=98&t=470089&page=5&pp=15

But if somebody can confirm this...

[marktsang]

hi arketip,
i also looked for the same thing for a long time -
as you say the only thing i found was subdividing by n interations and getting closest to one of these

i dont think you can do any other way

[arketip]

I have another question.
In this thread there is a solution for finding the intersection between a plane and a line. But do you know what is the best calculation for triangle intersection ?
I can calculate the intersection of the triangle plane and after calculate if this point is inside the triangle. I would like to find an optimized solution. Maybe with the dot product ?

[nenalata]

try this , arketip, if you can understand it, which i dont entirely



http://local.wasp.uwa.edu.au/~pbourke/geometry/linefacet/

[Anubis]

Is this correct for the intersection of three planes? I have found some situations where it does not seem to return the proper result (not when parallel)

fn threePlaneInter p1 p2 p3 selectPlanes: =
(
local d1 = (dot p1.dir p1.pos )
local d2 = (dot p2.dir p2.pos )
local d3 = (dot p3.dir p3.pos )

local num = ((d1 * (cross p2.dir p3.dir)) + (d2 * (cross p3.dir p1.dir)) + (d3 * (cross p1.dir p2.dir)))
local denom = (dot p1.dir (cross p2.dir p3.dir))

print ("num " + num as string)
print ("denom " + denom as string)
if (abs denom) < 0.0001 then (return undefined)

if selectPlanes == true then
(
selectMe = #(p1,p2,p3)
select selectMe
)
return (num / denom)
)

[Kameleon]

Here's a function for finding the circumcenter of any given triangle.


fn f_circumcenter p1 p2 p3 =
(
dp1=[0,0,0]
dp1.x=(p2.x-p1.x)
dp1.y=(p2.y-p1.y)
dp1.z=(p2.z-p1.z)
dp2=[0,0,0]
dp2.x=(p3.x-p1.x)
dp2.y=(p3.y-p1.y)
dp2.z=(p3.z-p1.z)
tnx = dp1.y * dp2.z - dp2.y * dp1.z
tny = dp1.z * dp2.x - dp2.z * dp1.x
tnz = dp1.x * dp2.y - dp2.x * dp1.y
sp1=[0,0,0]
sp1.x=(p1.x+p3.x)/2.0
sp1.y=(p1.y+p3.y)/2.0
sp1.z=(p1.z+p3.z)/2.0
dp1.x=(sp1.x-p1.x)
dp1.y=(sp1.y-p1.y)
dp1.z=(sp1.z-p1.z)
dp2.x=tnx
dp2.y=tny
dp2.z=tnz
nx = dp1.y * dp2.z - dp2.y * dp1.z
ny = dp1.z * dp2.x - dp2.z * dp1.x
nz = dp1.x * dp2.y - dp2.x * dp1.y
sp2=[0,0,0]
sp2.x=sp1.x+nx
sp2.y=sp1.y+ny
sp2.z=sp1.z+nz
sp3=[0,0,0]
sp3.x=(p1.x+p2.x)/2.0
sp3.y=(p1.y+p2.y)/2.0
sp3.z=(p1.z+p2.z)/2.0
dp1.x=(sp3.x-p1.x)
dp1.y=(sp3.y-p1.y)
dp1.z=(sp3.z-p1.z)
dp2.x=tnx
dp2.y=tny
dp2.z=tnz
nx = dp1.y * dp2.z - dp2.y * dp1.z
ny = dp1.z * dp2.x - dp2.z * dp1.x
nz = dp1.x * dp2.y - dp2.x * dp1.y
sx4=sp3.x+nx
sy4=sp3.y+ny
sz4=sp3.z+nz

ax=sp2.x-sp1.x
ay=sp2.y-sp1.y
az=sp2.z-sp1.z
bx=sx4-sp3.x
byy=sy4-sp3.y
bz=sz4-sp3.z
cx=sp3.x-sp1.x
cy=sp3.y-sp1.y
cz=sp3.z-sp1.z
qp1=[0,0,0]
qp1.x = cy * bz - byy * cz
qp1.y = cz * bx - bz * cx
qp1.z = cx * byy - bx * cy
qp2=[0,0,0]
qp2.x = ay * bz - byy * az
qp2.y = az * bx - bz * ax
qp2.z = ax * byy - bx * ay
dotp=qp1.x*qp2.x+qp1.y*qp2.y+qp1.z*qp2.z
lee=qp2.x*qp2.x+qp2.y*qp2.y+qp2.z*qp2.z
lee=sqrt lee
si=dotp/(lee * lee)
circumcenter=[0,0,0]
circumcenter.x=sp1.x+ax*si
circumcenter.y=sp1.y+ay*si
circumcenter.z=sp1.z+az*si

return circumcenter
)

[arketip]

I think that it is not the easiest solution. If you know a better means...

Here is an improvment of the function for calculating the intersection of 2 planes:

pA, nA : first plane, where p is a point and n is the normal of the plane
pB, nB : the second plane

fn PlanePlaneIntersection pA nA pB nB =
(
dir= cross nA nB
perp= cross nA dir
p= pA + (dot (pB-pA) nB) * perp / (dot perp nB)
ray p (normalize dir)
)

This is not checked in the code but if dir=[0,0,0] then the Planes are parallel...

[arketip]

Here is a function for calculating the intersection of 3 planes.

At present, that is just a Plane-Plane intersection followed by a line-Plane Intersection but I think It should exist a shorter way to calculate that.

fn PPPIntersection pA nA pB nB pC nC =
(
vD= cross nA nB
perp= cross nA vD
pD= pA + (dot (pB-pA) nB) * perp / (dot perp nB)
pD + (dot (pC-pD) nC)*vD / (dot vD nC)
)

@Anubis: I am curious to know if this code returns always the expected result (except of course if the planes are parallel.)

[Gravey]

i know that Kameleon already posted a function for the finding the circumcenter of a triangle on the previous page, but to be honest I have no idea how it works, so I have written a different solution using barycentric coordinates:fn circumcenter p1 p2 p3 =
(
BC = distance p2 p3
CA = distance p3 p1
AB = distance p1 p2

baryCoords = [ (BC^2*(CA^2+AB^2-BC^2)), (CA^2*(AB^2+BC^2-CA^2)), (AB^2*(BC^2+CA^2-AB^2)) ]
triArea = baryCoords.x + baryCoords.y + baryCoords.z
baryCoords /= triArea -- normalize the barycentric coordinates
-- substitute in for P = uA + vB + wC
baryCoords.x * p1 + baryCoords.y * p2 + baryCoords.z * p3
)
References:
Circumcenter @ Wikipedia (http://en.wikipedia.org/wiki/Circumcenter) for the barycentric circumcenter forumla
6th post on this page (http://groups.google.com/group/comp.graphics.algorithms/browse_thread/thread/aa9ce4fff33ee593?pli=1) for the formula to convert from barycentric to cartesian coordinates (P=uA+vB+wC)

[SyncViewS]

Hi guys,
I put together an essential library of geometrical calculations functions. They're approximately those already shown here a little polished, so this thread was back credited too. You can find them here (http://www.illusioncatalyst.com/mxs.php).

- Enrico

[alteredgene]

thank you everyone! great collection, found it super useful!

also thank you Syncviews for providing us with such informative site : )
will be referring to your site frequently!

[duke]

What about getting an intersection between 2 curves? Is it a matter of checking each step in the segment (as lines)?

Additionally, is there a way to do a fast check on a "spaghetti" of splines to narrow down possible intersections? I looked into some variation of delaunay (at least the radial part) but it's not really suited.

[denisT]

What about getting an intersection between 2 curves? Is it a matter of checking each step in the segment (as lines)?

Additionally, is there a way to do a fast check on a "spaghetti" of splines to narrow down possible intersections? I looked into some variation of delaunay (at least the radial part) but it's not really suited.

what curves are you asking about? Bezier? are you looking for a solution for 2D intersection? before my answering the question... why do you need it?

[Kameleon]

i know that Kameleon already posted a function for the finding the circumcenter of a triangle on the previous page, but to be honest I have no idea how it works, so I have written a different solution using barycentric coordinates:fn circumcenter p1 p2 p3 =
(
BC = distance p2 p3
CA = distance p3 p1
AB = distance p1 p2

baryCoords = [ (BC^2*(CA^2+AB^2-BC^2)), (CA^2*(AB^2+BC^2-CA^2)), (AB^2*(BC^2+CA^2-AB^2)) ]
triArea = baryCoords.x + baryCoords.y + baryCoords.z
baryCoords /= triArea -- normalize the barycentric coordinates
-- substitute in for P = uA + vB + wC
baryCoords.x * p1 + baryCoords.y * p2 + baryCoords.z * p3
)
References:
Circumcenter @ Wikipedia (http://en.wikipedia.org/wiki/Circumcenter) for the barycentric circumcenter forumla
6th post on this page (http://groups.google.com/group/comp.graphics.algorithms/browse_thread/thread/aa9ce4fff33ee593?pli=1) for the formula to convert from barycentric to cartesian coordinates (P=uA+vB+wC)

You just made me laugh my heart out. Your function just shows the difference between a guy that knows shit about math (me) and one who does (you) lol

Cheers!

[claude666]

came across this (http://www.realtimerendering.com/intersections.html#IV357) thought it could be useful in this thread. sorry if it's been linked to before