An example from Polygon Mesh Processing page 41
The original equation:
I❤️LA implementation:
`xᵢ` = T_*,1
`xⱼ` = T_*,2
`xₖ` = T_*,3
`n(T)` = (`xⱼ`-`xᵢ`)×(`xₖ`-`xᵢ`)/‖(`xⱼ`-`xᵢ`)×(`xₖ`-`xᵢ`)‖
where
T ∈ ℝ^(3×3)
I❤️LA compiled to C++/Eigen:
/*
`xᵢ` = T_*,1
`xⱼ` = T_*,2
`xₖ` = T_*,3
`n(T)` = (`xⱼ`-`xᵢ`)×(`xₖ`-`xᵢ`)/‖(`xⱼ`-`xᵢ`)×(`xₖ`-`xᵢ`)‖
where
T ∈ ℝ^(3×3)
*/
#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <iostream>
#include <set>
struct pmp_41ResultType {
Eigen::Matrix<double, 3, 1> xᵢ;
Eigen::Matrix<double, 3, 1> xⱼ;
Eigen::Matrix<double, 3, 1> xₖ;
Eigen::Matrix<double, 3, 1> n_left_parenthesis_T_right_parenthesis;
pmp_41ResultType(const Eigen::Matrix<double, 3, 1> & xᵢ,
const Eigen::Matrix<double, 3, 1> & xⱼ,
const Eigen::Matrix<double, 3, 1> & xₖ,
const Eigen::Matrix<double, 3, 1> & n_left_parenthesis_T_right_parenthesis)
: xᵢ(xᵢ),
xⱼ(xⱼ),
xₖ(xₖ),
n_left_parenthesis_T_right_parenthesis(n_left_parenthesis_T_right_parenthesis)
{}
};
pmp_41ResultType pmp_41(const Eigen::Matrix<double, 3, 3> & T)
{
Eigen::Matrix<double, 3, 1> xᵢ = T.col(1-1);
Eigen::Matrix<double, 3, 1> xⱼ = T.col(2-1);
Eigen::Matrix<double, 3, 1> xₖ = T.col(3-1);
Eigen::Matrix<double, 3, 1> n_left_parenthesis_T_right_parenthesis = ((xⱼ - xᵢ)).cross((xₖ - xᵢ)) / double((((xⱼ - xᵢ)).cross((xₖ - xᵢ))).lpNorm<2>());
return pmp_41ResultType(xᵢ, xⱼ, xₖ, n_left_parenthesis_T_right_parenthesis);
}
void generateRandomData(Eigen::Matrix<double, 3, 3> & T)
{
T = Eigen::MatrixXd::Random(3, 3);
}
int main(int argc, char *argv[])
{
srand((int)time(NULL));
Eigen::Matrix<double, 3, 3> T;
generateRandomData(T);
pmp_41ResultType func_value = pmp_41(T);
std::cout<<"return value:\n"<<func_value.n_left_parenthesis_T_right_parenthesis<<std::endl;
return 0;
}
I❤️LA compiled to Python/NumPy/SciPy:
"""
`xᵢ` = T_*,1
`xⱼ` = T_*,2
`xₖ` = T_*,3
`n(T)` = (`xⱼ`-`xᵢ`)×(`xₖ`-`xᵢ`)/‖(`xⱼ`-`xᵢ`)×(`xₖ`-`xᵢ`)‖
where
T ∈ ℝ^(3×3)
"""
import numpy as np
import scipy
import scipy.linalg
from scipy import sparse
from scipy.integrate import quad
from scipy.optimize import minimize
class pmp_41ResultType:
def __init__( self, xᵢ, xⱼ, xₖ, n_left_parenthesis_T_right_parenthesis):
self.xᵢ = xᵢ
self.xⱼ = xⱼ
self.xₖ = xₖ
self.n_left_parenthesis_T_right_parenthesis = n_left_parenthesis_T_right_parenthesis
def pmp_41(T):
T = np.asarray(T, dtype=np.float64)
assert T.shape == (3, 3)
xᵢ = T[:, 1-1]
xⱼ = T[:, 2-1]
xₖ = T[:, 3-1]
n_left_parenthesis_T_right_parenthesis = np.cross((xⱼ - xᵢ), (xₖ - xᵢ)) / np.linalg.norm(np.cross((xⱼ - xᵢ), (xₖ - xᵢ)), 2)
return pmp_41ResultType(xᵢ, xⱼ, xₖ, n_left_parenthesis_T_right_parenthesis)
def generateRandomData():
T = np.random.randn(3, 3)
return T
if __name__ == '__main__':
T = generateRandomData()
print("T:", T)
func_value = pmp_41(T)
print("return value: ", func_value.n_left_parenthesis_T_right_parenthesis)
I❤️LA compiled to MATLAB:
function output = pmp_41(T)
% output = pmp_41(T)
%
% `xᵢ` = T_*,1
% `xⱼ` = T_*,2
% `xₖ` = T_*,3
% `n(T)` = (`xⱼ`-`xᵢ`)×(`xₖ`-`xᵢ`)/‖(`xⱼ`-`xᵢ`)×(`xₖ`-`xᵢ`)‖
%
% where
%
% T ∈ ℝ^(3×3)
if nargin==0
warning('generating random input data');
[T] = generateRandomData();
end
function [T] = generateRandomData()
T = randn(3, 3);
end
assert( isequal(size(T), [3, 3]) );
x_i = T(:, 1);
x_j = T(:, 2);
x_k = T(:, 3);
n_T = cross((x_j - x_i), (x_k - x_i)) / norm(cross((x_j - x_i), (x_k - x_i)), 2);
output.x_i = x_i;
output.x_j = x_j;
output.x_k = x_k;
output.n_T = n_T;
end
I❤️LA compiled to LaTeX:
\documentclass[12pt]{article}
\usepackage{mathdots}
\usepackage[bb=boondox]{mathalfa}
\usepackage{mathtools}
\usepackage{amssymb}
\usepackage{libertine}
\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator*{\argmin}{arg\,min}
\usepackage[paperheight=8in,paperwidth=4in,margin=.3in,heightrounded]{geometry}
\let\originalleft\left
\let\originalright\right
\renewcommand{\left}{\mathopen{}\mathclose\bgroup\originalleft}
\renewcommand{\right}{\aftergroup\egroup\originalright}
\begin{document}
\begin{center}
\resizebox{\textwidth}{!}
{
\begin{minipage}[c]{\textwidth}
\begin{align*}
\textit{x\textsubscript{i}} & = \mathit{T}_{*, 1} \\
\textit{x\textsubscript{j}} & = \mathit{T}_{*, 2} \\
\textit{xₖ} & = \mathit{T}_{*, 3} \\
\textit{n(T)} & = \frac{\left( \textit{x\textsubscript{j}} - \textit{x\textsubscript{i}} \right) × \left( \textit{xₖ} - \textit{x\textsubscript{i}} \right)}{\left\|\left( \textit{x\textsubscript{j}} - \textit{x\textsubscript{i}} \right) × \left( \textit{xₖ} - \textit{x\textsubscript{i}} \right)\right\|_2} \\
\intertext{where}
\mathit{T} & \in \mathbb{R}^{ 3 \times 3 } \\
\\
\end{align*}
\end{minipage}
}
\end{center}
\end{document}
I❤️LA LaTeX output: