```python
from __future__ import print_function, division

import numpy as np
import math
```

Quadric transforms
===================
*Seth R Johnson*

*January 27, 2016*



```python
def rotate_x(frac):
    theta = 2 * frac * math.pi
    return np.array([[1, 0, 0],
                     [0, math.cos(theta), -math.sin(theta)],
                     [0, math.sin(theta), math.cos(theta)]])

def rotate_y(frac):
    theta = 2 * frac * math.pi
    return np.array([[math.cos(theta), 0, math.sin(theta)],
                     [0, 1, 0],
                     [-math.sin(theta), 0, math.cos(theta)]])

def rotate_z(frac):
    theta = 2 * frac * math.pi
    return np.array([[math.cos(theta), -math.sin(theta), 0.0],
                     [math.sin(theta), math.cos(theta), 0.0],
                     [0, 0, 1]])
```

Transform:
\begin{equation}
\mathbf{x}' \gets \mathbf{Rx} + \mathbf{t}
\end{equation}
becomes
\begin{equation}
\mathbf{\tilde x'}
\gets
\begin{bmatrix}
1 & 0 \\
\mathbf{t} & \mathbf{R} \\
\end{bmatrix}
\mathbf{\tilde x}
\end{equation}
where
\begin{equation}
\mathbf{\tilde x}
= 
\begin{bmatrix}
1 \\
x \\
y \\
z \\
\end{bmatrix}
\end{equation}


```python
def make_transform(rotation, translation):
    result = np.zeros((4,4))
    result[0,0] = 1
    result[1:,0] = translation
    result[1:,1:] = rotation
    return result
```

Inverse transform
\begin{equation}
\mathbf{R}^{-1}(\mathbf{x}' - \mathbf{t}) \to \mathbf{x}
\end{equation}
or
\begin{equation}
\mathbf{\tilde x}
\gets
\begin{bmatrix}
1 & 0 \\
-\mathbf{R}^{-1} \mathbf{t} & \mathbf{R}^{-1} \\
\end{bmatrix}
\mathbf{\tilde x}'
= \mathbf{\tilde R}^{-1} \mathbf{\tilde x}'
\end{equation}


```python
def invert_transform(transform):
    rinv = transform[1:,1:].T
    translation = transform[1:,0]
    result = np.zeros((4,4))
    result[0,0] = 1
    result[1:,0] = -rinv.dot(np.array(translation))
    result[1:,1:] = rinv
    return result
```

The regular quadric equation satisfies
\begin{equation}
\mathbf{\tilde x}^{T}\mathbf{Q}\mathbf{\tilde x} = 0
\end{equation}

Where the quadric equation 
\begin{equation}
ax^2 + by^2 + cz^2 + dxy + eyz + fzx + gx + hy + iz + j = 0
\end{equation}
becomes the symmetric matrix
\begin{equation}
\mathbf{Q} = \begin{bmatrix}
j   & g/2 & h/2 & i/2 \\
g/2 & a   & d/2 & f/2 \\
h/2 & d/2 & b   & e/2 \\
i/2 & f/2 & e/2 & c 
\end{bmatrix}
\end{equation}
It could also be written as the upper-triangular matrix
\begin{equation}
\mathbf{Q} = \begin{bmatrix}
j & g & h & i \\  
0 & a & d & f \\  
0 & 0 & b & e \\  
0 & 0 & 0 & c
\end{bmatrix},
\end{equation}
but this makes it harder to extract a general quadric equation from the result of transformations to $\mathbf{Q}$. Keeping the matrix symmetric makes life easier.


```python
def make_quadric(a,b,c,d,e,f,g,h,i,j):
    return np.array([[j  , g/2, h/2, i/2],
                     [g/2, a  , d/2, f/2],
                     [h/2, d/2, b  , e/2],
                     [i/2, f/2, e/2, c]])
```


```python
def quadric_to_str(q):
    assert q.shape == (4,4)
    assert np.allclose(q.T - q, 0), "Asymmetric: {}".format(q.T - q)
    bits = []
    for (i,j,lab) in ((1,1,'x^2'), (2,2,'y^2'), (3,3,'z^2'),
                      (1,2,'xy'),  (1,3,'xz'),  (2,3,'yz'),
                      (0,1,'x'),   (0,2,'y'),   (0,3,'z'),
                      (0,0,'')):
        val = q[i,j]
        if abs(q[i,j]) > 1e-14:
            if i != j: val *= 2
            bits.append(str(val) + lab)
    return (" + ".join(bits) + " = 0")
```


```python
def make_asymmetric_quadric(a,b,c,d,e,f,g,h,i,j):
    return np.array([[j, g, h, i],
                     [0, a, d, f],
                     [0, 0, b, e],
                     [0, 0, 0, c]])
def asymmetric_quadric_to_str(q):
    assert q.shape == (4,4)
    bits = []
    for (i,j,lab) in ((1,1,'x^2'), (2,2,'y^2'), (3,3,'z^2'),
                      (1,2,'xy'),  (1,3,'xz'),  (2,3,'yz'),
                      (0,1,'x'),   (0,2,'y'),   (0,3,'z'),
                      (0,0,'')):
        val = q[i,j]
        if abs(q[i,j]) > 1e-14:
            bits.append(str(val) + lab)
    return (" + ".join(bits) + " = 0")
```


```python
# Rotate by 45 degrees clockwise around the z axis, then shift right by 1.0; z shouldn't change
rot_and_shift = make_transform(rotate_z(.125), (1, 0, 1))
print(rot_and_shift)
print(rot_and_shift.dot((1,1,0,0)))
```

    [[ 1.          0.          0.          0.        ]
     [ 1.          0.70710678 -0.70710678  0.        ]
     [ 0.          0.70710678  0.70710678  0.        ]
     [ 1.          0.          0.          1.        ]]
    [ 1.          1.70710678  0.70710678  1.        ]



```python
# Transform multiplied by its inverse should be the identity matrix
invert_transform(rot_and_shift).dot(rot_and_shift)
```




    array([[ 1.,  0.,  0.,  0.],
           [ 0.,  1.,  0.,  0.],
           [ 0.,  0.,  1.,  0.],
           [ 0.,  0.,  0.,  1.]])




```python
# Plane at X = 2
plane = make_quadric(0,0,0, 0,0,0, 1,0,0, -2)
print(plane)
print(quadric_to_str(plane))
```

    [[-2.   0.5  0.   0. ]
     [ 0.5  0.   0.   0. ]
     [ 0.   0.   0.   0. ]
     [ 0.   0.   0.   0. ]]
    1.0x + -2.0 = 0



We want to operate on the transformed version
\begin{equation}
\mathbf{\tilde x'}^{T}\mathbf{Q}\mathbf{\tilde x'} = 0
\end{equation}
Substituting the inverse transform $\mathbf{\tilde R}^{-1}$, and expanding the transpose using $(A B)^T = B^T A^T$, this becomes
\begin{equation}
\mathbf{\tilde x}^{T}(\mathbf{\tilde R}^{-1})^T\mathbf{Q}\mathbf{\tilde R}^{-1}\mathbf{\tilde x} = 0
\end{equation}
which is equivalent to a new quadric $\mathbf{Q}'$ with coefficients
\begin{equation}
(\mathbf{\tilde R}^{-1})^T\mathbf{Q}\mathbf{\tilde R}^{-1}
\end{equation}




```python
def transform_quadric(quadric, transform):
    inverse = invert_transform(transform)
    return inverse.T.dot(quadric.dot(inverse))
```


```python
# Expected: plane like (.707)x + (.707)y = 2.707
quadric_to_str(transform_quadric(plane, rot_and_shift))
```




    '0.707106781187x + 0.707106781187y + -2.70710678119 = 0'



### C testing


```python
irs = invert_transform(rot_and_shift)
print(irs)
```

    [[ 1.          0.          0.          0.        ]
     [-0.70710678  0.70710678  0.70710678  0.        ]
     [ 0.70710678 -0.70710678  0.70710678  0.        ]
     [-1.          0.          0.          1.        ]]



```python
print(plane)
```

    [[-2.   0.5  0.   0. ]
     [ 0.5  0.   0.   0. ]
     [ 0.   0.   0.   0. ]
     [ 0.   0.   0.   0. ]]



```python
print(plane.dot(irs))
```

    [[-2.35355339  0.35355339  0.35355339  0.        ]
     [ 0.5         0.          0.          0.        ]
     [ 0.          0.          0.          0.        ]
     [ 0.          0.          0.          0.        ]]



```python
print(irs.T.dot(plane.dot(irs)))
```

    [[-2.70710678  0.35355339  0.35355339  0.        ]
     [ 0.35355339  0.          0.          0.        ]
     [ 0.35355339  0.          0.          0.        ]
     [ 0.          0.          0.          0.        ]]



```python

```

0671e7cbef0bdd3870bba702cb1b589dd9938fac
ipynb
Jupyter Notebook
doc/quadric-transform.ipynb
celeritas-project/orange-port
9aa2d36984a24a02ed6d14688a889d4266f7b1af
Apache-2.0 MIT
doc/quadric-transform.ipynb
celeritas-project/orange-port
9aa2d36984a24a02ed6d14688a889d4266f7b1af
Apache-2.0 MIT
doc/quadric-transform.ipynb
celeritas-project/orange-port
9aa2d36984a24a02ed6d14688a889d4266f7b1af
Apache-2.0 MIT
Qwen/Qwen-72B
1. YES
2. YES