Discussion on Methods for Calculating the Inner Points of a Quadrilateral

chuanzhen

Hi,
In the plane, there is a quadrilateral 'abcd' and a point P. The line connecting the equal percentage segments of line ab and line dc can fill the area of the entire polygon. (The 0.25, 0.5, and 0.75 percentage segments connecting line ab and line dc, respectively, are shown in the figure below) If we fill infinitely, then there must be a line segment passing through point P. How to know what percentage of line segment passes through point P.

i_mazlov

Hi @chuanzhen ,

Please note that according to our Support Procedures your question lies out of the scope of our support as it is solely an algorithmic question.

Regarding your question. There're so many open questions in your problem statement that directly affect the path you're going to follow. Hence, I will first try to ground at least some of them and then we'll see how it can be extended to your more complex scenario if there's a need for that. Feel free to elaborate on your question on a higher scale, if my thoughts go in a wrong direction.

First, let's assume your quadrilateral is actually a parallelogram (i.e. 1. lives in 2D-space; 2. opposite sides are parallel). That's quite a strong assumption, but we'll get back to it later. Then instead of having a line that intersects point P, we can think of it as a point P_ab on the line AB (because from that your can easily reconstruct your line as it will be parallel to the side lines).

Simplified case with parallelogramm:

Let's then instead working in a 2D space (i.e. working with point P_ab), simplify it further and use the corresponding single value 0.0 <= t_ab <= 1.0, considering that your line AB is defined as f(t) = A + t * (B - A), where t = [0...1], A and B were originally in R^2 (2-component vectors), but since we only work on the line AB, we can consider them being in R (and then will convert them back to R^2 using original points A and B in 2D space).

If you look at it now, we've reduced the task to the 1D space, where we have interval from 0.0 to 1.0 and some point t_ab in between that corresponds to your initial point P. Let's also for now consider interval [0...100] instead of [0...1], and let's for now consider only integral numbers on this interval, just for the sake of simplicity.

If we now get back to your original question, with this 1D space problem scope it was reduced to a simple task of searching for the greatest common divisor (GCD) between the two numbers: t_ab and 100. By the way, I assume here that you're looking for the greatest number that would satisfy your problem statement, because number 1 would always be a correct answer as well as some other numbers that are integral divisors of t_ab and 100). Additionally we can quickly make one small optimization step according to the euclidean algorithm and look for t_x = gcd(t_ab, 100 - t_ab) instead.

Generally speaking that's already the answer you were looking for, but here the interesting part begins. The biggest issue that you have here is the space your numbers are in. It's working just fine in Z (integers), but how do you handle R (real numbers)? Depending on your desired precision the task can be quite expensive for processing. I didn't do any reasearch about it, but the naive way would be to simply make a lossy conversion from R to Z and then back to R (although there might exist some interesting tricks for that).

Let's get back to our 0.0 <= t_ab <= 1.0 now. In this case, the easy way would be to scale t_ab to some big value and work with it in Z, e.g.

t_x = gcd(t_ab * 10000, (1 - t_ab) * 10000) / 10000.0

Since you now have t_ab, you can find back P_ab that lies on segment AB. For that you just need to apply your AB formula:

{P_ab}_k = A + k * t_x * (B - A)

, where k runs from 0 to N = 1.0 / t_x, and {P_ab}_k here is a set of points that lie on line AB. The only remaining step is to get your desired lines that are parallel to side lines AD and BC and go through the set {P_ab}. That was the simplified problem.

What'd happen now if we have an arbitrary quadrilateral? Actually, the core calculation would stay the same, the only issue would be to find t_ab, because the side lines are not parallel anymore. In this case, I think the easiest way would be to find a perspective transform M from your arbitrary quadrilateral ABCD to something that is easier to work with, e.g. a unit square A'B'C'D'. Once found, you just apply this transform to your original point p: p' = M * p and continue the same route as described with the parallelogram, because square also has parallel side lines.

Arbitrary quadrilateral:

Perspective transform:

Getting another step back to your original question and extending this approach further to 3D shouldn't be a problem, as you just need to transform your arbitrarily placed quadrilateral such that it is aligned with your coordinate system basis before doing any of the above mentioned calculations. For that you'd need to calculate cross product on 2 adjacent edges of your quadrilateral, make this basis orthonormal and transform coordinates to this new basis.

In case you also like to cover scenario with non-planar quadrilateral you would likely need to use some imagination for that. At least then the 'perspective transform' trick wouldn't work anymore and you would need to go into all the depths of solving the optimization problem for the minimal distance to your point P.

Please note that there's a convenient GeometryUtilsInterface that contains lots of useful implementations. For example,
Barycentric Coordinates functions will be helpful for the 'perspective transform' part
Projection functions will help you jumping between R^3 and R^2

Cheers,
Ilia

chuanzhen

@i_mazlov Thank you for your detailed reply. It has been very helpful to me. Due to being busy with other matters, I was unable to reply in a timely manner.
Before asking a question, the method I used was to continuously segment the quadrilateral and then find the block where point P is located. I will explain how it works in the image (it is approximate, not numerical)：

The perspective transformation you mentioned inspired me. If I could directly transform an arbitrary quadrilateral into a standard square a'b'c'd’. Then the coordinates of point P after transformation would be the desired result value. However, I encountered a problem during the transformation. I could only determine which triangle region of abcd point P is in (abcd contains two triangle regions, abc and acd), and then perform the transformation to map it to the standard square a'b'c'd‘.

The following video is a comparison of these two calculation methods. (Gray represents the original quadrilateral, green+red represents the grid with higher density, and the results are converted into weights for easy observation. The green grid uses continuous segmentation of the quadrilateral, while the red grid uses a method to determine which triangle is inside and then convert it to a standard quadrilateral.)

In the video, it is felt that the method of continuously segmenting quadrilaterals performs better during operation.

Is there a more detailed explanation about the perspective transformation? I am not sure how P is mapped to **P’**in the figure

Thanks for any help!

i_mazlov

Hi @chuanzhen ,

Sorry for the delayed answer.

I've just realized that perspective transform conversion that I've suggested to you earlier isn't going to work here, as it is not an affine transform, hence it doesn't preserve ratios on edges (and preservation of collinearity is not enough in your case).

However, once I had a fresh look at the problem here, I've also realized that even if some transform M would have worked here, it would have been such an overkill for this problem. What you're actually looking for is just a UV coordinates of your point P' that lies inside of an arbitrary quadrilateral, so there's even no need to transform anything here.

When you have your ABCD quadrilateral defined in R^3 space and some point Q' = <u, v>, defined in R^2: 2D plane of this quadrilateral, then you can find your point P' (defined in R^3) by doing a bilinear interpolation along edges of your quadrilateral. However, in your case you need to solve the inverse problem, when you have your P' and would like to find Q'. This is as easy as solving a properly defined quadratic equation on v (or u, doesn't really matter).

Luckily you can use CalculatePolygonPointST() in C++ (or GetPolyPointST() in Python) exactly for that. Please find attached scene that demonstrates exactly this.

sdk-quadrilateral.c4d

Cheers,
Ilia

Behavior:

The python tag code is:

import c4d

doc: c4d.documents.BaseDocument
op: c4d.BaseTag

def main() -> None:
    polyOp: c4d.BaseObject = op.GetObject()

    targetOp: c4d.BaseObject = polyOp.GetDown()
    lineUOp: c4d.BaseObject = targetOp.GetNext()
    lineVOp: c4d.BaseObject = lineUOp.GetNext()

    [A, B, C, D] = polyOp.GetAllPoints()

    P = targetOp.GetRelPos()
    
    hl = c4d.modules.hair.HairLibrary()
    U, V = hl.GetPolyPointST(P, A, B, C, D, True)

    AB, DC, AD, BC = B - A, C - D, D - A, C - B

    lineUOp.SetPoint(0, A + AB * U)
    lineUOp.SetPoint(1, D + DC * U)
    lineVOp.SetPoint(0, A + AD * V)
    lineVOp.SetPoint(1, B + BC * V)

In C++ the usage would be very similar:

PolygonObject* polyObjOp; // = static_cast<PolygonObject*>(op);
Vector* ptsPoly = polyObjOp->GetPointR();

Vector A = ptsPoly[0], B = ptsPoly[1], C = ptsPoly[2], D = ptsPoly[3];
Float U = 0, V = 0;
Vector P = pointOp->GetRelPos();

Bool inBoundary = maxon::GeometryUtilsInterface::CalculatePolygonPointST(P, A, B, C, D, true, U, V);

chuanzhen

@i_mazlov Thank you for your detailed reply, it was really helpful！
The following video shows the results, which work well. The blue polygon uses GetPolyPointST()

chuanzhen

Hi @i_mazlov
When using GetPolyPointST(), I encountered a problem where the positions of almost equal points (with only a slight difference in numerical accuracy) resulted in significantly different U and V results.
I am studying this page seriously, and it seems that some knowledge points are beyond my scope, so I need some time to explore.

The following are the position information of points a, b, c, d, and p. They only differ in accuracy, with u ≈ 0.25 and v ≈ 0.75 being the correct values .

-----The first group of points a, b, c, d, p**
Vector(-168.343, 0, 0) Vector(-33.643, -161.111, 0) Vector(168.343, 0, 0) Vector(33.643, 161.111, 0) Vector(16.822, 80.555, 0)
result
U:0.25000334744941866,V:0.7500002434367552

-----The second group of points a, b, c, d, p
a --
x: -168.34254079251306
y: 2.842170943040401e-14
z: 0.0
b --
-33.643095911746244
-161.11092598042546
0.0
c --
168.34254079251275
2.842170943040401e-14
0.0
d --
33.64309591174592
161.1109259804255
0.0
p --
16.82154795587286
80.55546299021273
0.0
result:
U:0.4999999999999992, V:0.7453163657465864

chuanzhen

Hi @i_mazlov
Through exploring these two pages，page1,page2, I have converted them into Python code Inverse_Bilinear_Interpolation(), which can replace GetPolyPointST(). I am not sure why GetPolyPointST() may obtain incorrect values in certain situations, but Inverse_Bilinear_Interpolation() does not go wrong in complex situations.

video use Inverse_Bilinear_Interpolation()

code:

import math


# algorithm_InverseBilinearInterpolation 逆双线性插值 的算法

def Wedge2D(v:c4d.Vector,w:c4d.Vector):
    return v[0] * w[1] - v[1] * w[0]

def Inverse_Bilinear_Interpolation(p:c4d.Vector,p0:c4d.Vector,p1:c4d.Vector,p3:c4d.Vector,p2:c4d.Vector):

    """
    注意对点顺序的定义
    Attention :  points Order

    p2 --- p3
    |      |
    |      |
    p0 --- p1

    """

    q = p - p0
    b1 = p1 - p0
    b2 = p2 - p0
    b3 = p0 - p1 - p2 + p3

    A = Wedge2D(b2,b3)
    B = Wedge2D(b3,q) - Wedge2D(b1,b2)
    C = Wedge2D(b1,q)
    # 求V - solve v
    if abs(A) < 0.00000001:
        v = -C / B
    else:
        discrim = B * B - 4.0 * A * C
        v = 0.5 * (-B - math.sqrt(discrim)) / A  # CCW
        #v = 0.5 * (-B + math.sqrt(discrim)) / A  # CW

    # 求 u - solve u
    denom = b1 + v * b3

    if abs(denom[0]) > abs(denom[1]):
        u = (q[0] - b2[0] * v) / denom[0]
    else:
        u = (q[1] - b2[1] * v) / denom[1]

    return u,v
# example use CPolygon point position a,b,c,d, and point p in quard(CPolygon)
# u,v = Inverse_Bilinear_Interpolation(p,a,b,c,d)

sample file: custom-quadrilateral.c4d

But for the code part, there are still some doubts about the difference between CCW and CW. When I understand what's going on, I will come back and reply. If any friend can explain, that's the best!

i_mazlov

Hi @chuanzhen,

sorry for the delayed answer.

Thank you for sharing your findings with us! This is a numerical issue in the GetPolyPointST() implementation. I've fixed it but it will only be available in one of our future releases.

Cheers,
Ilia

chuanzhen

@i_mazlov