Aplicatii ale geometriei în informatică - Laboratoare¶

Lab 1¶

În lista de probleme de mai jos, $\Delta ABC$ are vârfurile $A(1,1)$ , $B(4,1)$ , $C(2,3)$ .

Problema 2¶

Determinați imaginea triunghiului $ABC$ printr-o scalare uniformă de factor de scală $2$ relativ la punctul $Q(2,2)$ .

Soluție¶

Forma generală a matricei de scalare în jurul unui punct $Q$ :

S_{Q} (s_{x}, s_{y}) = (\begin{matrix} s_{x} & 0 & (1 - s_{x}) x_{Q} \\ 0 & s_{y} & (1 - s_{y}) y_{Q} \\ 0 & 0 & 1 \end{matrix})

$S_Q(s_x, s_y) = \begin{pmatrix} s_x & 0 & (1-s_x)x_Q \\ 0 & s_y & (1-s_y)y_Q \\ 0 & 0 & 1 \end{pmatrix}$

Concret, $Q=(2,2), s_x=s_y=2$ . Matricea transformării este:

T = S_{Q} (2, 2) = (\begin{matrix} 2 & 0 & - 2 \\ 0 & 2 & - 2 \\ 0 & 0 & 1 \end{matrix})

$T = S_Q(2,2) = \begin{pmatrix} 2 & 0 & -2 \\ 0 & 2 & -2 \\ 0 & 0 & 1 \end{pmatrix}$

Imaginea $\Delta ABC$ după scalare este:

[\begin{matrix} A^{'} & B^{'} & C^{'} \end{matrix}] = T \cdot [\begin{matrix} A & B & C \end{matrix}]

$\begin{bmatrix} A' & B' & C' \end{bmatrix} = T \cdot \begin{bmatrix} A & B & C \end{bmatrix}$

= (\begin{matrix} 2 & 0 & - 2 \\ 0 & 2 & - 2 \\ 0 & 0 & 1 \end{matrix}) \cdot (\begin{matrix} 1 & 4 & 2 \\ 1 & 1 & 3 \\ 1 & 1 & 1 \end{matrix}) = (\begin{matrix} 0 & 6 & 2 \\ 0 & 0 & 4 \\ 1 & 1 & 1 \end{matrix})

$= \begin{pmatrix} 2 & 0 & -2 \\ 0 & 2 & -2 \\ 0 & 0 & 1 \end{pmatrix} \cdot \begin{pmatrix} 1 & 4 & 2 \\ 1 & 1 & 3 \\ 1 & 1 & 1 \end{pmatrix} = \begin{pmatrix} 0 & 6 & 2 \\ 0 & 0 & 4 \\ 1 & 1 & 1 \end{pmatrix}$

Rezultă $A'(0,0)$ , $B'(6,0)$ , $C'(2,4)$ .

Problema 4¶

Determinați imaginea triunghiului $ABC$ printr-o forfecare de unghi $45^\circ$ , relativ la punctul $Q(2, 2)$ , în direcția vectorului $v(2, 1)$ .

Soluție¶

Forma generală a matricei de forfecare $Q$ , de unghi $\theta$ , în direcția versorului $\bar{\mathrm{v}}$ , este

S h e a r (Q, \bar{v}, θ) = (\begin{matrix} 1 - \tan θ \cdot v_{1} v_{2} & - \tan θ \cdot v_{1}^{2} & \tan θ \cdot v_{1} (v_{1} Q_{2} - v_{2} Q_{1}) \\ - \tan θ \cdot v_{2}^{2} & 1 + \tan θ \cdot v_{1} v_{2} & \tan θ \cdot v_{2} (v_{1} Q_{2} - v_{2} Q_{1}) \\ 0 & 0 & 1 \end{matrix})

$Shear(Q, \bar{\mathrm{v}}, \theta) = \begin{pmatrix} 1 - \tan \theta \cdot \mathrm{v}_1 \mathrm{v}_2 & - \tan \theta \cdot \mathrm{v}_1^2 & \tan \theta \cdot \mathrm{v}_1( \mathrm{v}_1 Q_2 - \mathrm{v}_2 Q_1) \\ - \tan \theta \cdot \mathrm{v}_2^2 & 1 + \tan \theta \cdot \mathrm{v}_1 \mathrm{v}_2 & \tan \theta \cdot \mathrm{v}_2( \mathrm{v}_1 Q_2 - \mathrm{v}_2 Q_1) \\ 0 & 0 & 1 \end{pmatrix}$

În particular, $Q=(2,2)$ , $\bar{\mathrm{v}} = \left( \frac{v1}{|v|}, \frac{v_2}{|v|} \right)=\left(\frac{2}{\sqrt{5}}, \frac{1}{\sqrt{5}} \right), \tan{\theta}=1$

Matricea transformării este:

T = S h e a r (Q, \bar{v}, 45^{\circ}) = (\begin{matrix} \frac{3}{5} & \frac{4}{5} & - \frac{4}{5} \\ - \frac{1}{5} & \frac{7}{5} & - \frac{2}{5} \\ 0 & 0 & 1 \end{matrix})

$T = Shear(Q, \bar{\mathrm{v}}, 45^\circ) = \begin{pmatrix} \frac{3}{5} & \frac{4}{5} & -\frac{4}{5}\\ -\frac{1}{5} & \frac{7}{5} & -\frac{2}{5} \\ 0 & 0 & 1 \end{pmatrix}$

Imaginea $\Delta ABC$ după forfecare este:

[\begin{matrix} A^{'} & B^{'} & C^{'} \end{matrix}] = T \cdot [\begin{matrix} A & B & C \end{matrix}]

$\begin{bmatrix} A' & B' & C' \end{bmatrix} = T \cdot \begin{bmatrix} A & B & C \end{bmatrix}$

= \frac{1}{5} (\begin{matrix} 3 & 4 & - 4 \\ - 1 & 7 & - 2 \\ 0 & 0 & 5 \end{matrix}) \cdot (\begin{matrix} 1 & 4 & 2 \\ 1 & 1 & 3 \\ 1 & 1 & 1 \end{matrix}) = (\begin{matrix} \frac{3}{5} & \frac{12}{5} & \frac{14}{5} \\ \frac{4}{5} & \frac{1}{5} & \frac{17}{5} \\ 1 & 1 & 1 \end{matrix})

$= \frac{1}{5}\begin{pmatrix} 3 & 4 & -4\\ -1 & 7 & -2 \\ 0 & 0 & 5 \end{pmatrix} \cdot \begin{pmatrix} 1 & 4 & 2 \\ 1 & 1 & 3 \\ 1 & 1 & 1 \end{pmatrix} = \begin{pmatrix} \frac{3}{5} & \frac{12}{5} & \frac{14}{5} \\ \frac{4}{5} & \frac{1}{5} & \frac{17}{5} \\ 1 & 1 & 1 \end{pmatrix}$

Rezultă $A'\left(\frac{3}{5},\frac{4}{5}\right)$ , $B'\left(\frac{12}{5},\frac{1}{5}\right)$ , $C'\left(\frac{14}{5},\frac{17}{5}\right)$

Problema 5¶

Aplicați forfecarea de la problema precedentă pătratului $ABCD$ , cu $A(0,0)$ , $B(3, 0)$ , $C(3, 3)$ și $D(0,3)$ .

Soluție¶

Aplicăm transformarea $T$ găsită la problema 4 vârfurilor pătratului $ABCD$ :

[\begin{matrix} A^{'} & B^{'} & C^{'} & D^{'} \end{matrix}] = T \cdot [\begin{matrix} A & B & C & D \end{matrix}]

$\begin{bmatrix} A' & B' & C' & D' \end{bmatrix} = T \cdot \begin{bmatrix} A & B & C & D \end{bmatrix}$

= [\begin{matrix} - \frac{4}{5} & 1 & \frac{17}{5} & \frac{8}{5} \\ - \frac{2}{5} & - 1 & \frac{16}{5} & \frac{19}{5} \\ 1 & 1 & 1 & 1 \end{matrix}]

$= \begin{bmatrix} -\frac{4}{5} & 1 & \frac{17}{5} & \frac{8}{5} \\ -\frac{2}{5} & -1 & \frac{16}{5} & \frac{19}{5} \\ 1 & 1 & 1 & 1 \end{bmatrix}$

Obținem $A'\left(-\frac{4}{5}, -\frac{2}{5}\right)$ , $B'\left(1, -1\right)$ , $C'\left(\frac{17}{5}, \frac{16}{5}\right)$ , $D'\left(\frac{8}{5}, \frac{19}{5}\right)$

Problema 7¶

Determinați imaginea triunghiului $ABC$ prin reflexia relativ la dreapta $2x + 3y − 5 = 0$ .

Soluție¶

Forma generală a matricei de reflexie față de o dreaptă $\Delta : ax+by+c=0$ este:

R_{Δ} = (\begin{matrix} \frac{b^{2} - a^{2}}{a^{2} + b^{2}} & - \frac{2 a b}{a^{2} + b^{2}} & - \frac{2 a c}{a^{2} + b^{2}} \\ - \frac{2 a b}{a^{2} + b^{2}} & - \frac{b^{2} - a^{2}}{a^{2} + b^{2}} & - \frac{2 b c}{a^{2} + b^{2}} \\ 0 & 0 & 1 \end{matrix})

$R_\Delta = \begin{pmatrix} \frac{b^2-a^2}{a^2+b^2} & -\frac{2ab}{a^2+b^2} & -\frac{2ac}{a^2+b^2} \\ -\frac{2ab}{a^2+b^2} & - \frac{b^2-a^2}{a^2+b^2} & -\frac{2bc}{a^2+b^2} \\ 0 & 0 & 1 \end{pmatrix}$

În particular, pentru dreapta $\Delta : 2x+3y-5=0$ , matricea transformării este:

T = R_{Δ} = (\begin{matrix} \frac{5}{13} & - \frac{12}{13} & \frac{20}{13} \\ - \frac{20}{13} & - \frac{5}{13} & \frac{30}{13} \\ 0 & 0 & 1 \end{matrix})

$T = R_\Delta = \begin{pmatrix} \frac{5}{13} & -\frac{12}{13} & \frac{20}{13} \\ -\frac{20}{13} & -\frac{5}{13} & \frac{30}{13} \\ 0 & 0 & 1 \end{pmatrix}$

Imaginea $\Delta ABC$ după transformare este:

[\begin{matrix} A^{'} & B^{'} & C^{'} \end{matrix}] = T \cdot [\begin{matrix} A & B & C \end{matrix}]

$\begin{bmatrix} A' & B' & C' \end{bmatrix} = T \cdot \begin{bmatrix} A & B & C \end{bmatrix}$

= (\begin{matrix} 1 & \frac{28}{13} & \frac{- 6}{13} \\ 1 & - \frac{23}{13} & \frac{- 9}{13} \\ 1 & 1 & 1 \end{matrix})

$= \begin{pmatrix} 1 & \frac{28}{13} & \frac{-6}{13} \\ 1 & -\frac{23}{13} & \frac{-9}{13} \\ 1 & 1 & 1 \end{pmatrix}$

Obținem $A'\left(1,1\right)$ , $B'\left(\frac{28}{13},-\frac{23}{13}\right)$ , $C'\left(-\frac{6}{13}, -\frac{9}{13}\right)$

Problema 10¶

Determinați imaginea triunghiului $ABC$ prin rotația cu $90^\circ$ în jurul punctului C, urmată de reflexia relativ la dreapta $AB$ .

Soluție¶

Forma generală a matricei de rotație în jurul unui punct $P$ cu unghi $\theta$ :

R o t_{P} (θ) = (\begin{matrix} \cos θ & - \sin θ & x_{P} (1 - \cos θ) + y_{P} \sin θ \\ \sin θ & \cos θ & - x_{P} \sin θ + y_{P} (1 - \cos θ) \\ 0 & 0 & 1 \end{matrix})

$Rot_P(\theta)=\begin{pmatrix} \cos\theta & -\sin\theta & x_P (1-\cos\theta)+ y_P \sin\theta \\ \sin\theta & \cos\theta & - x_P \sin\theta+ y_P (1-\cos\theta) \\ 0 & 0 & 1 \end{pmatrix}$

Calculăm imaginea $\Delta ABC$ după rotația cu $90^\circ$ în jurul punctului $C$ :

T = R o t_{C} (90^{\circ}) = (\begin{matrix} 0 & - 1 & 5 \\ 1 & 0 & 1 \\ 0 & 0 & 1 \end{matrix})

$T = Rot_C(90^\circ) = \begin{pmatrix} 0 & -1 & 5 \\ 1 & 0 & 1 \\ 0 & 0 & 1 \end{pmatrix}$

[\begin{matrix} A^{'} & B^{'} & C^{'} \end{matrix}] = T \cdot [\begin{matrix} A & B & C \end{matrix}] = (\begin{matrix} 4 & 4 & 2 \\ 2 & 5 & 3 \\ 1 & 1 & 1 \end{matrix})

$\begin{bmatrix} A' & B' & C' \end{bmatrix} = T \cdot \begin{bmatrix} A & B & C \end{bmatrix} = \begin{pmatrix} 4 & 4 & 2 \\ 2 & 5 & 3 \\ 1 & 1 & 1 \end{pmatrix}$

Rezultă $A'(4,2)$ , $B'(4,5)$ , $C'(2,3)$ .

Forma generală a matricei de reflexie față de o dreaptă $\Delta : ax+c=0$ este

R_{Δ} = (\begin{matrix} - 1 & 0 & - 2 \frac{c}{a} \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix})

$R_\Delta = \begin{pmatrix} -1& 0 & -2\frac{c}{a} \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}$

În cazul nostru, $\Delta = A'B'$ , $a=1$ , $c=-4$ , de unde matricea de transformare este:

T = R_{Δ} = (\begin{matrix} - 1 & 0 & 8 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix})

$T = R_\Delta = \begin{pmatrix} -1 & 0 & 8 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}$

Imaginea finală a $\Delta ABC$ după transformări este

[\begin{matrix} A_{10} & B_{10} & C_{10} \end{matrix}] = T \cdot [\begin{matrix} A^{'} & B^{'} & C^{'} \end{matrix}] = (\begin{matrix} 4 & 4 & 6 \\ 2 & 5 & 3 \\ 1 & 1 & 1 \end{matrix})

$\begin{bmatrix} A_{10} & B_{10} & C_{10} \end{bmatrix} = T \cdot \begin{bmatrix} A' & B' & C' \end{bmatrix} = \begin{pmatrix} 4 & 4 & 6 \\ 2 & 5 & 3 \\ 1 & 1 & 1 \end{pmatrix}$ Obținem triunghiul

Δ A_{10} B_{10} C_{10}

$\Delta A_{10}B_{10}C_{10}$ ,

A_{10} (4, 2)

$A_{10}(4,2)$ ,

B_{10} (4, 5)

$B_{10}(4,5)$ ,

C_{10} (6, 3)

$C_{10}(6,3)$

Lab 2 Info¶

Scrieți un program în Python care să determine imaginea unui poliedru printr-o reflexie față de un plan dat. Considerați atât planul dat prin ecuația generală (datele de intrare vor fi coeficienții ecuației), cât și planul dat printr-un punct și vectorul normal. În cel de-al doilea caz, scrieți mai întâi ecuația generala, plecând de la punct și vectorul normal. Programul trebuie să conțină următorii pași:

determinarea punctului de intersecție dintre plan și una dintre axe;
determinarea matricii translației care face planul sa treacă prin origine;
determinerea matricii rotației față de o axă care face ca normala la plan sâ fie paralelă cu un plan de coordonate;
determinarea matricii rotației față de altă axă care face ca planul să coincidă cu un plan de coordonate;
determinarea matricii reflexiei față de planul de la punctul precedent;
determinarea matricii rotației inverse rotației de la punctul 4;
determinarea matricii rotației inverse celei de la punctul 3;
determinarea matricii translației inverse celei de la punctul 2;
determinarea matricii transformării, prin înmulțirea celor 7 matrici precedente.
citirea vârfurilor poliedrului;
stabilirea matricii omogene a coordonatelor vârfurilor poliedrului transformat.

Soluție¶

Definim o clasă pentru matricea omogenă în 3D:

from __future__ import annotations
from math import sqrt, cos, sin


class Matrix4:
    def __init__(self, *args):
        self.items = [0]*16
        for i in range(min(16, len(args))):
            self.items[i]=args[i]

    def __getitem__(self, pos):
        i, j = pos
        return self.items[4 * i + j]

    def __setitem__(self, pos, value):
        i, j=pos
        self.items[4 * i + j]=value

    def __str__(self):
        s=""
        for i in range(4):
            s += str("".join(list(map(lambda x:f"{x:.3f}".rjust(8), self.items[4*i:4*i+4]))))+"\n"
        return s

    def _scalar_mul_(self, x:float) -> Matrix4: return Matrix4([x*i for i in self.items])

    def _matrix4_mul_(self, m:Matrix4) -> Matrix4:
        r = Matrix4()
        for i in range(4):
            for j in range(4):
                for k in range(4):
                    r[i,j]+=self[i,k]*m[k,j]
        return r

    def _vec4_mul(self, v:tuple[float]):
        r = [0]*4
        v = (v[0], v[1], v[2], 1)
        for i in range(4):
            for j in range(4):
                r[i] += self[i,j]*v[j]
        return tuple(r)

    def __mul__(self, other):
        if isinstance(other, float) or isinstance(other, int): return self._scalar_mul_(other)
        if isinstance(other, Matrix4): return self._matrix4_mul_(other)
        if isinstance(other, tuple): return self._vec4_mul(other)
        raise Exception("Unsupported operation: invalid multiplication type")

    def __rmul__(self, other):
        if isinstance(other, float) or isinstance(other, int): return self._scalar_mul_(other)
        if isinstance(other, Matrix4): return other._matrix4_mul_(self)
        raise Exception("Unsupported operation: invalid multiplication type")

    def apply(self, pts): return [self * p for p in pts]

Definim principalele tipuri de matrici de transformare:

def TranslationMatrix(p):
    x, y, z = p
    return Matrix4(1,0,0,x, 0,1,0,y, 0,0,1,z, 0,0,0,1)

def CosSinRotXMatrix(c,s): return Matrix4(1,0,0,0, 0,c,-s,0, 0,s,c,0, 0,0,0,1)
def RotationXMatrix(thetaX): return CosSinRotXMatrix(cos(thetaX), sin(thetaX))

def CosSinRotYMatrix(c,s): return Matrix4(c,0,s,0, 0,1,0,0, -s,0,c,0, 0,0,0,1)
def RotationYMatrix(thetaY): return CosSinRotYMatrix(cos(thetaY), sin(thetaY))

def CosSinRotZMatrix(c,s): return Matrix4(c,-s,0,0, s,c,0,0, 0,0,1,0, 0,0,0,1)
def RotationZMatrix(thetaZ): return CosSinRotZMatrix(cos(thetaZ), sin(thetaZ))

def ReflectYZMatrix(): return Matrix4(-1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1)

Câteva metode ajutătoare:

"""
    in p: planul (a,b,c,d)
    out: (point, v), point = punctul de intersectie, v = versor Ox, Oy sau Oz
"""
def intersectAxis(p):
    a,b,c,d = p
    if a!=0: return ((-d/a,0,0), (1,0,0))
    if b!=0: return ((0,-d/b,0), (0,1,0))
    if c!=0: return ((0,0,-d/c), (0,0,1))
    if d!=0: raise Exception("Runtime error: plan invalid d=0, d!=0")
    return ((0, 0, 0), (1, 0, 0)) # a=b=c=0, d=0 triat la input

def neg(p):
    x,y,z = p
    return (-x,-y,-z)

Definirea planului $ax+by+cz+d=0$ (a se vedea LiviuStefan_NeacsuMiclea_235_Lab2Info.py pentru citirea planului de la tastatura, in ambele variante precizate de cerinta):

plane = (2,1,3,-10) 
a,b,c,d = plane

Calcularea matricilor de transformare, conform pașilor precizați:

transforms = []
inv_transforms = []

print("Calculam intersectia cu una dintre axe")
i_point, i_axis = intersectAxis(plane)
print(f"Planul intersecteaza axa data de versorul {i_axis} in punctul {i_point}")

print("Determinam matricea de translatie care face planul sa treaca prin origine:")

if d!=0:
    translate_o = TranslationMatrix(neg(i_point))
    print(translate_o)
    transforms.append(translate_o)
    inv_transforms.append(TranslationMatrix(i_point))
else:
    print("Planul este deja in origine")

print("Determinam matricea de rotatie fata de o axa care face ca normala \n"
      "la plan sa fie paralela cu un plan de coordonate:")
print(">> Rotatie in jurul axei Ox pentru a face normala || cu planul xOy aka perp. pe axa Oz")

if c!=0:
    norm_x = sqrt(b**2+c**2)
    cos_tx = b/norm_x
    sin_tx = -c/norm_x
    rot_x = CosSinRotXMatrix(cos_tx, sin_tx)
    print(rot_x)
    transforms.append(rot_x)
    inv_transforms.append(CosSinRotXMatrix(cos_tx, -sin_tx))
else:
    rot_x = RotationXMatrix(0)
    print("Normala este deja paralela cu axa Oz")

print("Determinam matricea de rotatie fata de alta axa care face ca planul \n"
      "sa coincida cu un plan de coordonate")
print(">> Rotatie in jurul axei Oz pentru a face planul || cu planul yOz")

new_a, new_b,_,_ = rot_x.apply([(a,b,c)])[0]

if new_b != 0:
    norm_z = sqrt(new_a**2+new_b**2)
    cos_tz = new_a/norm_z
    sin_tz = -new_b/norm_z
    rot_z = CosSinRotZMatrix(cos_tz, sin_tz)
    print(rot_z)
    transforms.append(rot_z)
    inv_transforms.append(CosSinRotZMatrix(cos_tz, -sin_tz))
else:
    print("Planul este deja paralel cu axa Oz")


print("Planul se transforma in planul yOz")
print("Matricea de reflexie fata de YZ:")

refl_yz = ReflectYZMatrix()
print(refl_yz)
transforms.append(refl_yz)

print("Matricile inverse transformarilor planului (in ordinea inversa mentionarii lor)")

for m in inv_transforms[::-1]:
    print(str(m)+"\n")

transforms = (transforms+ inv_transforms[::-1])[::-1]

full_t = Matrix4(1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1)
for t in transforms:
    full_t *= t

print("Matricea transformarii compuse")
print(full_t)

print()

Calculam intersectia cu una dintre axe
Planul intersecteaza axa data de versorul (1, 0, 0) in punctul (5.0, 0, 0)
Determinam matricea de translatie care face planul sa treaca prin origine:
   1.000   0.000   0.000  -5.000
   0.000   1.000   0.000   0.000
   0.000   0.000   1.000   0.000
   0.000   0.000   0.000   1.000

Determinam matricea de rotatie fata de o axa care face ca normala 
la plan sa fie paralela cu un plan de coordonate:
>> Rotatie in jurul axei Ox pentru a face normala || cu planul xOy aka perp. pe axa Oz
   1.000   0.000   0.000   0.000
   0.000   0.316   0.949   0.000
   0.000  -0.949   0.316   0.000
   0.000   0.000   0.000   1.000

Determinam matricea de rotatie fata de alta axa care face ca planul 
sa coincida cu un plan de coordonate
>> Rotatie in jurul axei Oz pentru a face planul || cu planul yOz
   0.535   0.845   0.000   0.000
  -0.845   0.535   0.000   0.000
   0.000   0.000   1.000   0.000
   0.000   0.000   0.000   1.000

Planul se transforma in planul yOz
Matricea de reflexie fata de YZ:
  -1.000   0.000   0.000   0.000
   0.000   1.000   0.000   0.000
   0.000   0.000   1.000   0.000
   0.000   0.000   0.000   1.000

Matricile inverse transformarilor planului (in ordinea inversa mentionarii lor)
   0.535  -0.845   0.000   0.000
   0.845   0.535   0.000   0.000
   0.000   0.000   1.000   0.000
   0.000   0.000   0.000   1.000


   1.000   0.000   0.000   0.000
   0.000   0.316  -0.949   0.000
   0.000   0.949   0.316   0.000
   0.000   0.000   0.000   1.000


   1.000   0.000   0.000   5.000
   0.000   1.000   0.000   0.000
   0.000   0.000   1.000   0.000
   0.000   0.000   0.000   1.000


Matricea transformarii compuse
   0.429  -0.286  -0.857   2.857
  -0.286   0.857  -0.429   1.429
  -0.857  -0.429  -0.286   4.286
   0.000   0.000   0.000   1.000

Testarea transformării pe un poliedru generat care are vârfurile pe sfera unitate:

def gen_sphere(p, r):
    pts = []
    x,y,z=p
    for t in range(0, 4):
        at = 6.28 * t/4
        for p in range(0, 4):
            ap = 3.14 * p/4
            xx = x+r*cos(at)*cos(ap)
            yy = y + r * cos(at) * sin(ap)
            zz = z + r * sin(at)
            pts.append((xx,yy,zz))
    return pts

pts = gen_sphere((0, 0, 0), 1)

print("Matricea coordonatelor ale poliedrului initial")
for j in range(3):
    for i in range(len(pts)):
        print(f"{pts[i][j]:.2f}".rjust(5,' '), end=" ")
    print("")
for i in range(len(pts)):
    print(f"{1}".rjust(5, ' '), end=" ")
print("\n")

Matricea coordonatelor ale poliedrului initial
 1.00  0.71  0.00 -0.71  0.00  0.00  0.00 -0.00 -1.00 -0.71 -0.00  0.71 -0.00 -0.00 -0.00  0.00 
 0.00  0.71  1.00  0.71  0.00  0.00  0.00  0.00  0.00 -0.71 -1.00 -0.71  0.00 -0.00 -0.00 -0.00 
 0.00  0.00  0.00  0.00  1.00  1.00  1.00  1.00  0.00  0.00  0.00  0.00 -1.00 -1.00 -1.00 -1.00 
    1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1

pts_t = full_t.apply(pts)
print("Matricea coordonatelor poliedrului transformat")
for j in range(3):
    for i in range(len(pts_t)):
        print(f"{pts_t[i][j]:.2f}".rjust(5,' '), end=" ")
    print("")
for i in range(len(pts_t)):
    print(f"{1}".rjust(5, ' '), end=" ")
print("\n")

Matricea coordonatelor poliedrului transformat
 3.29  2.96  2.57  2.35  2.00  2.00  2.00  2.00  2.43  2.75  3.14  3.36  3.71  3.71  3.71  3.72 
 1.14  1.83  2.29  2.24  1.00  1.00  1.00  1.00  1.71  1.02  0.57  0.62  1.86  1.86  1.86  1.86 
 3.43  3.38  3.86  4.59  4.00  4.00  4.00  4.00  5.14  5.19  4.71  3.98  4.57  4.57  4.57  4.57 
    1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1

Reprezentarea grafică a transformării:

import matplotlib.pyplot as plt
import numpy as np

def plot_axes(ax):
    x, y, z = np.zeros((3, 1))
    u, v, w = np.array([5, 0, 0])
    ax.quiver(x, y, z, u, v, w, arrow_length_ratio=0.1, color='r')
    u, v, w = np.array([0, 5, 0])
    ax.quiver(x, y, z, u, v, w, arrow_length_ratio=0.1, color='g')
    u, v, w = np.array([0, 0, 5])
    ax.quiver(x, y, z, u, v, w, arrow_length_ratio=0.1, color='b')


def plot_points(ax, pts, m="o", c="r"):
    for p in pts:
        ax.scatter(p[0], p[1], p[2], marker=m, color=c)


def plot_surface(ax, pts):
    x = np.array([p[0] for p in pts])
    y = np.array([p[1] for p in pts])
    z = np.array([p[2] for p in pts])
    ax.scatter(x, y, z)


def plot_plane(fig, p):
    a, b, c, d = p
    if c != 0:
        x = np.linspace(-5, 5, 10)
        y = np.linspace(-5, 5, 10)
        X, Y = np.meshgrid(x, y)
        Z = (-d - a * X - b * Y) / c
    elif b != 0:
        x = np.linspace(-5, 5, 10)
        z = np.linspace(-5, 5, 10)
        X, Z = np.meshgrid(x, z)
        Y = (-d - a * X) / b
    elif a != 0:
        y = np.linspace(-5, 5, 10)
        z = np.linspace(-5, 5, 10)
        Y, Z = np.meshgrid(y, z)
        X = (-d - 0 * Y - 0 * Z) / a
    else:
        y = np.linspace(-5, 5, 10)
        z = np.linspace(-5, 5, 10)
        Y, Z = np.meshgrid(y, z)
        X = (0 - 0 * Y - 0 * Z)

    ax.plot_surface(X, Y, Z)

fig = plt.figure()
ax = fig.add_subplot(projection='3d')
plot_plane(ax, plane)

plot_points(ax, pts, c="r")
plot_points(ax, pts_t, c="b")

plot_axes(ax)
ax.view_init(-10, -30)
plt.show()

Lab 3 Info¶

Problema 1¶

Consideram curba Bézier cubică dată de punctele de control $b_0(−3, 1)$ , $b_1(−4, 4)$ , $b_2(4, 4)$ , $b_3(3, 1)$ .

Determinați punctele de control ale curbei, privită ca o curbă de gradul 4.

Reprezentați grafic pe același sistem de axe: curba, punctele de control și poligonul de control ale curbei privite ca o curbă de gradul 3, precum și punctele de control și poligonul de control ale curbei, privite ca o curbă de gradul 4.

Soluție¶

1.Pentru curba Bézier de gradul $n$ dată de punctele de control $b_0,\ldots,b_n$ , se poate construi curba identică de grad $n+1$ , dată de punctele de control $c_0,\ldots,c_{n+1}$ , astfel:

\begin{aligned} c_{0} & = b_{0} \\ c_{k} & = \frac{n + 1 - k}{n + 1} b_{k} + \frac{k}{n + 1} b_{k + 1}, \forall k = 1, \dots, n \\ c_{n + 1} & = b_{n} \end{aligned}

$\begin{align} c_0 &= b_0 \\ c_k &= \frac{n+1-k}{n+1}b_k + \frac{k}{n+1}b_{k+1}, \forall k=1, \ldots, n \\ c_{n+1} &= b_{n} \end{align}$

În particular, pentru $n=3$ , $b_0(−3, 1)$ , $b_1(−4, 4)$ , $b_2(4, 4)$ , $b_3(3, 1)$ , găsim:

\begin{aligned} c_{0} & = (- 3, 1) \\ c_{1} & = (- \frac{15}{4}, \frac{13}{4}) \\ c_{2} & = (0, 4) \\ c_{3} & = (\frac{15}{4}, \frac{13}{4}) \\ c_{4} & = (3, 1) \end{aligned}

$\begin{align} c_0 &= (-3, 1) \\ c_1 &= (-\frac{15}{4}, \frac{13}{4}) \\ c_2 &= (0, 4) \\ c_3 &= (\frac{15}{4}, \frac{13}{4}) \\ c_4 &= (3, 1) \end{align}$

2.Reprezentare grafică:

Problema 2¶

Considerăm curba Bézier de gradul 3 de la problema precedentă.

Divizați curba la valoarea $1/3$ a parametrului, determinând punctele de control și parametrizările pe intervalul $[0, 1]$ pentru fiecare dintre cele două segmente.

Reprezentați grafic pe același sistem de axe: cele doua segmente de curbă (cu culori diferite), punctele de control și poligonul de control pentru curba originală, precum și punctele de control și poligoanele de control pentru cele doua segmente obținute prin divizare.

Soluție¶

1.Pentru curba Bézier dată de punctele de control $b_0,\dots,b_n$ și un interval $[a,b]\subseteq[0,1]$ , curba Bézier care descrie curba segmentată pe intervalul $[a,b]$ are punctele de control

c_{k} = C (a)^{n - k} \cdot C (b)^{k} \cdot (\begin{matrix} b_{0} \\ \dots \\ b_{n} \end{matrix})

$c_k = C(a)^{n-k} \cdot C(b)^k \cdot \begin{pmatrix} b_0 \\ \cdots \\ b_n \end{pmatrix}$ Aplicăm pentru

n = 3

$n=3$ și

b_{0} (- 3, 1)

$b_0(−3, 1)$ ,

b_{1} (- 4, 4)

$b_1(−4, 4)$ ,

b_{2} (4, 4)

$b_2(4, 4)$ ,

b_{3} (3, 1)

$b_3(3, 1)$ , cu intervalele

[a, b] = [0, 1 / 3]

$[a,b]=[0,1/3]$ , respectiv

[a, b] = [1 / 3, 1]

$[a,b]=[1/3, 1]$ .

Rescriem ecuațiile lui $c_k$ din teoremă particularizate pentru $n=3$ :

\begin{aligned} c_{0} & = C (a)^{3} \cdot C (b)^{0} \cdot (\begin{matrix} b_{0} \\ b_{1} \\ b_{2} \\ b_{3} \end{matrix}) = C_{1} (a) C_{2} (a) C_{3} (a) (\begin{matrix} b_{0} \\ b_{1} \\ b_{2} \\ b_{3} \end{matrix}) \\ = [\begin{matrix} 1 - a & a \end{matrix}] \cdot [\begin{matrix} 1 - a & a & 0 \\ 0 & 1 - a & a \end{matrix}] \cdot [\begin{matrix} 1 - a & a & 0 & 0 \\ 0 & 1 - a & a & 0 \\ 0 & 0 & 1 - a & a \end{matrix}] \cdot (\begin{matrix} b_{0} \\ b_{1} \\ b_{2} \\ b_{3} \end{matrix}) \\ c_{1} & = C (a)^{2} \cdot C (b)^{1} \cdot (\begin{matrix} b_{0} \\ b_{1} \\ b_{2} \\ b_{3} \end{matrix}) = C_{1} (a) C_{2} (a) C_{3} (b) (\begin{matrix} b_{0} \\ b_{1} \\ b_{2} \\ b_{3} \end{matrix}) \\ = [\begin{matrix} 1 - a & a \end{matrix}] \cdot [\begin{matrix} 1 - a & a & 0 \\ 0 & 1 - a & a \end{matrix}] \cdot [\begin{matrix} 1 - b & b & 0 & 0 \\ 0 & 1 - b & b & 0 \\ 0 & 0 & 1 - b & b \end{matrix}] \cdot (\begin{matrix} b_{0} \\ b_{1} \\ b_{2} \\ b_{3} \end{matrix}) \\ c_{2} & = C (a)^{1} \cdot C (b)^{2} \cdot (\begin{matrix} b_{0} \\ b_{1} \\ b_{2} \\ b_{3} \end{matrix}) = C_{1} (a) C_{2} (b) C_{3} (b) (\begin{matrix} b_{0} \\ b_{1} \\ b_{2} \\ b_{3} \end{matrix}) \\ = [\begin{matrix} 1 - a & a \end{matrix}] \cdot [\begin{matrix} 1 - b & b & 0 \\ 0 & 1 - b & b \end{matrix}] \cdot [\begin{matrix} 1 - b & b & 0 & 0 \\ 0 & 1 - b & b & 0 \\ 0 & 0 & 1 - b & b \end{matrix}] \cdot (\begin{matrix} b_{0} \\ b_{1} \\ b_{2} \\ b_{3} \end{matrix}) \\ c_{3} & = C (b)^{1} \cdot C (b)^{2} \cdot (\begin{matrix} b_{0} \\ b_{1} \\ b_{2} \\ b_{3} \end{matrix}) = C_{1} (b) C_{2} (b) C_{3} (b) (\begin{matrix} b_{0} \\ b_{1} \\ b_{2} \\ b_{3} \end{matrix}) \\ = [\begin{matrix} 1 - b & b \end{matrix}] \cdot [\begin{matrix} 1 - b & b & 0 \\ 0 & 1 - b & b \end{matrix}] \cdot [\begin{matrix} 1 - b & b & 0 & 0 \\ 0 & 1 - b & b & 0 \\ 0 & 0 & 1 - b & b \end{matrix}] \cdot (\begin{matrix} b_{0} \\ b_{1} \\ b_{2} \\ b_{3} \end{matrix}) \end{aligned}

$\begin{align} c_0 &= C(a)^3 \cdot C(b)^0 \cdot \begin{pmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \end{pmatrix} = C_1(a)C_2(a)C_3(a)\begin{pmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \end{pmatrix} \\ &= \begin{bmatrix} 1-a & a \end{bmatrix} \cdot \begin{bmatrix} 1-a & a & 0 \\ 0 & 1-a & a \end{bmatrix} \cdot \begin{bmatrix} 1-a & a & 0 & 0 \\ 0 & 1-a & a & 0 \\ 0 & 0 & 1-a &a \end{bmatrix} \cdot \begin{pmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \end{pmatrix} \\ c_1 &= C(a)^2 \cdot C(b)^1 \cdot \begin{pmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \end{pmatrix} = C_1(a)C_2(a)C_3(b)\begin{pmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \end{pmatrix} \\ &= \begin{bmatrix} 1-a & a \end{bmatrix} \cdot \begin{bmatrix} 1-a & a & 0 \\ 0 & 1-a & a \end{bmatrix} \cdot \begin{bmatrix} 1-b & b & 0 & 0 \\ 0 & 1-b & b & 0 \\ 0 & 0 & 1-b & b \end{bmatrix} \cdot \begin{pmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \end{pmatrix} \\ c_2 &= C(a)^1 \cdot C(b)^2 \cdot \begin{pmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \end{pmatrix} = C_1(a)C_2(b)C_3(b)\begin{pmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \end{pmatrix} \\ &= \begin{bmatrix} 1-a & a \end{bmatrix} \cdot \begin{bmatrix} 1-b & b & 0 \\ 0 & 1-b & b \end{bmatrix} \cdot \begin{bmatrix} 1-b & b & 0 & 0 \\ 0 & 1-b & b & 0 \\ 0 & 0 & 1-b & b \end{bmatrix} \cdot \begin{pmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \end{pmatrix} \\ c_3 &= C(b)^1 \cdot C(b)^2 \cdot \begin{pmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \end{pmatrix} = C_1(b)C_2(b)C_3(b)\begin{pmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \end{pmatrix} \\ &= \begin{bmatrix} 1-b & b \end{bmatrix} \cdot \begin{bmatrix} 1-b & b & 0 \\ 0 & 1-b & b \end{bmatrix} \cdot \begin{bmatrix} 1-b & b & 0 & 0 \\ 0 & 1-b & b & 0 \\ 0 & 0 & 1-b & b \end{bmatrix} \cdot \begin{pmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \end{pmatrix} \\ \end{align}$

Efectuând calulele obținem:

pentru $[a,b]=[0,1/3]$ : $c_0\left(-3, 1\right)$ , $c_1\left(-\frac{10}{3}, 2\right)$ , $c_2\left(-\frac{8}{3}, \frac{8}{3}\right)$ , $c_3\left(-\frac{5}{3}, 3\right)$
pentru $[a,b]=[1/3,1]$ : $c_0\left(-\frac{5}{3}, 3\right)$ , $c_1\left(-\frac{1}{3}, \frac{11}{3}\right)$ , $c_2\left(\frac{11}{3}, 3\right)$ , $c_3\left(3,1\right)$

2.Reprezentare grafică:

Lab 4 Info¶

Problemă¶

Scrieți un program (în Python) care să determine un punct de pe o curbă B-spline cubică, în plan, cu 5 puncte de control, folosind algoritmul lui de Boor. Datele de intrare vor fi:

cele 5 puncte de control;

cele 9 noduri;

un $t \in [t_3,t_5]$ ,

iar rezultatul trebuie sa fie $r(t)$ .

Soluție¶

Definim o clasă pentru un punct în plan:

from future import __annotations__

class Point:
    def __init__(self,tuplexy):
        self.x=tuplexy[0]
        self.y=tuplexy[1]

    def __mul__(self, f:float):
        return Point((self.x*f, self.y*f))

    def __rmul__(self, f:float):
        return Point((self.x*f, self.y*f))

    def __add__(self, pt):
        return Point((self.x+pt.x, self.y+pt.y))

    def __str__(self):
        return f"({self.x:.3f}, {self.y:.3f})"

Algoritmul lui de Boor:

def de_boor(d,t,tx):
    n = 3
    N = 2
    r = 0
    for i in range(len(t)-1):
        if t[i] <= tx < t[i + 1]:
            r = i

    assert(n <= r)
    assert (r < n + N)

    dk = [Point((0,0))]*len(d)
    a = [0]*len(d)

    for i in range(r - n, r + 1): dk[i] = d[i]    

    for k in range(1,n+1):
        ndk = [Point((0,0))] * len(d)
        for i in range(r-n+k,r+1):
            a[i] = (tx-t[i])/(t[n+1+i-k]-t[i])
            ndk[i] = (1-a[i])*dk[i-1]+a[i]*dk[i]
        dk = ndk        
    return dk[r]

Datele de intrare (pentru exemplificare):

d = list(map(Point, [(-1, 17), (1, 2), (3, 4), (5, 6), (7, -8)]))
t = [0, 0, 0, 0, 0.5, 1, 1, 1, 1]

Lista punctelor în care vrem să calculăm valoarea curbei:

txs = [0, 0.2, 0.5, 0.7, 0.9]

print("Punctele de control:")
print(", ".join(list(map(str,d))))
print("Nodurile:")
print(t)

Punctele de control:
(-1.000, 17.000), (1.000, 2.000), (3.000, 4.000), (5.000, 6.000), (7.000, -8.000)
Nodurile:
[0, 0, 0, 0, 0.5, 1, 1, 1, 1]

Calculul valorilor $r(t)$ :

print("Valorile curbei in punctele date:")
for tx in txs:
    p = de_boor(d,t,tx)
    print(f"r({str(tx).rjust(3)}) = {p}")

Valorile curbei in punctele date:
r(  0) = (-1.000, 17.000)
r(0.2) = (0.984, 5.656)
r(0.5) = (3.000, 4.000)
r(0.7) = (4.264, 4.240)
r(0.9) = (5.912, -1.280)

Reprezentarea grafică:

pts = [de_boor(d, t, tx * 0.01) for tx in range(int(100*min(t)), int(100*max(t)))]
xs = [p.x for p in pts]
ys = [p.y for p in pts]
plt.plot(xs, ys)
pts = [de_boor(d, t, tx) for tx in txs]
xs = [p.x for p in pts]
ys = [p.y for p in pts]
plt.scatter(xs,ys)
for i in range(len(txs)):
    plt.annotate(f"t={txs[i]}", (xs[i],ys[i]))
plt.show()

Lab 5 Info¶

Problemă¶

Scrieți un program Python care să determine un punct de pe o suprafață Bézier triunghiulară, folosind algoritmul lui de Casteljau. Datele de intrare vor fi:

gradul $n$ al suprafeței;

$(n+1)(n+2)/2$ puncte din spațiu;

un triplet $(u_0,v_0,w_0)$ de numere reale cu suma $1$ ,

iar rezultatul trebuie să fie $r(u_0, v_0,w_0)$ .

Soluție¶

def iterate_idices(n):
    for i in range(n+1):
        for j in range(n+1-i):
            k = n-i-j
            yield (i,j,k)

def ix_add(i,e):
    i0,i1,i2 = i
    e0,e1,e2 = e
    return (i0+e0, i1+e1, i2+e2)

e1 = (1,0,0)
e2 = (0,1,0)
e3 = (0,0,1)

def pt_scale(u, p):
    x,y,z = p
    return (u*x,u*y,u*z)

def pt_add(*args):
    rx,ry,rz = 0,0,0
    for (x,y,z) in args:
        rx += x
        ry += y
        rz += z
    return (rx,ry,rz)

def deCasteljau(n, pts, u,v,w, save_intermediates = False):
    b = pts.copy()

    points = {}

    for r in range(1, n+1):
        new_b = {}
        for ix in iterate_idices(n-r):
            new_b[ix]=pt_add(
                pt_scale(u, b[ix_add(ix, e1)]),
                pt_scale(v, b[ix_add(ix, e2)]),
                pt_scale(w, b[ix_add(ix, e3)]),
            )
        if(save_intermediates): points = {**points, **b}
        b = new_b
    if(save_intermediates): points = {**points, **b}
    return b[(0,0,0)], points

n = 3

# punctele de control
in_pts = [
    (0, 3, 16), (1, 2, 9), (2, 1, 9), (3, 0, 16),
    (0, 2,  4), (1, 1, 5), (2, 0, 4),
    (0, 1,  1), (1, 0, 2),
    (0, 0, 5)
]

# (u, v, w)
coords = [
    (0.4, 0.2, 0.4),
    (0.5, 0.5, 0.0),
    (0.0, 0.8, 0.2),
    (0.7, 0.0, 0.3)
]

pts = {}
k = 0
for ix in iterate_idices(n):
    pts[ix] = in_pts[k]
    k += 1

print(f"n = {n}")
for k in pts.keys(): 
    print(f"b_{k[0]}{k[1]}{k[2]} = {pts[k]}")
for u0,v0,w0 in coords:
    print(f"u0,v0,w0 = {u0}, {v0}, {w0}")
    result, _ = deCasteljau(n, pts, u0, v0, w0)
    print(f"r{u0,v0,w0} = {tuple(map(lambda x:round(x,3), result))}")

n = 3
b_003 = (0, 3, 16)
b_012 = (1, 2, 9)
b_021 = (2, 1, 9)
b_030 = (3, 0, 16)
b_102 = (0, 2, 4)
b_111 = (1, 1, 5)
b_120 = (2, 0, 4)
b_201 = (0, 1, 1)
b_210 = (1, 0, 2)
b_300 = (0, 0, 5)
u0,v0,w0 = 0.4, 0.2, 0.4
r(0.4, 0.2, 0.4) = (0.6, 1.2, 5.072)
u0,v0,w0 = 0.5, 0.5, 0.0
r(0.5, 0.5, 0.0) = (1.5, 0.0, 4.875)
u0,v0,w0 = 0.0, 0.8, 0.2
r(0.0, 0.8, 0.2) = (2.4, 0.6, 12.64)
u0,v0,w0 = 0.7, 0.0, 0.3
r(0.7, 0.0, 0.3) = (0.0, 0.9, 3.344)

wireframe = {}  # se calculeaza niste puncte pe suprafata pentru plot
acc = 20
for i in range(acc + 1):
    for j in range(acc - i + 1):
        p, _ = deCasteljau(n, pts, i / acc, j / acc, (acc - i - j) / acc)
        wireframe[(i, j)] = p

fig = plt.figure()
ax = fig.add_subplot(projection='3d')

xs = [_[0] for _ in wireframe.values()]
ys = [_[1] for _ in wireframe.values()]
zs = [_[2] for _ in wireframe.values()]

surf = ax.plot_trisurf(xs, ys, zs, linewidth=0, zorder=10)
surf.set_facecolor((0,0,1,0.3))
fig.colorbar(surf)

for ix in pts:
    x,y,z = pts[ix]
    ax.scatter(x,y,z,color='b')
    ax.text(x,y,z,f"b{ix[0]}{ix[1]}{ix[2]}",size=10, zorder=1,color='k')

for u0, v0, w0 in coords:
    result, _ = deCasteljau(n, pts, u0, v0, w0)
    x, y, z = result
    ax.scatter(x, y, z, color='r')
    ax.text(x, y, z, f"r({u0},{v0},{w0})", size=10, zorder=1, color='k')

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.view_init(-15,45)
plt.show()