Polygon Mesh Processing page 41

An example from Polygon Mesh Processing page 41

The original equation:

placeholder

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: