Splines in 3d graphics, the most automated option

About a month ago, I began to learn Python from the book of Dawson and woke up deeply in the process of writing his game under pygame. TK was such that it seemed most promising to make a game with pseudo-three-dimensional graphics by stuffing the saved surfaces of 3d splines into sprites. I’ll write about the latter.



So, there are polygons (it is easiest to work with quadrangles) on which we want to stretch cubic surfaces so that they fit together quite smoothly - these surfaces are splines.







It is most convenient to present a spline for one polygon as a function





g(u,v)=v cdotg1(u)+(1v) cdotg2(u)+u cdotg3(v)+(1u) cdotg4(v)(1)







Here





g1,g2,g3,g4







- cubic polynomials satisfying some boundary conditions (in the figure below - light green and red curves, and derivative-initial conditions - lilac and blue vectors); u and v vary from 0 to 1.







In this interpretation, some degrees are lost (products of 2 and 3 degrees of parameters), but the polynomial coefficients can be found by specifying only 12 vectors of initial conditions (4 polygon points and derivative vectors with respect to u and v for each point). At the junctions, splines coincide if similar initial conditions are set for neighboring polygons (the derivative vectors must be collinear, the surface passes through the same points of the polygon).



One misfortune - the derivative with such a statement of the problem on the entire border may not coincide - there will be small artifacts at the junctions. You can think up 4 more conditions for correcting this and carefully add another term to the formula





u cdotv cdot(u1) cdot(v1) cdotg5(u,v)







which does not spoil the borders, but this is a topic for a separate article.



An alternative is the Bezier surface , but it proposes to set 16 parameters of incomprehensible (to me) mathematical meaning, i.e. it is assumed that the artist will work with his own hands. This is not very suitable for me, so a bicycle with crutches follows.



Coefficients (1) are most easily calculated by finding the inverse matrix to the matrix of values ​​and derivatives in the corners and multiplying by the input conditions (three times, for each coordinate). The matrix can be like this:



Hidden text
[[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], #g (0,0)

[1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0], #g (1,0)

[1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0], #g (0,1)

[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], #g (1,1)

[0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], # dg / du (0,0)

[0, 0, 0, 0, 1, 0, 0, 0, 2, 0, 3, 0], # dg / du (1,0)

[0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0], # dg / du (0,1)

[0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 3, 3], # dg / du (1,1)

[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], # dg / dv (0,0)

[0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1], # dg / dv (1,0)

[0, 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0], # dg / dv (0,1)

[0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 0, 1]] # dg / dv (1,1)



NumPy does an excellent job of this.



One question remains - where to get the vectors of derivatives. It is assumed that we must somehow choose them based on the position of neighboring (for the polygon) points and for reasons of smoothness.

For me personally, it was still completely counterintuitive how to deal with the length of the derivative vector. The direction is understandable, but the length?



As a result, the following algorithm was born:



1. At the first stage, some classification of points occurs. In the graph (which the polygon points and their connections specify), cycles of length 4 are searched and stored, and in addition, neighbors that are most suitable for the continuation of the edges of the polygon are recorded (it is determined in advance which edges correspond to a change in the parameter u and which correspond to v). Here is a piece of code that searches, sorts, and remembers neighbors for the 0th point of a cycle:



Hidden text
"""   : cykle[5] cykle[7] cykle[4]--cykle[0] --u-- cykle[1]-cykle[6] |v |v cykle[10]-cykle[3] --u-- cykle[2]-cykle[8] cykle[11] cykle[9] """ sosed = [] for i in range(0, N): if self.connects[i][ind] == 1 and i != cykle[0] and i != cykle[1] and i != cykle[3]: sosed += [i] #    ... if(len(sosed) < 2): continue #    else: #      cykle[0] vs0c3 = MinusVecs(self.dots[sosed[0]], self.dots[cykle[3]]) vs1c1 = MinusVecs(self.dots[sosed[1]], self.dots[cykle[1]]) #        vs0c1 = MinusVecs(self.dots[sosed[0]], self.dots[cykle[1]]) vs1c3 = MinusVecs(self.dots[sosed[1]], self.dots[cykle[3]]) Rvsc = [ScalMult(vs0c3,vs0c3),ScalMult(vs1c1,vs1c1),ScalMult(vs0c1,vs0c1),ScalMult(vs1c3,vs1c3) ] n1 = Rvsc.index(max(Rvsc)) if n1 < 2: cykle += [sosed[1],sosed[0]] #  u   v else: cykle += [sosed[0],sosed[1]]
      
      







Further, when constructing the spline, the derivative along the edge (for example, with respect to the parameter u) at the point of the polygon is selected based on the location of two points of the edge and one neighboring point (let it be the points vec1, vec2 and vec3; the point at which the derivative is sought is the second).







At first I tried to use the normalized vector vec3 - vec1 (like I applied the Lagrange theorem) for this role, but problems arose with the length of the vector of derivatives - making it a constant was a bad idea.



Lyrical digression:



For a small illustration of what the derivative vector is in the parametric version, we turn to a simple two-dimensional analogy - here is a piece of a circular arc:











g(u)=(Rsin(u pi/2),Rcos(u pi/2))











g(u)=(R pi/2 cdotcos(u pi/2),R pi/2 cdotsin(u pi/2))







Those. derivative modulus = R * pi / 2 and, generally speaking, depends on the size of the piece of the arc, which we set through the parameter.



What to do now with this? Leo Nikolaevich Tolstoy bequeathed to us that everything is round = good, therefore, if we set the derivative at a point as if we would like to build an arc there, we will get a smooth beautiful curve.



The end of the digression.



The second stage of the derivative search:



2. Through our three points vec1, vec2, vec3 we build a plane, in this plane we look for a circle that passes through all three points. The desired derivative will be directed along the tangent to the circle at vec2, and the modulus of the vector of the derivative should be equal to the product of the radius of the circle and the angle of the sector, which form two points of the polygon face (similar to our simple flat analogy from the lyrical digression).





g=R cdot phi







Everything seems to be simple, here is the function code, for example (NumPy is used again):



not a code, but ...
 def create_dd(vec1, vec2, tuda, vec3): """     1-2       3   == 1""" #0.         -         :( vecmult = VecMult(MinusVecs(vec2,vec1),MinusVecs(vec3,vec1)) #  vec2-vec1 x vec3 - vec1 if vecmult[0]*vecmult[0] + vecmult[1]*vecmult[1] + vecmult[2]*vecmult[2] < 0.000001: #  0 if tuda == 1: return NormVec(MinusVecs(vec3,vec2)) else: return NormVec(MinusVecs(vec1,vec2)) #1.   ,     #           2   M = np.zeros((3,3),float) C = np.zeros((3,1),float) M[0][0] = 2*vec2[0] - 2*vec1[0] M[0][1] = 2*vec2[1] - 2*vec1[1] M[0][2] = 2*vec2[2] - 2*vec1[2] C[0] = -vec1[0]*vec1[0] - vec1[1]*vec1[1] - vec1[2]*vec1[2] + vec2[0]*vec2[0] + vec2[1]*vec2[1] + vec2[2]*vec2[2] M[1][0] = 2*vec3[0] - 2*vec2[0] M[1][1] = 2*vec3[1] - 2*vec2[1] M[1][2] = 2*vec3[2] - 2*vec2[2] C[1] = -vec2[0]*vec2[0] - vec2[1]*vec2[1] - vec2[2]*vec2[2] + vec3[0]*vec3[0] + vec3[1]*vec3[1] + vec3[2]*vec3[2] #   ,   4     M[2][0] = vecmult[0] M[2][1] = vecmult[1] M[2][2] = vecmult[2] C[2] = vecmult[0]*vec1[0] + vecmult[1]*vec1[1] + vecmult[2]*vec1[2] # M *    = C Mobr = np.linalg.inv(M) Rvec = np.dot(Mobr,C) res = NormVec(VecMult(MinusVecs(Rvec,vec2), vecmult)) #    ka #R = math.sqrt((Rvec[0]-vec1[0])**2 + (Rvec[1]-vec1[1])**2 + (Rvec[2]-vec1[2])**2) #   = 3.14*2*asin(l/(2*R))/180 * R,      l  ,    l = 1.2*math.sqrt((vec1[0]-vec2[0])**2 + (vec1[1]-vec2[1])**2 + (vec1[2]-vec2[2])**2) if ScalMult(res, MinusVecs(vec2, vec1)) > 0: #     return MultVxC(res, tuda*l) else: return MultVxC(res, -tuda*l)
      
      







Well, in general, it all somehow works. For demonstration, I took a 5x5x5 cube and changed the position of the points with a randomizer:



GIF, 8Mb





All Articles