Triangle: Is a Point inside or outside?

Consider a triangle ABC specified by three non-collinear points A, B abd C. Any point P in the plane of the points can then be unqiquely expressed as:

P= uA+vB+wC

for some constants u, w, and w, where u+v+w = 1/

The triplet (u,v,w) correspodins to the barycentric coordinates of tthe points.

For triangle ABC, the baycentric coordinates of the cetrcies A, B and C are (1, 0, 0), (0,1,0) and (0,0,1) respectively.

In general, a point with barycentric coordinates (u,v,w) is inside (or on) the triangle if and only if \(0\le v, w,w \le 1\), or alternatively, if and only if \(0\le v \le 1\) and \(v+w \le 1\)

The abrycentric coordinates actually parametrrize the plane follows from P=uA+vB+wC realky just being a reformulation of 
P=A+v(B-A)+w(C-A), with v and w arbitraty as
P= A+v(B-A) +w(C-A)==(1-v-w)A+vB+wC
In the latter formulation, the two indepent diection vectors AB and AC form a coordinate system with origin A, allowing any point p in the plane to be paramterised in terms of v and w alone. Clearly barycenrtic coordinaets ia r edundant represntation in that the third component an be expressed in ters of the first two. 

To solve for the barycentric coordianets the expression reformulate the expression \(P= A+v(B-A)+w(C-A)\) as \(v(B-A)+w(C-A)=P-A\) and, for ease of analysis, let \(\mathbf v_0=B-A\), \(\mathbf v_1=C-A\) and \(\mathbf v_2=P-A\), so that the expression becomes:
\[v\mathbf v_0 + w\mathbf v_1=\mathbf v_2\]
Now dot product both sides with both \(\mathbf v_0\) and \(\mathbf v_1\)
\[(v\mathbf v_0 + w\mathbf v_1)\cdot \mathbf v_0=\mathbf v_2\cdot \mathbf v_0 \]
\[(v\mathbf v_0 + w\mathbf v_1)\cdot \mathbf v_1=\mathbf v_2\cdot \mathbf v_1 \]
Since the dot product is distributive over addition, the above expressions are equivalent to:
\[v(\mathbf v_0\cdot \mathbf v_0)+w(\mathbf v_1\cdot \mathbf v_0)=\mathbf v_2\cdot \mathbf v_0\]
\[v(\mathbf v_0\cdot \mathbf v_1)+w(\mathbf v_1\cdot \mathbf v_1)=\mathbf v_2\cdot \mathbf v_1\]
The above are set of linear equations and hence can be solved with the Cramer's rule.


Cramer's rule for 2x2 system


For the linear system:
\[\left\{{\begin{matrix}a_{1}x+b_{1}y&={\color {red}c_{1}}\\a_{2}x+b_{2}y&={\color {red}c_{2}}\end{matrix}}\right.\]
which in matrix format is
\[\begin{bmatrix}a_1&b_1\\a_2&b_2 \end{bmatrix} \begin{bmatrix}x\\y\end{bmatrix}= \begin{bmatrix}{\color {red}c_{1}}\\{\color {red}c_{2}}\end{bmatrix}\]
Assume \(a_1b_2 − b_1a_2\) nonzero. Then, with help of determinants, x and y can be found with Cramer's rule as
\[\begin{aligned}x&={\frac {\begin{vmatrix}{\color {red}{c_{1}}}&b_{1}\\{\color {red}{c_{2}}}&b_{2}\end{vmatrix}}{\begin{vmatrix}a_{1}&b_{1}\\a_{2}&b_{2}\end{vmatrix}}}={{\color {red}c_{1}}b_{2}-b_{1}{\color {red}c_{2}} \over a_{1}b_{2}-b_{1}a_{2}},\quad y={\frac {\begin{vmatrix}a_{1}&{\color {red}{c_{1}}}\\a_{2}&{\color {red}{c_{2}}}\end{vmatrix}}{\begin{vmatrix}a_{1}&b_{1}\\a_{2}&b_{2}\end{vmatrix}}}={a_{1}{\color {red}c_{2}}-{\color {red}c_{1}}a_{2} \over a_{1}b_{2}-b_{1}a_{2}}\end{aligned}\]


Javascript Implementation

let v0 = b - a, v1 = c - a; v2 = p - a;
let d00 = v0.dot(v0);
let d10 = v1.dot(v0);
let d01 = v0.dot(v1);
let d11 = v1.dot(v1);
let d20 = v2.dot(v0);
let d21 = v2.dot(v1);
let D = d00*d11-d10*d01;
let v = (d11*d20 - d01*d21) / D;
let w = (d00}d21 - d01*d20) / D;
let u = 1 - v - w;

As dot product is commutative, d10 and d01 are equal. I have calculated them separately for ease of comprehension.




時刻: 0 件のコメント:

Line and Line Segment: Which side of Line does a Point lie?

We want to determine which side point \(P=(x,y)\) lies with respect to the line that passes through points \(A=(x_1, y_1)\) and \(B=(x_2,y_2)\)



Solution


Calculate:
\[d=(x−x_1)(y_2−y_1)−(y−y_1)(x_2−x_1)\]
If d>0 then it lies on the "positive side" of the line (ie in the direction that n is pointing)
If d<0 then the point lies on the other side of the line,
If d=0 then the point lies exactly on the line.

Geometric interpretation



Line AB can be defined as:
 \[(x_2,y_2)−(x_1,y_1)=(x_2−x_1,y_2−y_1)\]
The orthogonal (perpendicular) direction to that line will be (flip the x and y and negate one component i.e (y, -x)):
\[\mathbf n= ((y_2−y_1),−(x_2−x_1))\]
One possible vector going from the line to the Point \(P\) is:
\[D=P−A=(x−x_1, y−y_1)\]

Dot product between n and D


The expression for d above is actually the dot product between \(\mathbf n\) and \(D\).

We know that for the dot product between 2 vectors:
  • If angle between the 2 vectors is less than 90 degrees, value of dot product is positive.
  • If angle between the 2 vectors more than 90 degrees, negative.
  • If the vectors are perpendicular to each other, dot product is zero.
Which explains the solution above

Cross or Perp-Dot product between lines AB and AP


Another interpretation is to express the problem as a perp-dot product, ie calculating the signed area of the parallelogram formed by vectors AB and AP.

Perp-dot-product between two vectors \(a\) and \(b\) is defined as:
\[\mathbf a\perp \mathbf b=\mathbf a^\perp\cdot \mathbf b=a_xb_y-a_yb_x=\begin{vmatrix} a_x &a_y \\ b_x &b_y \end{vmatrix}\]

Referring to the diagram above, we are interested in the perp-dot product between \(\vec {AB}\) and \(\vec {AP}\).
\[\begin{align} d&=\vec {AB}\perp \vec {AP}=\vec {AB}^\perp\cdot \vec {AP}\\ &= (x_2-x_1,y_2-y_1)^\perp\cdot (x-x_1,y-y_1)\\ &=(y_2-y_1, x_1-x_2)\cdot (x-x_1,y-y_1)\\ &=(y_2-y_1)(x-x_1)+(x_1-x_2)(y-y_1)\\ &=(x-x_1)(y_2-y_1)-(y-y_1)(x_2-x_1) \end{align}\]

hence same as the above.

Closest Point on Plane to Point


Referring back to the original diagram, the blue line can be interpreted as a plane, defined by the normal \(\mathbf n\) and either of the points \(A\) or \(B\), ie all points \(X\) on the plane satisfy the equation \(\mathbf n\cdot (X-A)=0\) (or \(\mathbf n\cdot (X-B)=0\)).

Let closest point R on the plane to P be the orthogonal projection of P onto the plane, obtained by moving P perpendicularly (with repsect to \(\mathbf n\)) towards the plane. That is, \(R=P-t\mathbf n\) for some value of t. Inserting this expression for R into the plane equation and solving for t gives:
\[\begin{aligned} &\mathbf n\cdot((P-t\mathbf n)-A)=0\\ &\mathbf n\cdot P-t(\mathbf n\cdot \mathbf n)-\mathbf n\cdot A=0\\ &\mathbf n\cdot(P-A)=t(\mathbf n\cdot \mathbf n)\\ &t=\mathbf n\cdot(P-A)/(\mathbf n\cdot \mathbf n)\\ &t=\mathbf n\cdot D/(\mathbf n\cdot \mathbf n) \end{aligned}\]
(Notice that the numerator of the last expression is the same as the original solution!)

Taking this expression further, substituting this expression for \(t\) in \(R=P-t\mathbf n\) gives the projection point \(R\) as simply:
\[R=P-(\mathbf n\cdot(P-A)/(\mathbf n\cdot \mathbf n) )\mathbf n\] 
Hence it is clear that \(t=\mathbf n\cdot(P-A)\) correspods to the signed distance of \(P\) from the plane in units of the length of \(\mathbf n\).

If \(t\) is positive, \(P\) is 'in front' of the plane, and if negative \(P\) is 'behind' the plane.


Useful References



時刻: 0 件のコメント:

Numerical Integration: Euler methods

Introduction


When creating physics simulations, a mathematical model that describes the motion of the bodies - equations of motion - must be created. Then those equations of motion must be solved.

Numerical integration is the branch of numerical analysis that deals with the solution of differential equations. Fundamentally, we need an interation scheme to simulate motion because it involves solving Newton's second law of motion which is a a diferential equation.

As a reminder, Newton's second law is \(F = ma\).

Acceleration is actually a derivative of velocity \(a= dv/dt\), so that we can rewrite the second law as:
\[F=m\frac{dv}{dt}\]
This is so called a first order differential equation with respect to velocity because it involves the first derivative of velocity.

By recalling that \(v= dx/dt\), and expressing force \(F\) as a function of position and velocity, we can rewrite this as:
\[F(x,t)=m\frac{d^2x}{dt^2}\]
This is now a second order differential equation with respect to displacement, as it relates the second derivative of the position to the force applied.

In principle, the precededing differential equation can be solved by integrating analytically or numerically to yield the diplacement \(x\) and velocity \(v\) as a function of time. Specifically you can integrate the first equation to arrive at velocity, then integrate the velocity to arrive at displaycement.

Analytical solution


For simple equations of motion, you can use calculus and other mathematical techniques to find the exact solution. This is called the analytic solution. It is also referred to as a closed form solution.

Numerical solution


Most physics simulations (even relatively simple ones) are too complicate to be able to find an analytic solution. Instead numerical analysis is used to find an approximate solution.

Problem Statement


Let an initial value problem be specified as follows:
\[v= \frac {dx}{dt}=f(x, t),\quad x(t_{0})=x_{0}\]
The displacement \(x\) is an unknown function of time \(t\), which we would like to approximate.

We are told that velocity \(v=\frac{dx}{dt}\), the rate at which \(x\) changes, is a function of time \(t\) and of \(x\) itself.

At the initial time \(t_0\), the corresponding \(x\) value is \(x_0\). The function \(f\) and the initial conditions \(t_0,\, y_0\) are given.

Given \(x_n\) and \(t_n\) we want to calculate an approximation for \(x_{n+1}\) at a brief time later, \(t_{n+1}=(t_n+h)\), where \(h\) is the step-size, ie the "brief time".

The above is indeed what a physics simulation needs to solve, using numerical integration.

Popular Integration Methods


There are many numerical integration methods. In this post we will look at probably the simplest of such methods called the Euler method, including a number of its variants.
  1. Explicit Euler method (often referred to simply as Euler method)
  2. Implicit Euler method (also referred to as Backward Euler)
  3. Semi-Implicit method (Symplectic Euler method)

Euler method


The simplest method is to use finite difference approximations. 

From any point on a curve, you can find an approximation of a nearby point on the curve by moving a short distance along a line tangent to the curve.

\(x′\) is replaced by the finite difference approximation\[x'(t)\approx {\frac {x(t+h)-x(t)}{h}}\]
which when re-arranged yields the following formula
\[x(t+h)\approx x(t)+hx'(t)\]
which gives:
\[x(t+h)\approx x(t)+hf(x(t), t)\]

This formula is usually applied in the following way.

We choose a step size h, and we construct the sequence \(t_{0},t_{1}=t_{0}+h,t_{2}=t_{0}+2h\),... We denote by \(x_{n}\) a numerical estimate of the exact solution \(x(t_{n})\). These estimates can be calculated by the following recursive scheme
\[x_{n+1}=x_{n}+hf(x_{n},t_{n})\]
or,
\[x_{n+1}=x_{n}+k_1\]
where \(k_1=hf(x_{n},t_{n})\)

This Euler method is an example of an explicit method. This means that the new \(x_{n+1}\) is defined in terms of things that are already known, like \(x_n\).

Sometimes our differential equation which requires solving is something like below where we start with the second derivative, like in Newton's second law (acceleration is function g of x and t):
\[F=ma=m\frac{d^2x}{dt^2}=mg(x,t)\]
For this second derivative ODE we can apply Euler's method by defining:
\[v=\frac{dx}{dt}\]
which menas our differential equation becomes
\[\frac{dv}{dt}=g(x_n,t_n)\]
which is solved by using the Euler method:
\[v_{n+1}=v_n+hg(x_n,t_n)\]
This gives us two Euler equations:
\[x_{n+1}=x_{n}+hf(x_{n},t_{n})\]
\[v_{n+1}=v_n+hg(x_n,t_n)\]
Writing this in terms of variable names more familiar to games programming gives:
\begin{aligned}x_{n+1}&=x_n+v_n\cdot \Delta t\\v_{n+1}&=v_n+a_{n}\cdot \Delta t\end{aligned}

This is the forward Euler method, in contrast with a similar but different backward Euler method described later.

Error term


\(x(t)\) can be expressed as follows by means of Taylor's expansion.
\[x(t+h)=x(t)+hx'(t)+\frac{1}{2!}h^2x''(t)+\frac{1}{3!}h^3x'''(t)+\cdots\]
You can see that the Euler method is simply the above equation with exerything after the first two terms truncated.

Euler formula would be correct if all derivatives beyond the first were zero, otherwise there is an error. The error is described as \(O(h^2)\).
\[x_{n+1}=x_n+k_1+O(h^2)\]

Energy of the system


I will not delve into the details here, but using forward Euler integration has a tendency to "add" energy, thereby violating the energy conversation principle. For physics simulations this is a problem (for a very good explanation, see here). 

Implicit Euler method


This is the backward Euler method which uses the following approximation.\[y'(t)\approx {\frac {y(t)-y(t-h)}{h}}\]
which, following the same derivation process as above gives:\[y_{n+1}=y_{n}+hf(t_{n+1},y_{n+1})\]Writing this in terms of variable names more familiar to games programming gives:\[\begin{aligned}v_{n+1}&=v_n+a_{n+1}\cdot\Delta t\\x_{n+1}&=x_n+v_{n+1}\cdot \Delta t\end{aligned}\]
This relies on knowing the future value of acceleration \(a_{n+1}\); hence "implicit". This tend not to be used for games physics simulations.

Semi-Implicit Euler method


This is probably the most popular method, used by a lot of game physics engines.
\[\begin{aligned} v_{n+1}&=v_n+a_{n}\cdot \Delta t\\  x_{n+1}&=x_n+v_{n+1}\cdot \Delta t \end{aligned}\]
時刻: 0 件のコメント:
ラベル: ,
登録: コメント (Atom)

Numerical Integration: 2nd and 4th order Runge-Kutta schemes

Second-order Runge-Kutta scheme (RK2) There are several variations of the second order Runge-Kutta schemes of numerical integration. They ar...