Computer Graphics with OpenGL

OpenGL Loaders: [ESSENTIAL]

• Libraries that abstracts OpenGL function pointers loading in a platform-independent way.
• GLEW - OpenGL Extension Wrangler [MOST USED]
  • "The OpenGL Extension Wrangler Library (GLEW) is a cross-platform open-source C/C++ extension loading library. GLEW provides efficient run-time mechanisms for determining which OpenGL extensions are supported on the target platform. OpenGL core and extension functionality is exposed in a single header file. GLEW has been tested on a variety of operating systems, including Windows, Linux, Mac OS X, FreeBSD, Irix, and Solaris."
• GLAD - [MOST-USED] Multi-Language GL/GLES/EGL/GLX/WGL Loader-Generator based on the official specs.
• GitHub - cginternals/glbinding
  • "A C++ binding for the OpenGL API, generated using the gl.xml specification."
• GitHub - anholt/libepoxy
  • "Epoxy is a library for handling OpenGL function pointer management for you."
• GitHub - imakris/glatter
  • "An OpenGL loading library, with support for GL, GLES, EGL, GLX and WGL"
• Galogen OpenGL Loader Generator
  • "Galogen is an OpenGL loader generator. Given an API version and a list of extensions, Galogen will produce corresponding headers and code that load the exact OpenGL entry points you need. The produced code can then be used directly by your C or C++ application, without having to link against any additional libraries."
• GitHub - SFML/SFML-glLoadGen
  • Customized glLoadGen for SFML

Window System Abstraction

Libraries for window systems, event handling and OpenGL context abstraction: [ESSENTIAL]

• Abstract platform-specific window system and event handling.
• GLFW [BEST] [MOST-USED]
  • C library that provides graphics windows for OpenGL, Vulkan, OpenGL ES and deals with event handling.
• SDL (Simple Direct Media Layer) [BEST] [MOST-USED]
  • Cross-platform C library that provides windows and event handling for many computer graphics APIs such as OpenGL, Vulkan and DirectX. SLD also has facilities for dealing with audio, joystick, CD-ROM, network and threads.
• SFML (Simple and Fast Multimedia Library) [MOST-USED]
  • "SFML provides a simple interface to the various components of your PC, to ease the development of games and multimedia applications. It is composed of five modules: system, window, graphics, audio and network."
• GLUT (FreeGlut) - OpenGL Utility Toolkit
  • Deals with window creation, OpenGL initialization, event handling and so on.
  • Docs: https://www.glfw.org/documentation.html

Graphics Math Libraries

OpenGL Math and Computer Graphics Math: [ESSENTIAL]

• GLM (OpenGL Mathematics Library) [MOST-USED]
  • Source code: https://github.com/g-truc/glm
  • Header-only C++ library that provides classes for computer graphics mathematics such as: 2D, 3D and homogeneous coordinate vector; 2D, 3D and homogeneous coordinate transformation matrices; quaternions and subroutines for computing camera, perspective or orthogonal transformation matrices.
• GitHub - Kazade/kazmath - A C math library targeted at games
  • "Kazmath is a simple 3D maths library written in C. It was initially coded for use in my book, Beginning OpenGL Game Programming - Second edition, but rapidly gained a life of its own. Kazmath is now used by many different projects, and apparently is used in 25% of the worlds mobile games (yeah, I don't believe it either - but it's used in Cocos2d-x)."
• GitHub - recp/cglm - Highly Optimized Graphics Math (glm) for C
• See: GitHub - chunkyguy/Math-Library-Test - A comparison of the various major math libraries for speed and ease of use.

Image Loading for texture

• GLI - OpenGL Image
  • "OpenGL Image (GLI) is a header only C++ image library for graphics software. GLI provides classes and functions to load image files (KTX and DDS), facilitate graphics APIs texture creation, compare textures, access texture texels, sample textures, convert textures, generate mipmaps, etc."
• smoked-herring/sail - Squirrel Abstract Image Library
  • "SAIL is a format-agnostic cross-platform image decoding library providing rich APIs, from one-liners to complex use cases with custom I/O sources. It enables a client to read and write static, animated, multi-paged images along with their meta data and ICC profiles."
  • Note: Supports APNG, BMP, GIF, JPEG, PNG and TIFF image formats.
• stb_image.h - single-file header-only C library for loading images from several file formats including, jpeg, bmp, tga, ppm, pgm, gif and so on.
• SOIL - Simple OpenGL Image Library
  • "SOIL is a tiny C library used primarily for uploading textures into OpenGL. It is based on stb_image version 1.16, the public domain code from Sean Barrett (found here). It has been extended to load TGA and DDS files, and to perform common functions needed in loading OpenGL textures. SOIL can also be used to save and load images in a variety of formats."
• freeimage
  • "FreeImage is an Open Source library project for developers who would like to support popular graphics image formats like PNG, BMP, JPEG, TIFF and others as needed by today's multimedia applications. FreeImage is easy to use, fast, multithreading safe, compatible with all 32-bit or 64-bit versions of Windows, and cross-platform (works both with Linux and Mac OS X)."
• gldraw
  • "With glraw you can preconvert your texture assets and load them without the need of any image library. The generated raw files can easily be read. For this, glraw also provides a minimal Raw-File reader that you can either source-copy or integrate as C++ library into your project. Image to OpenGL texture conversion can be done either by glraws command line interface, e.g., within an existing tool-chain, or at run-time with glraw linked as asset library (requires linking Qt)."
• lodepng
  • "LodePNG is a PNG image decoder and encoder, all in one, no dependency or linkage to zlib or libpng required. It's made for C (ISO C90), and has a C++ wrapper with a more convenient interface on top."
• libpng
  • "libpng is the official PNG reference library. It supports almost all PNG features, is extensible, and has been extensively tested for over 23 years. The home site for development versions (i.e., may be buggy or subject to change or include experimental features) is https://libpng.sourceforge.io/, and the place to go for questions about the library is the png-mng-implement mailing list. libpng is available as ANSI C (C89) source code and requires zlib 1.0.4 or later (1.2.5 or later recommended for performance and security reasons). The current public release, libpng 1.6.37, fixes the use-after-free security vulnerability noted below, as well as an ARM NEON memory leak in the palette-to-RGB(A) expansion code (png_do_expand_palette())."
• libspng
  • "libspng (simple png) is a C library for reading and writing Portable Network Graphics (PNG) format files with a focus on security and ease of use. It is licensed under the BSD 2-clause “Simplified” License."
• IJG - Indepedent JPEG Group
  • "IJG is an informal group that writes and distributes a widely used free library for JPEG image compression. The first version was released on 7-Oct-1991. The current version is release 9d of 12-Jan-2020. This is a stable and solid foundation for many application's JPEG support. You can find our original code and some supporting documentation in the directory files. There is a Windows format package in zip archive format jpegsr9d.zip and a Unix format package in tar.gz archive format jpegsrc.v9d.tar.gz. A collection of modified versions with adaptions and error fixes for system maintenance is available on jpegclub.org in the directory support."

Asset/Object Loaders - Mesh Importing Libraries [ESSENTIAL]

Object loader / Asset import libraries:

• Assimp - The Open-Asset-Import Library [MOST-USED]
  • "The Open Asset Import Library (short name: Assimp) is a portable Open-Source library to import various well-known 3D model formats in a uniform manner. The most recent version also knows how to export 3d files and is therefore suitable as a general-purpose 3D model converter. See the feature-list."
  • Note: Allows importing blender-generated models/assets in OpenGL.
  • Repository: https://github.com/assimp/assimp
  • See: https://www.khronos.org/opengl/wiki/Tools/Open_Asset_Import
• OBJ-Loader
  • "OBJ Loader is a simple, header only, .obj model file loader that will take in a path to a file, load it into the Loader class object, then allow you to get the data from each mesh loaded. This will load each mesh within the model with the corresponding data such as vertices, indices, and material. Plus a large array of vertices, indices and materials which you can do whatever you want with."
• tinyobjloader/tinyobjloader
  • "Tiny but powerful single file wavefront obj loader written in C++03. No dependency except for C++ STL. It can parse over 10M polygons with moderate memory and time. tinyobjloader is good for embedding .obj loader to your (global illumination) renderer ;-)."
• codelibs/libdxfrw
  • "C++ library to read and write DXF/DWG (Autocad file format) files. - libdxfrw is a free C++ library to read and write DXF files in both formats, ascii and binary form. Also can read DWG files from R14 to the last V2015. It is licensed under the terms of the GNU General Public License version 2 (or at you option any later version)."

Text and Font Rendering

• The FreeType Project
  • Brief: "It is written in C, designed to be small, efficient, highly customizable, and portable while capable of producing high-quality output (glyph images) of most vector and bitmap font formats."
• OGLFT: OpenGL-FreeType Library
  • Brief: "This C++ library supplies an interface between the fonts on your system and an OpenGL or Mesa application. It uses the excellent FreeType library to read font faces from their files and renders text strings as OpenGL primitives."
• GitHub - vallentin/glText - [HEADER-ONLY-LIBRARY]
  • Brief: "glText is a simple cross-platform single header text rendering library for OpenGL. glText requires no additional files (such as fonts or textures) for drawing text, everything comes pre-packed in the header."
• GitHub - MartinPerry/OpenGL-Font-Rendering
  • Brief: "Rendering UNICODE fonts with OpenGL This library is still work-in-progress. This is a working beta version."
• libdrawtext - OpenGL text rendering library
  • Brief: "Libdrawtext uses freetype2 for glyph rasterization. If you would rather avoid having freetype2 as a dependency, you can optionally compile libdrawtext without it, and use pre-rendered glyphmaps. Glyphmaps can be generated by the included font2glyphmap tool, or by calling dtx_save_glyphmap."
• GitHub - codetiger/Font23D - Convert any text to a 3d mesh using any font style
  • Brief: "Font23D is a C++ library for creating a 3d mesh of any Text in the given True type font."

C++ Wrappers

• OGplus - Self described as C++ Wrapper for modern OpengL.
  • Brief: "OGLplus is a header-only library which implements a thin object-oriented facade over the OpenGL® (version 3 and higher) C-language API. It provides wrappers which automate resource and object management and make the use of OpenGL in C++ safer and easier."
  • Repository: https://github.com/matus-chochlik/oglplu2
• Oglwrap
  • Brief: "Oglwrap is a lightweight, cross-platform, object-oriented, header-only C++ wrapper for modern (2.1+) OpenGL, that focuses on preventing most of the trivial OpenGL errors, and giving as much debug information about the other errors, as possible."
• Globjects - globjects is a cross-platform C++ wrapper for OpenGL API objects
  • Brief: "globjects provides object-oriented interfaces to the OpenGL API (3.0 and higher). It reduces the amount of OpenGL code required for rendering and facilitates coherent OpenGL use by means of an additional abstraction layer to glbinding and GLM. Common rendering tasks and processes are automated and missing features of specific OpenGL drivers are partially simulated or even emulated at run-time."

Form-based Graphical User Interface

• Dear imgui
  • Brief: "Dear ImGui is a bloat-free graphical user interface library for C++. It outputs optimized vertex buffers that you can render anytime in your 3D-pipeline enabled application. It is fast, portable, renderer agnostic and self-contained (no external dependencies)."
  • See: An introduction to Dear Imgui Library

Non-categorized / Miscellaneous

• bkaradzic/bgfx - "Bring Your Own Engine/Framework"
  • Brief: Library for providing a unified interface to many GPU accelerated graphics APIs, including OpenGL, OpenGL ES, WebGL, WebGPU (W3C consoritum), DirectX (Microsoft), Metal (Apple) and Vulkan.
  • Note: This library also abstracts the shading language, a subset of GLSL OpenGL, that works with all mentioned GPU accelerated rendering backends.
  • Documentation:
    • https://bkaradzic.github.io/bgfx/bgfx.html
    • https://bkaradzic.github.io/bgfx/examples.html
  • Community:
    • https://gitter.im/bkaradzic/bgfx
  • See:
    • Shader Debugging for BGFX Rendering Engine (Kai Wang)
• GLSDK - Unofficial OpenGL Software Development Kit
  • Brief: "The Unofficial OpenGL Software Development Kit is a collection of libraries and utilities that will help you get started working with OpenGL. It provides a unified, cross-platform build system to make compiling the disparate libraries easier. Many of the components of the SDK are C++ libraries. Each component of the SDK specifies the terms under which they are distributed. All licenses used by components are approximately like the MIT license in permissivity. The parts of the SDK responsible for maintaining the build, as well as all examples, are distributed under the MIT License."
• NXPmicro/gtec-demo-framework
  • "A multi-platform framework for fast and easy demo development. The framework abstracts away all the boilerplate & OS specific code of allocating windows, creating the context, texture loading, shader compilation, render loop, animation ticks, benchmarking graph overlays etc. Thereby allowing the demo/benchmark developer to focus on writing the actual 'demo' code. Therefore demos can be developed on PC or Android where the tool chain and debug facilities often allows for faster turnaround time and then compiled and deployed without code changes for other supported platforms. The framework also allows for ‘real’ comparative benchmarks between the different OS and windowing systems, since the exact same demo/benchmark code run on all of them."
  • Supported Operating Systems: Android NDK; Linux with various windowing systems (Yocto); Ubuntu 18.04; Windows 7+

WebGL Companion Libraries: (OpenGL on the WEB)

• glMatrix - GLM math library ported to Javascript. Computer graphics math library, similar to OpenGL GLM math library.
• Math.GL - "math.gl is JavaScript math library focused on geospatial and 3D use cases, designed as a composable, modular toolbox. math.gl provides a core module with classic vector and matrix classes, and a suite of optional modules implementing various aspects of geospatial and 3D math. While the math.gl is highly optimized for use with the WebGL and WebGPU APIs, math.gl itself has no WebGL dependencies."
• GLM-JS - "glm-js is an experimental JavaScript implementation of the OpenGL Mathematics (GLM) C++ Library."
• Three.JS - High level wrapper library around WebGL, that provides many high level features that includes: camera objects; scene graphs; geometry; perspective; forward kinematics, inverse kinematics; graphics math library containing, matrices, vectors, quaternions; image loaders; animation and more.
• phoria.js - [NOT WEBGL] "JavaScript library for simple 3D graphics and visualisation on a HTML5 canvas 2D renderer. It does not use WebGL. Works on all HTML5 browsers, including desktop, iOS and Android."
  • Note: It does use WebGL, but it is an interesting codebase about computer graphics algorithms implementation.

• https://free3d.com/3d-models/obj-file
• common-3d-test-models
  • "This is a repository containing common 3D test models in original format with original source if known. While a similar list exists on wikipedia, it does not host the actual models and is incomplete."
• https://github.com/pichiliani/ModelsOBJ
  • "ModelsOBJ - .OBJ files with free 3D models (CC license or royalty free)"
• Gist - cube.obj
  • Gist containing a cube 3D model obj file with vertices and normals.
• Teapot.obj - Teapot common 3D test model with normals.

Note: 3D models can be created using applications such as Blender, MeshlLab, Autodesk Maya, Solidworks and so on.

See also:

• List of common 3D test models
• Smoothing .obj 3D models (Guru Mulay)

Obj File Documentation:

• https://docs.fileformat.com/3d/obj/
• Weverfron OBJ File Format - Library of Congress

• OpenGL Loading Library - OpenGL Wiki / Khronos Group
  • Brief: "An OpenGL Loading Library is a library that loads pointers to OpenGL functions at runtime, core as well as extensions. This is required to access functions from OpenGL versions above 1.1 on most platforms. Extension loading libraries also abstracts away the difference between the loading mechanisms on different platforms."
• Load OpenGL Functions - OpenGL Wiki / Khronos Group
  • Brief: "Loading OpenGL Functions is an important task for initializing OpenGL after creating an OpenGL context. You are strongly advised to use an OpenGL Loading Library instead of a manual process. However, if you want to know how it works manually, read on."
• Image Libraries - OpenGL Wiki / Khronos Group
• When do I need to use an OpenGL function loader? - Stack Overflow
• KeyJ's Blog : Blog Archive » Modern OpenGL with lcc-win32, the hard way
• Loading OpenGL without GLEW
• Using OpenGL With SDL - LibSDL
• Tutorial1: Creating a Cross Platform OpenGL 3.2 Context in SDL (C / SDL) - OpenGL Wiki

Command Line Shortcuts for Troubleshooting

The following applications can be accessed either from command line (cmd.exe shell) or via the shortcut Windows Key + R.

• $ dxdiag => Tool for DirectX and OpenGL diagnostics.
• $ devmgmt.msc => Device manager shortcut.
• $ msinfo32 => View details about hardware.
• $ systeminfo => View summarized information about operating system, hardware, RAM memory, disk space and patches. Note: It only works from command line.

GPU details and OpenGL Implementation Information

• OpenGL Extensions Viewer 6 | realtech VR
  • Graphical application that allows viewing details about Vulkan and OpenGL implementation at the current machine.
  • Brief: "A reliable software which displays useful information about the current OpenGL 3D accelerator and new Vulkan 3D API. This program displays the vendor name, the version implemented, the renderer name and the extensions of the current OpenGL 3D accelerator."
• GPU Caps View - Geeks3D
  • Brief: "A new version of GPU Caps Viewer is available. GPU Caps Viewer is a graphics card / GPU information and monitoring utility that quickly describes the essential capabilities of your GPU including GPU type, amount of VRAM , OpenGL, Vulkan, OpenCL and CUDA API support level."
• GPU Test
  • Brief: "GpuTest is a cross-platform (Windows, Linux and Max OS X) GPU stress test and OpenGL benchmark. GpuTest comes with several GPU tests including some popular ones from Windows'world (FurMark or TessMark)."

Setting Default GPU on Windows

• How to update Video Card Driver and Set Default GPU

Install OpenGL development dependencies

Install OpenGL development dependencies on Ubuntu or Debian-like distributions:

$ sudo apt-get install -y libgl1-mesa-dev libglu1-mesa-dev freeglut3-dev  libgl1-mesa-dev libxrandr-dev
$ sudo apt-get install -y libxinerama-dev libxcursor-dev libxi-dev

Information about OpenGL and graphics card

Vendor:

 $ >> glxinfo | grep -E "vendor"
server glx vendor string: SGI
client glx vendor string: Mesa Project and SGI
OpenGL vendor string: Intel

Rendering:

 $ >> glxinfo | grep -E "rendering"
direct rendering: Yes

Version:

 $ >> glxinfo | grep -E "version"
server glx version string: 1.4
client glx version string: 1.4
GLX version: 1.4
    Max core profile version: 4.6
    Max compat profile version: 4.6
    Max GLES1 profile version: 1.1
    Max GLES[23] profile version: 3.2
OpenGL core profile version string: 4.6 (Core Profile) Mesa 20.3.2
OpenGL core profile shading language version string: 4.60
OpenGL version string: 4.6 (Compatibility Profile) Mesa 20.3.2
OpenGL shading language version string: 4.60
OpenGL ES profile version string: OpenGL ES 3.2 Mesa 20.3.2
OpenGL ES profile shading language version string: OpenGL ES GLSL ES 3.20
    GL_EXT_shader_implicit_conversions, GL_EXT_shader_integer_mix,

Display:

 $ >> glxinfo | grep -E "display"
name of display: :1
display: :1  screen: 0

Computer graphics math is based on vector algebra, linear algebra and affine transforms. Those concepts are not exclusive to OpenGL, they are essential and universal to all computer graphics APIS - Application Programming Interfaces, including OpenGL, DirectX, Metal, Vulkan, WebGL and html5 canvas.

Given 2 3D vectors \(\vec{A} = [ x_A, y_A, z_A ]\) and \(\vec{B} = [ x_B, y_B, z_B ]\) , the following properties can be defined:

Vector Sum

Vector Difference / Subtraction

Vector Norm

The norm or magnitude of A, \(|| \vec{A} ||\) is given by:

Normalized Vector / Unit Vector

A normalized vector, is a vector with the same direction than a given vector, but with norm equal to one. A normalized vector of A can be computed as:

Where:

• \(\hat{i}\), \(\hat{j}\), \(\hat{k}\) are the unit vectors of axis X, Y and Z respectively.

Distance between two position (point) vectors

Given two vectors A and B which represent the position relative to the coordinate system origin or point. (0, 0, 0). The distance between A and B, or the length of the vector difference between A and B is determined by the following relation.

Element-wise Product

The element-wise product between two vectors is a vector which the components are the product between the two vectors components. As there is no standard mathematical notation for element-wise product between two vectors, the symbol \(\odot\) will be used for denoting this operation and avoiding confusion with dot product or vector product. This notation is useful for describing illumination (lighting) calculations, which most books and introduction materials present without a clear notation.

Dot Product (a.k.a - Scalar Product)

The vector dot product is:

• Where:
  • \(|| \vec{A} ||\) is the norm of vector A
  • \(|| \vec{B} ||\) is the norm of vector B
• The vectors \(\vec{A}\) and \(\vec{B}\) are orthogonal, the angle between them is 90 degrees, when the dot product is zero.

The dot product has the property:

And also:

Where:

• \(\theta\) is the angle between vectors A and B.

If the vectors \(\vec{A}\) and \(\vec{B}\) are expressed as column matrices, the product can be computed as a matrix multiplication.

• In the following equations. \(A^T\) is the transpose of matrix A.

Cross Product (a.k.a - Vector product)

The cross product between two vectors \(\vec{A}\) and \(\vec{B}\) results in a vector which is perpendicular (orthogonal) to both A and B.

Where \([A]_{\times}\) is the cross product matrix.

The cross product have the following properties:

• Where:
  • \(\vec{n}\) is unit vector with the same direction as the cross product vector.
  • \(\theta\) is the angle between the two underlying vectors.

Vector Triple Product

Relation between vectors

Orthogonal (perpendicular) vectors:

Two vectors \(\vec{A}\) and \(\vec{B}\) are orthogonal, the angle between them is 90 degrees or PI/2 radians, if their dot product is zero.

Parallel vectors

Two vectors \(\vec{A}\) and \(\vec{B}\) are parallel when the angle between them are zero or 180 degrees. When this happens, then:

Same direction

Let the unit vector that points in the same direction of a vector \(\vec{A}\) be \([\vec{A}]_u\), where \([ \cdot ]_u\) is a operator that determines the unit vector of a vector.

A pair of vector \(\vec{A}\) and \(\vec{B}\) have the same direction both are parallel and the angle between them are zero. When this happens, both have the same unit vector or then angle between the are zero:

Or:

Two vectors are parallel and have opposite direction (angle between them is 180 degrees) if the following relation holds.

Angle between two vectors

The angle between two 3D vectors can be determined by using their dot product and cross product.

Affine transforms, which are represented by matrices, are a particular class of linear transforms which preserves ratios between distances, colinearity and parallelism. Affine transforms has many applications that includes, computer graphics, computer vision, image processing, CAD (Computer Aided Design) and robotics.

Outline of properties preserved by affine transforms:

• Ratios between distances.
• Colinearity
  • => Points in the same line, remains in the same line, after the transform was applied.
• Parallelism
  • => Lines that are parallels, remains parallel.

Types of geometric linear transforms:

• Affine Transforms
  • => Preserves ratios, colinearity and parallelism. Some affine transforms are: identity, translation, rotation and scaling.
• Projection Transforms
  • => They not preserve parallelism. However, they preserve colinearity. Affine transforms are a particular case of projection transforms.
  • => Some projection transforms are:
    • Orthogonal view transform
    • Projection view transform
• Rigid body transforms
  • => Are a particular case of affine transforms. The set of rigid body transforms comprises: identity, translation and rotation. This set of transforms does not include shear and scaling.
  • => Use cases: Computer graphics, Newtonian mechanics, robotics, aerospace and many other cases.

Some affine transforms are:

• Identity (Identity matrix) => Represents no transform, all points or vertices remain the same.
• Translation
• Scaling
• Reflection
• Rotation
• Shear

Further Reading

General:

• CMU - Transformations
• Cornell CS4620 - 2D Geometric Transformations
• AML170 CAD - 3D Transformations - Lecture 6
• CS_CS3100 (Introduction to) Computer Graphics

Coordinate Systems:

• Coordinate system - Wikipedia
  • Brief: "n geometry, a coordinate system is a system that uses one or more numbers, or coordinates, to uniquely determine the position of the points or other geometric elements on a manifold such as Euclidean space."
• Homogeneous coordinates - Wikipedia
  • Brief: "Homogeneous coordinates have a range of applications, including computer graphics and 3D computer vision, where they allow affine transformations and, in general, projective transformations to be easily represented by a matrix."

Rotation Matrix and right-hand-rule:

• Right-hand rule - Wikipedia
  • Brief: "In mathematics and physics, the right-hand rule is a common mnemonic for understanding orientation of axes in three-dimensional space."
• Rotation matrix - Wikipedia
  • Brief: "In linear algebra, a rotation matrix is a transformation matrix that is used to perform a rotation in Euclidean space."
• Rotation formalisms in three dimensions - Wikipedia
  • Brief: "In geometry, various formalisms exist to express a rotation in three dimensions as a mathematical transformation. In physics, this concept is applied to classical mechanics where rotational (or angular) kinematics is the science of quantitative description of a purely rotational motion."
• Axes conventions - Wikipedia
  • Brief: "In ballistics and flight dynamics, axes conventions are standardized ways of establishing the location and orientation of coordinate axes for use as a frame of reference."
• Euler angles - Wikipedia
  • Brief: "The Euler angles are three angles introduced by Leonhard Euler to describe the orientation of a rigid body with respect to a fixed coordinate system."
• Yaw, Pitch, Roll angles - Aircraft principal axes - Wikipedia
  • Brief: Rotation angles convention for aircrafts.
• Celestial coordinate system - Wikipedia
  • Brief: "In astronomy, a celestial coordinate system (or celestial reference system) is a system for specifying positions of satellites, planets, stars, galaxies, and other celestial objects relative to physical reference points available to a situated observer (e.g. the true horizon and north cardinal direction to an observer situated on the Earth's surface)"
• Conversion between quaternions and Euler angles - Wikipedia
  • Brief: "Spatial rotations in three dimensions can be parametrized using both Euler angles and unit quaternions. This article explains how to convert between the two representations. Actually this simple use of "quaternions" was first presented by Euler some seventy years earlier than Hamilton to solve the problem of magic squares."

Affine Tranform Matrices in SVG, Html5 Canvas and WebGL:

• Taming 2D Transforms in Html5 Canvas
• Html5 Canvas: Matrix Transforms
• Html5 Rocks - WebGL Transforms (WebGL - OpenGL for the web)
• MDN - SVG Transform matrices (Mozilla)
• W3C - 7 Coordinate Systems, Transformations and Units (SVG)
• SVG Transformations with affine matrices (Paul Cowan)
• Examples using matrix transforms in SVG

Affine Transform Matrices in DirectX (Direct3D):

• MSDN - Transforms (Direct3D 9)
• MSDN - World Transform
• MSDN - View Transform
• MSDN - Projection Transform
• MSDN - Viewport and Clipping
• D3DXMatrixAffineTransformation function (D3dx9math.h)

Affine Transform Matrices in Java AWT:

• Oracle - Class AffineTransform
• Java 2D Graphics, Simple Affine Transforms

Affine Transform Matrices in non-categorized APIs

• BRL Cad - Object Motion and Matrix Manipulation
• Affine Transformation Tutorial - GeoTools OSGeo
• FreeCAD - Macro MatrixTransform
• A complete overview of matrix transformations in the SOLIDWORKS API (SolidWorks CAD for mechanical engineering)

General form of 2D affine transforms

2D affine transforms, including translation, rotation and scaling can be represented by matrices with the following format, that transforms homogeneous coordinates from one coordinate system to another. Homogeneous coordinates, are 2D or 3D coordinates with an extra pseudo-coordinate, often designated by 'w', for representing translations affine transforms in the same way as rotations and scaling.

This affine transform performs the coordinate transformation from the coordinate frame C2 (which axis are x, y), which could an object local-space coordinate, to coordinate frame C1 (which axis are x', y') which could be a world-coordinate system. The matrix transforms homogenous coordinates, which are coordinates with an extra parameter w = 1 for allowing translation transformation to be expressed in the same way as rotation transformations.

The multiplication between two affine transforms also results in a affine transform. Consider two affine transforms \(A = a_{ij}\) and \(B = b_{ij}\). The product between these two affine transforms is also an affine transform.

Where:

• \(c_{11} = a_{11} \cdot b_{11} + a_{12} \cdot b_{21}\)
• \(c_{12} = a_{11} \cdot b_{12} + a_{12} \cdot b_{22}\)
• \(c_{21} = a_{21} \cdot b_{11} + a_{22} \cdot b_{21}\)
• \(c_{22} = a_{21} \cdot b_{12} + a_{22} \cdot b_{22}\)
• \(c_{13} = a_{11} \cdot b_{13} + a_{12} \cdot b_{23} + a_{13}\)
• \(c_{23} = a_{21} \cdot b_{13} + a_{22} \cdot b_{23} + a_{23}\)

2D Canonical Affine Transforms

Identity matrix

• Causes no coordinate transformation. The result coordinate system and vertices remain at the same position. The identity matrix is also a reasonable initial default value for the model matrix, view matrix and projection matrix.

Translation:

• Translate the coordinate system from one position into another and translate vertices.

Scaling:

• Increases or decreases object size. This transformation allows resizing an object without sending the vertices multiple times to the GPU in the case of a GPU accelerated graphics API.
• Where: \(s_x\) is the scale for X axis and \(s_y\) is the scale for the y axis.

Shear:

Rotation around Z axis:

• Where \(\theta\) is the angle of counterclockwise rotation around Z axis. A negative angle results in a rotation in the opposite direction.

The window-to-viewport affine transform maps coordinates from a world-space window to a viewport (physical display coordinates). The world-space window defines what region of the world-space will be viewed. It is specified using the world-space coordinates \(x_{wmin}\), \(x_{wmax}\), \(y_{wmin}\) and \(y_{wmax}\). The viewport is a rectangle within the graphics display window to where the world-window coordinates will be mapped to. The viewport window is defined by the coordinates: \(x_{vmin}\), \(x_{vmax}\), \(y_{vmin}\), \(y_{vmax}\) - given in screen device-coordinates with origin at the graphics display bottom left corner.

Note: In OpenGL this transform can be obtained using the subroutines glViewport() and glm::ortho().

Applications

• Draw in a limited area of screen. (Note: it is often possible to define multiple viewport windows.)
• Draw using user-defined coordinates and make the drawing independent of the screen size.
• Draw multiple views of the world-space.
• Draw charts (a.k.a curve plotting).
• Draw multiple charts on the same screen.

Parts:

• World-Window (a.k.a clipping window) => Defines what the user wants to see from the world-space. Define by: \(x_{wmin}\), \(x_{wmax}\), \(y_{wmin}\) and \(y_{wmax}\).
  • \(x_{wmin}\) - Minimum X axis coordinate from window-space that can be viewed.
  • \(x_{wmax}\) - Maximum X axis coordinate from window-space that can be viewed.
  • \(y_{wmin}\) - Minimum y axis coordinate from window-space that can be viewed.
  • \(y_{wmax}\) - Maximum y axis coordinate from window-space that can be viewed.
• Viewport => Defines where the user wants to see the world-window within the graphics display window (a.k.a canvas).
  • \(x_{vmin}\) => \(0 \leq x_{vmin} \leq w\)
  • \(x_{vmax}\) => \(0 \leq x_{vmin} \leq w\)
  • \(y_{vmin}\) => \(0 \leq y_{vmin} \leq h\)
  • \(y_{vmax}\) => \(0 \leq y_{vmax} \leq h\)
• Default values of viewport:
  • \(x_{vmin} = 0\)
  • \(x_{vmax} = w\)
  • \(y_{vmin} = 0\)
  • \(y_{vmax} = h\)

Where:

• h - graphics display screen height, often in pixels.
• w - graphics display screen width, often in pixelss

The window-to-viewport transform matrix can be computed as:

Where:

Coordinates from world-space can be mapped to the screen-device space by applying the affine transform \(T_{ W \rightarrow V}\) to the world-space coordinates, designated by \(x_w\), \(y_w\).

Finally, the coordinates \(x_v\) and \(y_v\) can be determined in a more human-friendly way with the following expression:

Image distortion

The world-space propotions will only be preserved when the following predicate holds. For instance, a square in the world-space will look like a rectangle if the rations \(s_x\) and \(s_y\) are not equal.

In orther words,

Upper-left coordinate system

If the graphics API has a default upper-left coordinate system, where the origin is at the upper-left corner of the display window and Y axis is pointing downwards. The viewport matrix transform becomes:

Then, the \(x_v\), \(y_v\) coordinates can be computed as:

References and further reading:

• Khronos Group - Viewport transform
  • "The viewport transform defines the transformation of vertex positions from NDC space to window space. These are the coordinates that are rasterized to the output image. The viewport is defined by a number of viewport parameters. These parameters are set by these functions: glViewport; glDepthRange; glDepthRangef."
• MCA - 301 - Computer Graphics
• Computer Graphics Viewing Objective
• Windows and viewports
• Drawing and coordinate system

Notations for homogenous coordinate

Point-Vector - for denoting postion, location, point in space or vertex:

• It represents a point or a vertex in space, a coordinate relative to the origin of the coordinate system. Some valid affinte transform operations are translation, rotation, shear and scaling.

Vector - for denoting force, acceleration, or difference between two point-vectors, and so on:

• It represents entities, such as difference between two point-vectors, force, speed, angular speed, acceleration and so on, with a magnitude and direction. It does not make sense to apply translation and scaling to those entities, as a result the pseudo coordinate w is zero for this case.

3D Affine Transforms

3D affine transforms can be represented by matrices with the following format.

An affine transformation maps the coordinate system (X, Y, Z), which could be an object local coordinate system, to the coordinate system (X', Y', Z'), which could be a world coordinate system. In a similar manner to 2D homogeneous coordinates, an extra pseudo-coordinate w = 1 is added for expressing translations transformations as matrix multiplications, in the same way as rotations.

Combination between affine transforms

Consider two affine transforms \(A = a_{ij}\) and \(B = b_{ij}\). The product between those two transform is also an affine transform.

Rotation affine transforms

Rotation affine transforms represents rotation around some axis or direction have the following form:

Properties of Rotation Matrix Affine Transforms:

• Rotation matrices are orthogonal matrices and have the following properties:

• Where:
  • R is a rotation matrix
  • \(R^T\) is the transpose of the rotation matrix R.
  • \(R^{-1}\) is the inverse of the rotation matrix R.

3D Canonical Affine Transforms

Translation:

• Where: \(T^{-1}\) is the inverse transform.

Translation transform applied to an homogeneous vector (x, y, z, w = 1):

Scaling:

• Where: \(s_x\), \(s_y\), \(s_z\) are the scale factors for axis x, y, z and \(S^{-1}\) is the inverse transform (inverse matrix).

Scaling transform applied to a homogeneous coordinate vector:

Rotation around x axis

• \(R_x^{-1} = R_x^T\) => The inverse is equal to the transpose.

Transform \(R_x\) applied to a homogeneous vector (X, Y, Z, w = 1).

Rotation around y axis

• \(R_y^{-1} = R_y^T\)

Transform \(R_y\) applied to a homogeneous vector:

Rotation around z axis

• \(R_z^{-1} = R_z^T\)

Transform \(R_z\) applied to a homogeneous vector:

Rotation around arbirtary axis

• Formula 1 : The axis is defined by unit vector \(\hat{n}\)
  • References: (Chrobotics), (Ben-Ari - Weizmann) - derived from quaternion equations.

Where:

• \(\hat{n}\) is a unit vector that designates the direction, thus: \(\hat{n} = 1\)
• \(\hat{n} = [x_n \quad y_n \quad z_n]^T\)
  • \(x_n\), \(y_n\), \(z_n\) are the components of vector \(\hat{n}\).
• \(a = \cos( \theta / 2)\)
• \(b = x_n \cdot \sin (\theta / 2)\)
• \(c = y_n \cdot \sin (\theta / 2)\)
• \(d = z_n \cdot \sin (\theta / 2)\)

Rotation around arbirtary axis

• Reference: (Steve Rotenberg - Computer Animation)

Where:

• \(\theta\) is the rotation angle around the unit vector \(\hat{n}\)
• \(C = \sin \theta\)
• \(S = \cos \theta\)
• \(\hat{n}\) is a unit vector that designates the direction, thus: \(\hat{n} = 1\)
• \(\hat{n} = [x_n \quad y_n \quad z_n]^T\)
  • \(x_n\), \(y_n\), \(z_n\) are the components of vector \(\hat{n}\).

Testing formula in Sympy - Python CAS (Computer Algebra System):

import sympy

# t represents the angle theta
x, y, z, t = symbols('x y z t')
C, S = symbols("C S")

row1 = [ x**2 + C * (1 - x**2), x * y * (1 - C) - z * S, x * z * (1 - C)  + y * S ]
row2 = [ x * y * (1 - C) + z * S, y**2 + C * (1 - y**2), y * z * (1 - C) - x * S ]
row3 = [ x * z * (1 - C) - y * S,  y * z * (1 - C) +  x * S, z**2 + C * (1 - z**2) ]

m = Matrix([row1, row2, row3])

# C = cos(t)
# S = sin(t)

In [7]: m
Out[7]:
⎡   ⎛     2⎞    2                                          ⎤
⎢ C⋅⎝1 - x ⎠ + x     -S⋅z + x⋅y⋅(1 - C)  S⋅y + x⋅z⋅(1 - C) ⎥
⎢                                                          ⎥
⎢                       ⎛     2⎞    2                      ⎥
⎢S⋅z + x⋅y⋅(1 - C)    C⋅⎝1 - y ⎠ + y     -S⋅x + y⋅z⋅(1 - C)⎥
⎢                                                          ⎥
⎢                                           ⎛     2⎞    2  ⎥
⎣-S⋅y + x⋅z⋅(1 - C)  S⋅x + y⋅z⋅(1 - C)    C⋅⎝1 - z ⎠ + z   ⎦

# --- Determine rotation matrix around Z axis (Particular case) ------#
#

In [8]: m.subs({x: 0, y: 0, z: 1})
Out[8]:
⎡C  -S  0⎤
⎢        ⎥
⎢S  C   0⎥
⎢        ⎥
⎣0  0   1⎦

In [11]: rotZ = m.subs({x: 0, y: 0, z: 1, C: cos(t), S: sin(t) })

In [12]: rotZ
Out[12]:
⎡cos(t)  -sin(t)  0⎤
⎢                  ⎥
⎢sin(t)  cos(t)   0⎥
⎢                  ⎥
⎣  0        0     1⎦


#---- Determine rotation matrix around X axis (Particular case) ------#
#
In [14]: rotX = m.subs({x: 1, y: 0, z: 0, C: cos(t), S: sin(t) })

In [15]: rotX
Out[15]:
⎡1    0        0   ⎤
⎢                  ⎥
⎢0  cos(t)  -sin(t)⎥
⎢                  ⎥
⎣0  sin(t)  cos(t) ⎦

#---- Determine rotation matrix around Y axis (Particular case) ------#
#
In [16]: rotY = m.subs({x: 0, y: 1, z: 0, C: cos(t), S: sin(t) })

In [17]: rotY
Out[17]:
⎡cos(t)   0  sin(t)⎤
⎢                  ⎥
⎢   0     1    0   ⎥
⎢                  ⎥
⎣-sin(t)  0  cos(t)⎦

The Rodrigues' Rotation formula, named after the mathematician Olinde Rodrigues, allows rotating a vector in a 3D space, given an axis unit vector and a rotation angle around this axis. Variants of this formula can be used for determining the rotation matrix around any axis and for computing the axis-angle equivalent of a rotation matrix. This formula has wide variety of applications, including computer graphics, games, robotics, mechanical engineering and aerospace design.

Consider the following picture that contains a vector \('\boldsymbol{v}\) which is the position vector \(\boldsymbol{v}\) rotated by an angle \(\theta\) around an axis designated by the unite vector \(\boldsymbol{n}\). The vector \(\boldsymbol{v}_p\) is the projection of vector \(\boldsymbol{v}\) onto the vector \(\boldsymbol{n}\) and the vector \(\boldsymbol{v}_r\) is the rejection of vector \(\boldsymbol{v}\).

According to the previous picture it is possible to find the following expressions. Where r is the radius of the circle, which is the same as the lenght of CA.

\(\boldsymbol{v}_p\) and \(\boldsymbol{v}_r\) that

The projection vector \(\boldsymbol{v_p}\) (OC) can be determined using the following identity:

The cross product between vectors \(\boldsymbol{n}\) and vector \(\boldsymbol{v}\) is computed as:

Knowing that \(\| \boldsymbol{v}_r = \| \boldsymbol{v} \| \sin \alpha\), the following can be determined:

The vector \(\boldsymbol{n} \times (\boldsymbol{n} \times \boldsymbol{v})\) is orthogonal to both \(\boldsymbol{n}\) and \(\boldsymbol{n} \times \boldsymbol{v}\) and has opposite direction to the vector \(\boldsymbol{v_r}\). The norm/magnitude of this cross product vector is determined by using the cross product identity.

Since the vector \(\boldsymbol{n} \times (\boldsymbol{n} \times \boldsymbol{v})\) has the same magnitude as the vector \(\boldsymbol{v}_r\) and opposite direction, it is possible to find that.

The projection vector \(\boldsymbol{v}_p\) can also be computed in terms of rejection vector \(\boldsymbol{v}_r\) by using the expression \(\boldsymbol{v} = \boldsymbol{v}_p + \boldsymbol{v}_r\).

The unit vector \(\boldsymbol{e}_a\) that has the same direction as the vector \(\boldsymbol{v}_r\) and it can computed as:

The unit vector \(\boldsymbol{e}_b\) (direction CD), the same direction of vector \(\boldsymbol{n} \times \boldsymbol{v}\) is orthogonal to \(\boldsymbol{n}\), \(\boldsymbol{v}\) and \(\boldsymbol{v}_r\) is determined by normalizing the vector \(\boldsymbol{n} \times \boldsymbol{v}\).

By projecting the rotated rejection vector \('\boldsymbol{v}_r\) on the vectors \(\boldsymbol{e}_a\) and \(\boldsymbol{e}_b\), the next expression is found as follow:

The rotated vector \(\boldsymbol{v}'\) can be found by adding the rotated rejection vector \(\boldsymbol{v}_r'\) and the projection vector \(\boldsymbol{v}_p\), replacing the vector \(\boldsymbol{v}_p\) with \(\boldsymbol{v} + \boldsymbol{n} \times (\boldsymbol{n} \times \boldsymbol{v})\) and also replacing the \(\boldsymbol{e}_a\) and \(\boldsymbol{e}_b\) with its previous values.

Finally, by further simplyfing the previous expression, it possible to find the Rodriguez's formula for rotation:

The cross product can be written as a matrix multiplication by using the cross product matrix \([\boldsymbol{u}]_{\times}\) operator defined as:

Cross-product represented as matrix multiplication.

The Rodriguez's formula for rotation can be formulated as a matrix multiplication by replacing cross-product operations with matrix multiplication where \(\boldsymbol{I}\) is a 3 x 3 identity matrix.

So, the rotation matrix R such that \(\boldsymbol{v}' = R \boldsymbol{v}\) can be found as the following expression, where \(N = [\boldsymbol{n}]_{\times}\) is the cross product matrix of the rotation axis vector and \(\boldsymbol{I}\) is a 3 x 3 identity matrix. In computer graphics, it is more efficient to compute the rotation matrix only once and use to rotate all vertices of an object on the GPU side instead of applying the Rodriguez's formula to every object vertex.

Deriving the Rotation Matrix in Wolfram Mathemtica

The general rotation matrix R can be easily found by using a CAS - Computer Algebra System such as Wolfram Mathematica, SymPy (Python) library or GNU Maxima. In this case, it is used Wolfram Mathematica for determining this matrix due to its ease of use and its rule-rewriting systems that makes it easier to manipulate math formulas and generate LaTex, Mathml, ascii representation, C, C++ or Fortran code. An alternative CAS is mathics, a free and open source alternative implementation of Wolfram programming language, built on top of Python's sympy library for symbolic math computation.

Wolfram Mathematica or Mathics session:

(* Define the rotation axis *)
n = {nx, ny, nz}

(* Turns a vector in to its cross product matrix *)
cross[u_] := {{0, -u[[3]], u[[2]] }, { u[[3]], 0, -u[[1]]}, {-u[[2]], u[[1]], 0} }


In[6]:= I3 = IdentityMatrix[3]

Out[6]= {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}

In[7]:= I3 // MatrixForm

Out[7]//MatrixForm= 1   0   0

                    0   1   0

                    0   0   1

 In[8]:= xN = cross[n] (* Cross product matrix of rotation axis n*)

 Out[8]= {{0, -nz, ny}, {nz, 0, -nx}, {-ny, nx, 0}}

 In[9]:= xN // MatrixForm

 Out[9]//MatrixForm= 0     -nz   ny

                     nz    0     -nx

                     -ny   nx    0


 In[10]:= R = I3 + S * xN + (1 - C) * (xN . xN)

                            2     2
 Out[10]= {{1 + (1 - C) (-ny  - nz ), (1 - C) nx ny - nz S, 

                                                                     2     2
 >     (1 - C) nx nz + ny S}, {(1 - C) nx ny + nz S, 1 + (1 - C) (-nx  - nz ), 

 >     (1 - C) ny nz - nx S}, {(1 - C) nx nz - ny S, (1 - C) ny nz + nx S, 

                       2     2
 >     1 + (1 - C) (-nx  - ny )}}

 In[11]:= 
 In[11]:= R = I3 + S * xN + (1 - C) * MatrixPower[xN, 2] (* Note: a^2 - does not work with matrices *)

                            2     2
 Out[11]= {{1 + (1 - C) (-ny  - nz ), (1 - C) nx ny - nz S, 

                                                                     2     2
 >     (1 - C) nx nz + ny S}, {(1 - C) nx ny + nz S, 1 + (1 - C) (-nx  - nz ), 

 >     (1 - C) ny nz - nx S}, {(1 - C) nx nz - ny S, (1 - C) ny nz + nx S, 

                       2     2
 >     1 + (1 - C) (-nx  - ny )}}



 rotX = { nx -> 1, ny -> 0, nz -> 0,  C -> Cos[ \[Theta] ], S -> Sin[ \[Theta] ] };
 rotY = { nx -> 0, ny -> 1, nz -> 0,  C -> Cos[ \[Theta] ], S -> Sin[ \[Theta] ] };
 rotZ = { nx -> 0, ny -> 0, nz -> 1,  C -> Cos[ \[Theta] ], S -> Sin[ \[Theta] ] };


 (* --------------- Rotation Matrix around X axis  -------------------------------*)

 In[26]:= Rx = R /. rotX

 Out[26]= {{1, 0, 0}, {0, Cos[θ], -Sin[θ]}, {0, Sin[θ], Cos[θ]}}

 In[27]:= Rx // MatrixForm

 Out[27]//MatrixForm= 1          0          0

                      0          Cos[θ]    -Sin[θ]

                      0          Sin[θ]    Cos[θ]

 In[28]:= Rx // MatrixForm // TeXForm

 Out[28]//TeXForm= \left(
                   \begin{array}{ccc}
                    1 & 0 & 0 \\
                    0 & \cos (\theta ) & -\sin (\theta ) \\
                    0 & \sin (\theta ) & \cos (\theta ) \\
                   \end{array}
                   \right)

 (* --------------- Rotation Matrix around Y axis  -------------------------------*)

 In[31]:= Ry = R /. rotY

 Out[31]= {{Cos[θ], 0, Sin[θ]}, {0, 1, 0}, {-Sin[θ], 0, Cos[θ]}}

 In[32]:= Ry // MatrixForm

 Out[32]//MatrixForm= Cos[θ]    0          Sin[θ]

                      0          1          0

                      -Sin[θ]   0          Cos[θ]

 In[33]:= Ry // MatrixForm // TeXForm

 Out[33]//TeXForm= \left(
                   \begin{array}{ccc}
                    \cos (\theta ) & 0 & \sin (\theta ) \\
                    0 & 1 & 0 \\
                    -\sin (\theta ) & 0 & \cos (\theta ) \\
                   \end{array}
                   \right)

 (* --------------- Rotation Matrix around Z axis  -------------------------------*)

 In[34]:= Rz = R /. rotZ

 Out[34]= {{Cos[θ], -Sin[θ], 0}, {Sin[θ], Cos[θ], 0}, {0, 0, 1}}

 In[35]:= Rz // MatrixForm 

 Out[35]//MatrixForm= Cos[θ]    -Sin[θ]   0

                      Sin[θ]    Cos[θ]    0

                      0          0          1

 In[36]:= Rz // MatrixForm // TeXForm

 Out[36]//TeXForm= \left(
                   \begin{array}{ccc}
                    \cos (\theta ) & -\sin (\theta ) & 0 \\
                    \sin (\theta ) & \cos (\theta ) & 0 \\
                    0 & 0 & 1 \\
                   \end{array}
                   \right)

 (* -------------- General Rotation Matrix --------------------------------------- *)

In[41]:= Clear[n]
In[42]:= Rp = R /. { nx -> Subscript[n, x], ny -> Subscript[n, y], nz -> Subscript[n, z] };

In[44]:= Rp // MatrixForm // TeXForm

Out[44]//TeXForm= 
  \left(
  \begin{array}{ccc}
   (1-C) \left(-n_y^2-n_z^2\right)+1 & (1-C) n_x n_y-S n_z & (1-C) n_x n_z+S
     n_y \\
   (1-C) n_x n_y+S n_z & (1-C) \left(-n_x^2-n_z^2\right)+1 & (1-C) n_y n_z-S
     n_x \\
   (1-C) n_x n_z-S n_y & (1-C) n_y n_z+S n_x & (1-C)
     \left(-n_x^2-n_y^2\right)+1 \\
  \end{array}
  \right)


 In[52]:= R2 = R /. {-ny^2 - nz^2 -> nx^2 - 1, -nx^2 - nz^2 -> ny^2 - 1, -nx^2 - ny^2 -> nz^2 - 1 }

                                2
 Out[52]= {{1 + (1 - C) (-1 + nx ), (1 - C) nx ny - nz S, 

                                                                         2
 >     (1 - C) nx nz + ny S}, {(1 - C) nx ny + nz S, 1 + (1 - C) (-1 + ny ), 

 >     (1 - C) ny nz - nx S}, {(1 - C) nx nz - ny S, (1 - C) ny nz + nx S, 

                           2
 >     1 + (1 - C) (-1 + nz )}}


 In[53]:= Rp2 = R2 /. { nx -> Subscript[n, x], ny -> Subscript[n, y], nz -> Subscript[n, z] };

 In[54]:= Rp2 

                                2
 Out[54]= {{1 + (1 - C) (-1 + n  ), (1 - C) n  n  - S n , 
                               x             x  y      z

                                                                         2
 >     S n  + (1 - C) n  n }, {(1 - C) n  n  + S n , 1 + (1 - C) (-1 + n  ), 
          y            x  z             x  y      z                     y

 >     -(S n ) + (1 - C) n  n }, 
            x             y  z

                                                                          2
 >    {-(S n ) + (1 - C) n  n , S n  + (1 - C) n  n , 1 + (1 - C) (-1 + n  )}}
            y             x  z     x            y  z                     z


 In[55]:= Rp2 // MatrixForm // TeXForm

 Out[55]//TeXForm= 
    \left(
    \begin{array}{ccc}
     (1-C) \left(n_x^2-1\right)+1 & (1-C) n_x n_y-S n_z & (1-C) n_x n_z+S n_y
       \\
     (1-C) n_x n_y+S n_z & (1-C) \left(n_y^2-1\right)+1 & (1-C) n_y n_z-S n_x
       \\
     (1-C) n_x n_z-S n_y & (1-C) n_y n_z+S n_x & (1-C) \left(n_z^2-1\right)+1
       \\
    \end{array}
    \right)

The overall rotation matrix around an arbitrary axis is given by the following expression where \(n_x\), \(n_y\) and \(n_z\) are components of the rotation axis vector \(\boldsymbol{n}\) (unitary vector); \(\theta\) is the rotation angle; \(C = \sin \theta\); and \(S = \sin \theta\).

The general rotation matrix can also be expressed in this form that is obtained by replacing \(n_x^2 - 1 = - (n_y^2 + nz^2)\), \(n_y^2 - 1 = - (n_x^2 + nz^2)\) and \(n_z^2 - 1 = - (n_x^2 + ny^2)\).

Rotation matrix around X axis.

Rotation matrix around Y axis.

Rotation matrix around Z axis.

Rodriguez Formula in Alternative Format

From previous equation, the unit vector \(\boldsymbol{e}_a\) can be found as:

Then, rotated vector \(\boldsymbol{v}'\) can be found by replacing the values of previous equations in the next equation.

Alternative Format of Roriguez Formula: Finally, the alternative format of Rodriguez Formula for rotation can be determined as follow.

It can be found that the previous expression can be turned into a matrix multiplication as in shown in the next equation where \(a^T\) is the tranpose of matrix a; I 3x3 identiy matrix; \(N = [\boldsymbol{b}]_{\times}\) is the cross product matrix of rotation axis vector \(\boldsymbol{n}\). In this case, the vectors \(\boldsymbol{v}\), \(\boldsymbol{v}'\) and \(\boldsymbol{n}\) are treated as 3 x 1 (3 by 1) column matrices.

So, the general rotation matrix R that represents a rotation around an arbitrary axis such that \(\boldsymbol{v}' = R \boldsymbol{v}\) can be found as:

The rotation matrix can be determined from this expression by using Wolfram Mathematica computer algebra system.

In[1]:= n = { {nx}, {ny}, {nz} }

Out[1]= {{nx}, {ny}, {nz}}

In[2]:= Dimensions[n]

Out[2]= {3, 1}

In[3]:= n // MatrixForm

Out[3]//MatrixForm= nx

                    ny

                    nz

In[4]:= P = n . Transpose[n]

            2                           2                           2
Out[4]= {{nx , nx ny, nx nz}, {nx ny, ny , ny nz}, {nx nz, ny nz, nz }}

In[5]:= P // MatrixForm

Out[5]//MatrixForm=   2
                    nx      nx ny   nx nz

                              2
                    nx ny   ny      ny nz

                                      2
                    nx nz   ny nz   nz


In[6]:= I3 = IdentityMatrix[3]

Out[6]= {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}

In[7]:= SN = {{0, -nz, ny}, {nz, 0, -nx}, {-ny, nx, 0}};

(* C - represents cos(t), S - represents sin(t) *)
In[8]:= R = P + (I3 - P) * C + SN * S 


In[9]:= R // MatrixForm

Out[9]//MatrixForm= 

>     2            2
    nx  + C (1 - nx )        nx ny - C nx ny - nz S   nx nz - C nx nz + ny S

                               2            2
    nx ny - C nx ny + nz S   ny  + C (1 - ny )        ny nz - C ny nz - nx S

                                                        2            2
    nx nz - C nx nz - ny S   ny nz - C ny nz + nx S   nz  + C (1 - nz )


rotX = { nx -> 1, ny -> 0, nz -> 0,  C -> Cos[ \[Theta] ], S -> Sin[ \[Theta] ] };
rotY = { nx -> 0, ny -> 1, nz -> 0,  C -> Cos[ \[Theta] ], S -> Sin[ \[Theta] ] };
rotZ = { nx -> 0, ny -> 0, nz -> 1,  C -> Cos[ \[Theta] ], S -> Sin[ \[Theta] ] };


In[14]:= Rx = R /. rotX

Out[14]= {{1, 0, 0}, {0, Cos[θ], -Sin[θ]}, {0, Sin[θ], Cos[θ]}}


In[15]:= Rx // MatrixForm

Out[15]//MatrixForm= 1          0          0

                     0          Cos[θ]    -Sin[θ]

                     0          Sin[θ]    Cos[θ]


In[16]:= Ry = R /. rotY

Out[16]= {{Cos[θ], 0, Sin[θ]}, {0, 1, 0}, {-Sin[θ], 0, Cos[θ]}}

In[17]:= MatrixForm[Ry]

Out[17]//MatrixForm= Cos[θ]    0          Sin[θ]

                     0          1          0

                     -Sin[θ]   0          Cos[θ]


In[18]:= Rz = R /. rotZ

Out[18]= {{Cos[θ], -Sin[θ], 0}, {Sin[θ], Cos[θ], 0}, {0, 0, 1}}

In[19]:= Rz // MatrixForm

Out[19]//MatrixForm= Cos[θ]    -Sin[θ]   0

                     Sin[θ]    Cos[θ]    0

                     0          0          1


In[27]:= ClearAll[n]

In[28]:= toSubscr =  { nx -> Subscript[n, x], ny -> Subscript[n, y], nz -> Subscript[n, z]};

In[29]:= Rp = R /. toSubscr;


In[31]:= Rp

             2            2
Out[31]= {{n   + C (1 - n  ), n  n  - C n  n  - S n , 
            x            x     x  y      x  y      z

                                                          2            2
>     S n  + n  n  - C n  n }, {n  n  - C n  n  + S n , n   + C (1 - n  ), 
         y    x  z      x  z     x  y      x  y      z   y            y

>     -(S n ) + n  n  - C n  n }, 
           x     y  z      y  z

                                                           2            2
>    {-(S n ) + n  n  - C n  n , S n  + n  n  - C n  n , n   + C (1 - n  )}}
           y     x  z      x  z     x    y  z      y  z   z            z


In[33]:= Rp // MatrixForm // TeXForm

Out[33]//TeXForm= 
   \left(
   \begin{array}{ccc}
    C \left(1-n_x^2\right)+n_x^2 & -C n_x n_y-S n_z+n_x n_y & -C n_x n_z+S
      n_y+n_x n_z \\
    -C n_x n_y+S n_z+n_x n_y & C \left(1-n_y^2\right)+n_y^2 & -C n_y n_z-S
      n_x+n_y n_z \\
    -C n_x n_z-S n_y+n_x n_z & -C n_y n_z+S n_x+n_y n_z & C
      \left(1-n_z^2\right)+n_z^2 \\
   \end{array}
   \right)

Finally, from Wolfram mathematica session, the rotation matrix can be determined as:

The camera view transform matrix transforms world-space coordinates to camera-space coordinates. On old OpenGL versions, this matrix was can be constructed through the OpenGL routine gluLookAt or the GLM (OpenGL Mathematics) library routine glm::lookAt.

The GLM function glm::lookAt() has the following type signature.

glm::mat4 glm::lookAt(
                // Eye vector - camera position in world-space coordinates
                glm::vec3 const& eye

                // At vector - point to where camera is looking at
                // in world-space coordinates.
              , glm::vec3 const& at

                // Up vector - camera's orientation. Often set to vertical axis
                // or Y axis up = (x = 0, y = 1, z = 0) by default.
              , glm::vec3 const& up
              );

  // --------- Usage ------------------------//

  ... ... ... ... ...

glm::mat4 camera_View_transform = glm::lookAt(  eye_camera_current_position
                                              , point_to_where_camera_is_looking_at
                                              , up_vector
                                              );

  // Get shader ID of shader's view transform uniform variable.
  const GLint u_view_transform = glGetUniformLocation(prog, "u_view_transform");

  // Set View tranform shader uniform variable
  glUniformMatrix4fv(u_view_transform, 1, GL_FALSE, glm::value_ptr(camera_View_transform) );

This function has the following algorithm:

And the unit vectors, X, Y and Z are computed in the following manner:

The vector \(\vec{Z}\) is the opposite direction, represented by a unit vector, to where the camera is looking at. If this vector were an input, and was directly set by the calling code, rotating that camera viewing direction would require just rotating this vector.

Where:

• \(\vec{forward} = \vec{at} - \vec{eye}\) => Direction to where the camera is looking at (a.k.a pointing at).
• INPUTS:
  • \(\vec{eye}\) (Eye vector) => Camera's position in world-space coordinates.
  • \(\vec{at}\) (At vector) => Point in world-space coordinate where the camera is lookin at.
  • \(\vec{up}\) (Up vector) => Camera orientation, often set by default to Y axis.
• Intermediate Outputs:
  • \(\vec{X} = [X_x \quad X_y \quad X_z]\)
  • \(\vec{Y} = [Y_x \quad Y_y \quad Y_z]\)
  • \(\vec{Z} = [Z_x \quad Z_y \quad Z_z]\)
• Notation Used:
  • \(\left\| \vec{vector} \right\|\) - means the vector norm
  • \(\vec{A} \times \vec{B}\) - means vector cross product
  • \(\vec{A} \cdot \vec{B}\) - means vector dot product.

The previous matrix turns world-space into camera-space coordinates in the following mode:

Where:

• The right-side (4 x 1) column vector (input) is the world-space coordinate.
• The left-side (4 x 1) column vector (ouput) is the camera-space coordinate.

Pointing camera to specific direction

The function lookAt() points the camera view to a specific point in the space at vector in world-space coordinates. The following function could be defined for pointing the camera at specific direction (vector or normalized vector) instead of a point-vector.

glm::mat4
lookAt_direction(glm::vec3 const& eye, glm::vec3 const& forward, glm::vec3 const& up  )
{
    glm::vec3 at = eye + forward;
    return lookAt(eye, at, up);
}

Where:

• eye - (point-verctor) Camera position in world-space coordinates
• forward - (vector) Direction to where camera is looking at.
• up (vector) Camera's orientation

Algorithm Validation

Testing GLM (OpenGL Mathematic Library) in CERN's root REPL:

// -------------- Case  1 --------------------- //
//
auto eye = glm::vec3(3.0, 10.0, 20.0);
auto at  = glm::vec3(50.0, 25.0, 10.0);
auto up  = glm::vec3(0.0, 1.0, 0.0);

root [108] auto view_matrix = glm::lookAt(eye, at, up);
root [109]
root [109] show_matrix("view = ", view_matrix)

 [MATRIX] view =  =
   0.208  -0.000   0.978 -20.186
  -0.291   0.955   0.062  -9.912
  -0.934  -0.298   0.199   1.808
   0.000   0.000   0.000   1.000
root [110]

// ------------ Case 2 --------------------------//
auto eye = glm::vec3(50.0, -10.0, 0.0);
auto up = glm::vec3(0.0, 1.0, 0.0);
auto at = glm::vec3(25.0, 5.615, -200.0);
auto view_matrix = glm::lookAt(eye, at, up);

root [114]
root [114] show_matrix("view = ", view_matrix)

 [MATRIX] view =  =
   0.992   0.000  -0.124 -49.614
   0.010   0.997   0.077   9.491
   0.124  -0.077   0.989  -6.956
   0.000   0.000   0.000   1.000

Implement glm::lookAt() algorithm in Julia programming language:

# Return vector v normalized =>> Returns column vector
normalize(v) = begin
   x, y, z = v
   mag = sqrt(x * x + y * y + z * z)
   [ (x / mag) (y / mag) (z / mag)]'
end

# Vector cross product =>> Returns column vector
function cross(va, vb)
   x, y, z = va
   T = [ 0 -z y ; z 0 -x; -y x 0]
   return T * vb
end

# Vector Dot product via vector column X transpose multiplication
dot(va, vb) = va' * vb

# Note: The inputs must be column vectors
function lookAt(eye, at, up)
   # The vector Z is the direction to where the camera is looking at.
   Z = normalize( eye - at      )
   X = normalize( cross(up, Z)  )
   Y = normalize( cross(Z, X )  )
   t = [ dot(eye, X)  dot(eye, Y)  dot(eye, Z) ]

   # Define identity matrix
   matrix = [ 1.0 0.0 0.0 0.0 ; 0.0 1.0 0.0 0.0 ; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0]

   # Assign vector X to matrix's first row
   matrix[1, 1:3] = X
   # Assign vector Y to matrix's second row
   matrix[2, 1:3] = Y
   # Assign vector Z to matrix's third row
   matrix[3, 1:3] = Z
   # Assign forth (last) column to vector t
   matrix[1:3, 4] = -t
   return matrix
end

Test implementation:

#  --------- Test Case 1 --------------------------#
# =================================================#
eye = [3.0 10.0 20.0  ]'
at  = [50.0 25.0 10.0 ]'
up  = [0.0   1.0 0.00 ]' # Y axis

julia> eye
3×1 LinearAlgebra.Adjoint{Float64,Array{Float64,2}}:
  3.0
 10.0
 20.0

# Camera's view transform matrix
julia> Tview = lookAt(eye, at, up)
4×4 Array{Float64,2}:
 -0.208108  0.0       -0.978106   20.1864
 -0.291457  0.954572   0.062012   -9.91159
  0.933672  0.297981  -0.198654  -52.1466
  0.0       0.0        0.0         1.0

# Camera's view transform matrix rounded with 3 decimal digits
julia> map(x -> round(x, digits=3), Tview)
4×4 Array{Float64,2}:
  0.208   0.0    0.978  -20.186
 -0.291   0.955  0.062   -9.912
 -0.934  -0.298  0.199    1.808
  0.0     0.0    0.0      1.0

#  --------- Test Case 2 --------------------------#
# =================================================#
eye = [50.0 -10.0 0.0 ]'
up  = [0.0 1.0 0.0 ]'
at  = [25.0 5.615 -200.0]'

Tview = lookAt(eye, at, up)

julia> map(x -> round(x, digits=3), Tview)
4×4 Array{Float64,2}:
 0.992   0.0    -0.124  -49.614
 0.01    0.997   0.077    9.491
 0.124  -0.077   0.989   -6.956
 0.0     0.0     0.0      1.0

References

• CSE 167: Introduction to Computer Graphics - Lecture 3 - Coordinate Systems.
• Viewing and Modelling - CS354

The camera's projection transform matrix turns camera-space coordinates into clip-space coordinates, which are normalized NDC (Normalized-Device Coordinates). Any coordinate in the clip-space coordinate out of the range (-1.0 to 1.0) in the three axis X, Y and Z will not be shown at the screen, they will be clipped. It is worth noting that some projection matrices are not affine transforms.

Ortographic Projection Transform - (CSE167 - page 25) ; (Microsoft- WGL)

The orthographic perspective transformation preserves parallel lines and provides no sense of depth. This type of projection is suitable for charts, engineering drawing, engineering blueprints, CADs (Computer Aided Design) views and so on.

• Computed with subroutines:
  • glm::ortho (xmin, xmax, ymin, ymax, zNear, zFar) => GLM math library
  • glm::ortho (left, right, bottom, top, near, far)
  • glOrtho(left, right, bottom, top, near, far) => Legacy OpenGL

Where:

• \(x_{min}\) (left)
• \(x_{max}\) (right)
• \(y_{min}\) (bottom)
• \(y_{max}\) (top)
• \(z_{near}\) (near) => Distance from the camera coordinate system to the nearest clipping plane.
• \(z_{far}\) (far) => Distance from the camera coordinate system to the farthest clipping plane.

Observations:

• The axis Z points the opposite direction to where the camera is looking at.
• Any coordinate outside the range (\(x_{min}\), \(x_{max}\)), (\(y_{min}\), \(y_{max}\)) and (\(z_{far}\), \(z_{near}\)) will not be visible.
• Those parameters are in the camera's coordinate system.

If a camera looking at the negative direction of Z axi is positioned at (x = 0, y = 0, z = 0) in world-coordinate origin. The ortographic projection for OpenGL's canonical view volume is the identity matrix.

• \(x_{min} = -1\), \(x_{min} = +1\)
• \(y_{min} = -1\), \(z_{min} = +1\)
• \(z_{near} = +1\)
• \(z_{far} = -1\)

Ortographic Projection for 2D drawing

The ortographic projection can be used 2D drawing by setting the coordinate Z of any vertex to zero, \(z_{near} = -1\) and \(z_{far} = +1\). Then, the ortographic projection becomes:

If the transform is applied to any vertex with z = 0, the following outcome is achieved.

Perspective projection matrix (with FOV - Field-of-View)

• This matrix is generated by the following subroutines:
  • glm::perspective (FOV, aspect_ratio, zNear, zFar) => GLM library
  • glPerspective (FOV, aspect_ratio, zNear, zFar) => Old OpenGL.

The matrix elements are:

Where:

• k - is the window aspect ration, the ratio between the window width and its height. So:
  • \(k = \text{width} / \text{height}\)
• \(\theta\) - is the FOV (Field-Of-View angle)
• \(zN\) (zNear) - Distance to near plane.
• \(zF\) (zFar) - Distance to far plane.
• Only camera-space coordinates within the range \(zN\) to \(zF\) will be displayed on the screen. Anything out of this range will be clipped.

This matrix transform coordinates in the following manner:

Where:

• The right-side (4 x 1) column vector (input) is the camera-space coordinate.
• The left-side (4 x 1) column vector (ouput) is the clip-space coordinate. Any clip-space (NDC coordinates) coordinate out of the range -1.0 to 1.0 on each axis will not be displayed on the screen.

Perspective Projection Matrix (defined from View frustum) - (CSE167),

The perspective projection matrix provides a sense of depth, objects that are far from the camera will look like smaller and objects near the camera will look like smaller. Note: parallel lines are not preserved by this transform and also this transform is not affine.

• This matrix can be computed by using subroutines:
  • glm::frustum (Xmin, Xmax, Ymin, Ymax, Zmin, Zmax)
  • glm::frustum (left, right, bottom, top, near, far)

Aspect ratio and field of view (FOV) angle:

Where:

• \(x_{min}\) => (left) - Minimum X coordinate
  • Coordinate for the left clipping plane.
• \(x_{max}\) => (right) - Maximum X coordinate
  • Coordinate for the right clipping plane.
• \(y_{min}\) => (bottom) - Minimum Y coordinate
  • Coordinate for the bottom clipping plane.
• \(y_{max}\) => (top) - Maximum Y coordinate
  • Coordinate for the top clipping plane.
• \(z_{near}\) => (near)
  • Minimum Z coordinate - distance to the near depth clipping plane. The \(z_{near}\) value should never be zero, in this case, the \(z_{near}\) value should as close as possible to zero.
• \(z_{far}\) => (far)
  • Maximum Z coordinate - distance to the far depth clippling plane.

Notes:

• The camera is looking at the opposite direction of Z axis.
• Any point in camera-space coordinate outside of the previous range will not be visible on the screen.
• The parameters \(z_{far}\) and \(z_{near}\) must be positive.

Perspective and orthogonal transform matrix formula validation - in Julia Language

Running GLM math library in CERN's Root RPL:

root [0] .I .

#include <glm/glm.hpp>
#include <glm/gtc/matrix_transform.hpp>
#include <glm/gtc/type_ptr.hpp>
#include <glm/gtx/string_cast.hpp>
#include <glm/gtx/io.hpp>

std::cout << std::setprecision(5);

// Note: GLM matrices are stored in Column-major order
void show_matrix(const char* label, glm::mat4 const& m){
    std::cout << "\n [MATRIX] " << label << " = " << '\n';
    std::cout << std::fixed << std::setprecision(5);
    for(size_t i = 0; i < 4; i++)
    {
        for(size_t j = 0; j < 4; j++)
        {
            std::cout << std::setw(10) << m[j][i];
        }
        std::cout << '\n';
    }
}


float xmin, xmax, ymin, ymax, zmin, zmax;

// --------- Test Case 1 --------------------------//
//
root [58] xmin = 4.0, xmax = 10.0, ymin = 0.1, ymax = 25.0, zmin = 8.0, zmax = 20.0;

root [59] glm::mat4 T = glm::frustum(xmin, xmax, ymin, ymax, zmin, zmax);

root [60] std::cout << " T = " << T << '\n'
 T =
[[    2.667,    0.000,    2.333,    0.000]
 [    0.000,    0.643,    1.008,    0.000]
 [    0.000,    0.000,   -2.333,  -26.667]
 [    0.000,    0.000,   -1.000,    0.000]]


// ------- Test Case 2 ------------------------------//

root [93] xmin = -20.0, xmax = 50.0, ymin = -10.0, ymax = 25.0, zmin = 0.1, zmax = 100.0;
root [94]
root [94] glm::mat4 T = glm::frustum(xmin, xmax, ymin, ymax, zmin, zmax)
(glm::mat4 &) @0x7f73c3d5d490
root [95]
root [95] show_matrix("T", T)

 [MATRIX] T =
   0.00286   0.00000   0.42857   0.00000
   0.00000   0.00571   0.42857   0.00000
   0.00000   0.00000  -1.00200  -0.20020
   0.00000   0.00000  -1.00000   0.00000
root [96]


// ------- Test Case 3 - Ortographics Matrix  -----------//

root [96] xmin = -20.0, xmax = 50.0, ymin = -10.0, ymax = 25.0, zmin = 0.1, zmax = 100.0;
root [97]
root [97] glm::mat4 T = glm::ortho(xmin, xmax, ymin, ymax, zmin, zmax)
(glm::mat4 &) @0x7f73c3d5d4d0
root [98]
root [98] show_matrix("T", T)

 [MATRIX] T =
   0.02857   0.00000   0.00000  -0.42857
   0.00000   0.05714   0.00000  -0.42857
   0.00000   0.00000  -0.02002  -1.00200
   0.00000   0.00000   0.00000   1.00000

Test in Julia programming language:

 function frustum(xmin, xmax, ymin, ymax, zmin, zmax)
      a11 = 2 * zmin / (xmax - xmin)
      a22 = 2 * zmin / (ymax - ymin)
      a13 = (xmax + xmin) / (xmax - xmin)
      a23 = (ymax + ymin) / (ymax - ymin)
      a33 = -(zmax + zmin) / (zmax - zmin)
      a34 = -2 * zmax * zmin / (zmax - zmin)
      T = [ a11 0 a13 0 ; 0 a22 a23 0; 0 0 a33 a34; 0 0 -1 0]
      return T
  end

  function ortho(xmin, xmax, ymin, ymax, zmin, zmax)
        a11 = 2 / (xmax - xmin)
        a22 = 2 / (ymax - ymin)
        a33 = -2 / (zmax - zmin)
        a14 = -(xmax + xmin ) / (xmax - xmin)
        a24 = -(ymin + ymax) / (ymax - ymin)
        a34 = -(zmax + zmin) / (zmax - zmin)
        T = [ a11 0 0 a14; 0 a22 0 a24 ; 0 0 a33 a34; 0 0 0 1]
        return T
  end


 ## --------- Test Case 1 --------------------------##

 julia> xmin = 4.0; xmax = 10.0; ymin = 0.1; ymax = 25.0; zmin = 8.0; zmax = 20.0;

 julia> frustum(xmin, xmax, ymin, ymax, zmin, zmax)
 4×4 Array{Float64,2}:
  2.66667  0.0       2.33333    0.0
  0.0      0.64257   1.00803    0.0
  0.0      0.0      -2.33333  -26.6667
  0.0      0.0      -1.0        0.0

 ## --------- Test Case 2 --------------------------##

julia> xmin = -20.0; xmax = 50.0; ymin = -10.0; ymax = 25.0; zmin = 0.1; zmax = 100.0;

julia> frustum(xmin, xmax, ymin, ymax, zmin, zmax)
4×4 Array{Float64,2}:
 0.00285714  0.0          0.428571   0.0
 0.0         0.00571429   0.428571   0.0
 0.0         0.0         -1.002     -0.2002
 0.0         0.0         -1.0        0.0

## -------- Test Case 3 - Ortographic Projection ---##

julia> xmin = -20.0; xmax = 50.0; ymin = -10.0; ymax = 25.0; zmin = 0.1; zmax = 100.0;

julia> ortho(xmin, xmax, ymin, ymax, zmin, zmax)
4×4 Array{Float64,2}:
 0.0285714  0.0         0.0      -0.428571
 0.0        0.0571429   0.0      -0.428571
 0.0        0.0        -0.02002  -1.002
 0.0        0.0         0.0       1.0

References

• GlFrustum - Khronos Group (Documentation of function glFrustum())
  • https://www.khronos.org/registry/OpenGL-Refpages/gl2.1/xhtml/glFrustum.xml
• Projection and viewing - Computer Graphics - CSE167 - Lecture 4 [BEST]
  • https://cseweb.ucsd.edu/classes/wi18/cse167-a/lec4.pdf
  • glFrustum() - perspective transform formula located at page 20.
• View Frustum Clipping - University of Texas - CS354 [BEST]
  • https://www.cs.utexas.edu/~fussell/courses/cs354-fall2015/lectures/lecture9.pdf
  • glFrustum() - formula located at page
• glOrtho - Khronos Grouo (Documentation of function glOrtho())
  • https://www.khronos.org/registry/OpenGL-Refpages/gl2.1/xhtml/glOrtho.xml
• Ortographic Projection
  • https://en.wikipedia.org/wiki/Orthographic_projection
• Computergrafik - Matthias Zwicker Universität Bern - Herbst 2016
  • https://www.cs.umd.edu/~zwicker/courses/computergraphics/03_Projection.pdf
• Matrix4x4.Ortho - Unity 3D engine
  • https://docs.unity3d.com/ScriptReference/Matrix4x4.Ortho.html
• Microsoft WGL Docs - glFrustum() [BEST]
  • https://docs.microsoft.com/en-us/windows/win32/opengl/glfrustum
• Microsoft WGL Docs - glOrtho() [BEST]
  • https://docs.microsoft.com/en-us/windows/win32/opengl/glortho
• CSE - 420 Computer Graphics - California State University [BEST]
  • http://cse.csusb.edu/tongyu/courses/cs420/notes/viewing2.php

Concatenating / Combinating of affine transforms

Multiple affine transforms can be combined or concatenated as a single affine transform, just by multiplying matrices. Considering the following sequence of 3D affine transforms, such as:

1. scale(3, 3, 3) - scale by 3 on all axis.
2. translate(x = 25, y = 10, z = 20) - translate to point (25, 30, 2)
3. rotX(30) - rotate 30 degrees around X axis
4. rotZ(54) - rotate by 54 degrees around Z axis.

Let V = [ x y z 1 ]' to be (4 x 1 - 4 rows and 1 column) a column vector containing the current vertice homogeneous coordinates.

Those transforms can be combined in the following way:

• Node: the symbol (·) means multiplication.

// Where:
//  =>  scale(3, 3, 3) is the scaling affine transform matrix.
//  =>  Va is the vertex coordinate after applying rotZ(54)
//
Va = rotZ(54) · V

// Where:
// => translate(25, 10, 20) is the translation affine transfom matrix.
// => Vb is the result, vertex coodinate after the transform was applied.
Vb = rotX(30) · Va

Vc = translate(25, 10, 20) · Vb
Vd = scale(3, 3, 3) · Vc

// V' (V - prime) - overall transformation result (end result)
V' = Vd

By performing algebraics substitution, the overall transform result becomes:

STEP 1:
 Vb = rotX(30) · Va
 Vb = rotX(30) · ( rotZ(54) · Va )

STEP 2:
  Vc = translate(25, 10, 20) · Vb
  Vc = translate(25, 10, 20) ·  ( rotX(30) · ( rotZ(54) · Va ) )

STEP 3:
  Vd = scale(3, 3, 3) · Vc
  Vd = scale(3, 3, 3) · ( translate(25, 10, 20) ·  rotX(30) ·  rotZ(54) · Va )

STEP 4:  - End result
  V' = Vd
  V' = scale(3, 3, 3) · translate(25, 10, 20) ·  rotX(30) ·  rotZ(54) · Va 

The final vertex position can be stated as:

   V' = scale(3, 3, 3) · translate(25, 10, 20) ·  rotX(30) ·  rotZ(54) · Va 

OR:

   V' =  [ scale(3, 3, 3) · translate(25, 10, 20) ·  rotX(30) ·  rotZ(54) ] · V

All the sequence of transformations can be combined as a single transform, designated by T, which is the multiplication between all transforms in reverse order that they were applied. The last transform is multiplied first and first transform is multiplied last. In the same that matrix multiplication is not commutative, the operation result depends on the order that of matrices, transform combinations are not commutative, they also depends on the order that transforms were applied.

V' = T · V

Where:
   T =  scale(3, 3, 3) · translate(25, 10, 20) ·  rotX(30) ·  rotZ(54)

Multiple transforms can be combined as:

Tranform_combined = Transform[1] · Transform[2] · ... · Transform[N-1] · Transform[N]

Combining 2D canonical transforms in Sympy - Python CAS

Sympy - Python CAS (Computer Algebra System) can be used for debugging and testing affine transforms in a symbolic way.

Define 2D canonical affine transforms:

from sympy import *

x,  y  = symbols("x y")
dx, dy = symbols("dx dy")
sx, sy = symbols("sx sy")
t      = symbols("t")  # t represents the rotation angle theta (θ)


# Column vector representing - generic Vertex (aka point) in homogeneous coordinate
#
v = Matrix([x, y, 1])

>>> pprint(v)
⎡x⎤
⎢ ⎥
⎢y⎥
⎢ ⎥
⎣1⎦

# Translation of (dx, dy) coordinate transform
>>> Translate = Matrix([[1, 0, dx], [0, 1, dy], [0, 0, 1]])

>>> pprint(Translate)
⎡1  0  dx⎤
⎢        ⎥
⎢0  1  dy⎥
⎢        ⎥
⎣0  0  1 ⎦


# Scale affine transform
>>> Scale = Matrix([[sx, 0, 0], [0, sy, 0], [0, 0, 1]])

>>> pprint(Scale)
⎡sx  0   0⎤
⎢         ⎥
⎢0   sy  0⎥
⎢         ⎥
⎣0   0   1⎦


# Rotation around Z axis
>>> RotZ  = Matrix([[ cos(t), -sin(t), 0], [ sin(t), cos(t), 0], [0, 0, 1]])

>>> pprint( RotZ )
⎡cos(t)  -sin(t)  0⎤
⎢                  ⎥
⎢sin(t)  cos(t)   0⎥
⎢                  ⎥
⎣  0        0     1⎦

Apply single transforms to vertex v.

# Apply translation affine transform to vertex V
#-------------------------------------------
>>> pprint( Translate * v )
⎡dx + x⎤
⎢      ⎥
⎢dy + y⎥
⎢      ⎥
⎣  1   ⎦

# Apply scale to vertex v
#-------------------------------------------

>>> Scale * v
Matrix([
[sx*x],
[sy*y],
[   1]])

>>> pprint( Scale * v )
⎡sx⋅x⎤
⎢    ⎥
⎢sy⋅y⎥
⎢    ⎥
⎣ 1  ⎦

>>> pprint( (Scale * v).subs({sx: 10, sy: 10}) )
⎡10⋅x⎤
⎢    ⎥
⎢10⋅y⎥
⎢    ⎥
⎣ 1  ⎦


>>> pprint( (Scale * v).subs({sx: 10, sy: 20}) )
⎡10⋅x⎤
⎢    ⎥
⎢20⋅y⎥
⎢    ⎥
⎣ 1  ⎦


# Apply rotation to vertex v
#------------------------------

>>> v_rot = RotZ * v

>>> pprint( v_rot )
⎡x⋅cos(t) - y⋅sin(t)⎤
⎢                   ⎥
⎢x⋅sin(t) + y⋅cos(t)⎥
⎢                   ⎥
⎣         1         ⎦

>>> pprint( v_rot.subs(t, 0) )
⎡x⎤
⎢ ⎥
⎢y⎥
⎢ ⎥
⎣1⎦

# Rotate by 90 degrees (pi / 2 radians)
>>> pprint( v_rot.subs(t, pi / 2) )
⎡-y⎤
⎢  ⎥
⎢x ⎥
⎢  ⎥
⎣1 ⎦

# Rotate by 60 degreees (pi / 3)
>>> pprint( v_rot.subs(t, pi / 3) )
⎡x   √3⋅y⎤
⎢─ - ────⎥
⎢2    2  ⎥
⎢        ⎥
⎢√3⋅x   y⎥
⎢──── + ─⎥
⎢ 2     2⎥
⎢        ⎥
⎣   1    ⎦

Combine 2 tranforms in multiple orders and check the outcome.

# Combine transforms: 1st - translation; 2nd - scale.
#--------------------------------------------------

  >>> T = Scale * Translate

  >>> pprint(T)
  ⎡sx  0   dx⋅sx⎤
  ⎢             ⎥
  ⎢0   sy  dy⋅sy⎥
  ⎢             ⎥
  ⎣0   0     1  ⎦


  # Apply to vertex
  >>> pprint( T * v )
  ⎡dx⋅sx + sx⋅x⎤
  ⎢            ⎥
  ⎢dy⋅sy + sy⋅y⎥
  ⎢            ⎥
  ⎣     1      ⎦


# Combine transforms: 1st - scale; 2nd - translation.
#--------------------------------------------------

   >>> T = Translate * Scale

   >>> pprint( T )
   ⎡sx  0   dx⎤
   ⎢          ⎥
   ⎢0   sy  dy⎥
   ⎢          ⎥
   ⎣0   0   1 ⎦

   # Apply to vertex
   >>> pprint(T * v)
   ⎡dx + sx⋅x⎤
   ⎢         ⎥
   ⎢dy + sy⋅y⎥
   ⎢         ⎥
   ⎣    1    ⎦


# Combine transforms: 1st - rotation; 2nd - translation.
#--------------------------------------------------

   >>> T = Translate * RotZ

   >>> pprint( T )
   ⎡cos(t)  -sin(t)  dx⎤
   ⎢                   ⎥
   ⎢sin(t)  cos(t)   dy⎥
   ⎢                   ⎥
   ⎣  0        0     1 ⎦


   # Apply to vertex
   >>> pprint( T * v )
   ⎡dx + x⋅cos(t) - y⋅sin(t)⎤
   ⎢                        ⎥
   ⎢dy + x⋅sin(t) + y⋅cos(t)⎥
   ⎢                        ⎥
   ⎣           1            ⎦

# Combine transforms: 1st - translation ; 2nd - rotation
#------------------------------------------------------

   >>> T = RotZ * Translate

   >>> pprint( T )
   ⎡cos(t)  -sin(t)  dx⋅cos(t) - dy⋅sin(t)⎤
   ⎢                                      ⎥
   ⎢sin(t)  cos(t)   dx⋅sin(t) + dy⋅cos(t)⎥
   ⎢                                      ⎥
   ⎣  0        0               1          ⎦


   # Apply to vertex
   >>> pprint( T * v )
   ⎡dx⋅cos(t) - dy⋅sin(t) + x⋅cos(t) - y⋅sin(t)⎤
   ⎢                                           ⎥
   ⎢dx⋅sin(t) + dy⋅cos(t) + x⋅sin(t) + y⋅cos(t)⎥
   ⎢                                           ⎥
   ⎣                     1                     ⎦

Combine transforms in the following order: 1st - rotation; 2nd - scaling; 3rd - translation:

>>> T = Translate * Scale * RotZ

>>> pprint( T )
⎡sx⋅cos(t)  -sx⋅sin(t)  dx⎤
⎢                         ⎥
⎢sy⋅sin(t)  sy⋅cos(t)   dy⎥
⎢                         ⎥
⎣    0          0       1 ⎦

>>> pprint( T * v )
⎡dx + sx⋅x⋅cos(t) - sx⋅y⋅sin(t)⎤
⎢                              ⎥
⎢dy + sy⋅x⋅sin(t) + sy⋅y⋅cos(t)⎥
⎢                              ⎥
⎣              1               ⎦

In a 3D scene, the following typical transformations happens to the vertex coordinates of some object to be rendered. The model matrix transform, \(T_m\), converts the object vertex coordinates from the local-space (a.k.a object-space) to world-space, then the view-matrix, \(T_v\), converts vertex world-space coordinates to the camera-space, finally the projection matrix transform, \(T_p\), turns the camera-space coordinates into clip-space coordinates (NDC - Normalized Device Coordinates) within the range -1.0 to 1.0 on each axis. Any NDC coordinate out of the range -1.0 to 1.0 will not be visible on the canvas, graphics window screen.

             +---------+              +---------+ Space      +---------------+  Clip-Space
 Local       |         |  World       |         |  Camera    |               |  (NDC - Coordinates)
Space        |   Model |  Space       |  View   |  Space     |  Projection   | 
   +------->>+  Matrix +------>>----->+ Matrix  +---------->>+   Matrix      +------->>>
 V           |         |   Vw         |         |   Vc       |               |  Vndc
             | Tm      |              |  Tv     |            |    Tp         |
             +---------+              +---------+            +---------------+


   Tm - (4x4) Model matrix transform 
   Tv - (4x4) View matrix transform 
   Tp - (4x4) Projection matrix transform. 

   V    - Vertex coordinate in local space 
   Vw   - Vertex coordinate in world-space
   Vc   - Vertex coordinate in camera-space coordinate (a.k.a eye-space)
   Vndc - Vertex coordinate in clip-space coordinate (NDC)

• Transform from local-space coordinates to world-space coordinates
  • Where:
    • \(T_m\) - is the model-matrix, often comprised of rotation, scaling and translation concatened canonical transforms.
    • v - Vertex coordinate in local-space.
      • \(v = \begin{bmatrix} x & y & z & 1 \end{bmatrix}^T\)
    • \(v_w\) - Vertex coordinate in world-space
      • \(v_w = \begin{bmatrix} x_w & y_w & z_w & 1 \end{bmatrix}^T\)

• Transform from world-space coordinates to camera-space coordinates (a.k.a eye space).
  • \(T_v\) - is the view matrix, computed using the eye vector, camera position; up vector - camera orientation; at point-vector, point which the camera is looking at.

• Transform from camera-space to clip-space coordinates.
  • Where:
    • \(T_p\) - is the projection matrix which can be either the perspective projection matrix or ortographic projection matrix.
    • \(v_{\text{ndc}}\) - is the vertex NDC coordinate.

• Overall scene transform from local-space to clip-space coordinates.

• MVP transform - The overall transform is often designated as Model-View-Projection transform ( \(T_{mvp}\) ).

• The Model-View transform is the product between the model and view transforms.

• Normal transform matrix ( \(T_{normal}\) ) is used for illumination calculations, such as Phong illumination model, and is computed as the transpose of model-view matrix inverse. This transform maps vertex normal vector in the local-space coordinates to camera-space coordinates, also known as eye-space coordinates.
  • Where:
    • n - is the normal vector in local-space coordinates.
      • \(n = \begin{bmatrix} x_n & y_n & z_n & 0 \end{bmatrix}^T\)
    • \(n_c\) is the normal vector in camera-space (eye-space) coordinates.

References

• Chapter 8 - OpenGL – A 3D Graphics API 8.1 - Pipeline architecture
  • https://people.cs.clemson.edu/~dhouse/courses/817/notes/OpenGL-pipeline.pdf
• CS315 Lab 3: 3D Transformations
  • https://www.cs.uregina.ca/Links/class-info/315/WebGL2/Lab3/
• Appendix F - Homogeneous Coordinates and Transformation Matrices / OpenGL Programming Guide
  • https://www.glprogramming.com/red/appendixf.html
• GLSL Programming/Vertex Transformations - Wikibooks
  • https://en.wikibooks.org/wiki/GLSL_Programming/Vertex_Transformations
• OpenGL Normal Vector Transformation - Song ho Ahn
  • http://www.songho.ca/opengl/gl_normaltransform.html
• The normal matrix - Light House 3D
  • https://www.lighthouse3d.com/tutorials/glsl-12-tutorial/the-normal-matrix/
• Normal Vector Transform
  • https://www.cs.upc.edu/~robert/teaching/idi/normalsOpenGL.pdf

A 2D scene can be generated as particular case of a 3D scene by setting the vertex Z coordinate to always zero, replacing the view transform matrix with identity matrix and setting the projection matrix as the ortographic projection matrix. This approach only works if the graphics API default coordinate system is the same as the OpenGL one, where the Y and X axis are on the screen and the Z axis positive direction goes from the screen towards the viewer (computer user).

By following the previous approach for a 2D scene, the vertex transformation becomes:

The model matrix is a combination of translation affine transforms with z set as zero, scaling affine transforms with Z set as zero and rotations around Z axis.

Some matrices that can be used for building the model matrix are:

• 2D translation

• 2D scaling

• 2D rotation (rotation around Z axis)

Let A and B be two reference frames in euclidean three dimensional space \(\mathbb{R}^3\) whose orthonormal basis vectors are \(\{\hat{a}_1,\hat{a}_2, \hat{a}_2 \}\) and \(\{\hat{b}_1, \hat{b}_2, \hat{b}_2 \}\) where \(\hat{a}_1\) corresponds to the X axis in the frame A, \(\hat{a}_2\) corresponds to the Y axis in the frame A, and \(\hat{a}_3\) corresponds to Z axis in frame A. The basis vectors satisfy the following properties, since they are orthonormal vectors, which means that they are unitary orthogonal (aka - also known as perpendicular) vectors.

A DCM - Direction Cosine Matrix, sometimes also called Direct Cosine Matrix, is a matrix that describes the attitude, also known as orientation, of a reference frame with respect to another reference frame. In this case \(\text{DCM}_{BA}\) is the matrix that describes the orientation of a rotated reference frame B with respect to reference frame A, in other words, the orientation of frame B relative to frame A, as shown below where the matrix \(\text{DCM}_{BA} = c_{ij}\) is expresses the basis vectors of frame B in terms of basis vectors of frame A.

It is noticiable that \(c_{ij} = \hat{b}_i \cdot \hat{a}_j\), for instance

By expressing the dot product \(\hat{b}_i \cdot \hat{a}_j\) in terms of the cosine between vector angles, it is possible to find the next relation.

Where:

• \(\cos(\hat{b}_i, \hat{a}_j)\) is the cosine of the angle between the vectors \(\hat{b}_i\) and \(\hat{a}_j)\).
• \(\theta_{ij}\) is the angle between vectros \(\hat{b}_i\) and \(\hat{a}_j\) .

So, the previous expression in matrix form becomes:

So, it can be found that the matrix \(\text{DCM}_{BA}\) is given by the following expression where the first formula expresses the matrix in terms of coefficients \(c_{ij}\), the second formula expresses the matrix in terms of the dot product between reference frames basis vectors; the third states is defined in terms of the cosine between basis vectors of reference frames A and B.

The \(\text{DCM}_{BA}\) could also be expressed in more concise way in this form.

Where:

• \(^A\hat{b}_1^T\) is the transpose of frame A coordinates of vector \(\hat{b}_1\). In other words, it is projection of \(\hat{b}_1\) onto all basis vectors of reference frame A. Note: all vectors are assumed to be column vectors by default.
• \(^A\hat{b}_2^T\) is the transpose of frame A coordinates of vector \(\hat{b}_2\).
• \(^A\hat{b}_3^T\) is the transpose of frame A coordinates of vector \(\hat{b}_3\).

Coordinate Transformation between frames - A point-vector \(\boldsymbol{r}\) representing a position in the space can be expressed in both frames A and B as linear combination of basis vectors.

The coordinates of the vector \(\boldsymbol{r}\) in referece frame B can be found by taking the dot product between the vector and each basis vectors of frame B.

So, the it can be found that the previous expression can be reduced to the following expression and the matrix \(\text{DCM}_{BA}\) is the rotation matrix \({}^{B}_{A}R = \text{DCM}_{BA}\) that maps coordinates from frame A to frame B.

Where:

• \(^B\boldsymbol{r}\) - is the coordinate of position vector in frame B (rotated frame)
• \(^A\boldsymbol{r}\) - is the coordinate of position vector in frame A.

Opposite Transformation - It is often desirable to find the rotation matrix that maps coordinates from roted frame B to the frame A. It can be acomplished by taking the dot product between the previous position vector \(\boldsymbol{r}\) and each basis vector of frame A.

The previous expression can be simplified as:

It is the same as a multiplication by the transpose of the direct cosine matrix.

So, it is possible to find that the rotation matrix that maps coordinates from frame B to frame A is the tranpose of the DCM - direct cosine matrix \(\text{DCM}_{BA}\).

The rotation matrix \({}^{A}_{B}R\) can be quick determined using the following formula where each column of the rotation matrix is a basis vector of the rotated frame B expressed in terms of basis vectors of frame A. For instance, the first column of this rotation matrix is a column vector that contains the coordinates of basis vector \(\hat{b}_1\) in reference frame A, which is the same as the projection of this basis vector onto each basis vectors of frame A.

Example - Rotation around Z axis

Given two frames A, whose orthornmal basis vectors are \({a_1, a_2, a_3}\), and B whose basis vectors are \({b_1, b_2, b_3}\). The reference frame A is static (innertial reference frame) and the reference frame B is the frame A rotated counterclockwise around the Z axis by an angle \(\theta\).

APPROACH 1 => Finding the Direction Cosine Matrix \(\text{DCM}_{BA}\) by using dot product or angle between vectors

Find the the dot products between basis vectors.

• \(\hat{b}_1 \cdot \hat{a}_1 = \cos( \hat{b}_1, \hat{a}_1) = \cos \theta\)
• \(\hat{b}_1 \cdot \hat{a}_2 = \cos( \hat{b}_1, \hat{a}_2) = \cos(90 - \theta) = \sin \theta\)
• \(\hat{b}_1 \cdot \hat{a}_4 = \cos( \hat{b}_1, \hat{a}_3) = \cos(90) = 0\)
• \(\hat{b}_2 \cdot \hat{a}_1 = \cos( \hat{b}_2, \hat{a}_1) = \cos (90 + \theta) = - \sin \theta\)
• \(\hat{b}_2 \cdot \hat{a}_2 = \cos( \hat{b}_2, \hat{a}_2) = \cos \theta\)
• \(\hat{b}_2 \cdot \hat{a}_3 = \cos( \hat{b}_2, \hat{a}_3) = \cos 90 = 0\)
• \(\hat{b}_3 \cdot \hat{a}_1 = \cos( \hat{b}_3, \hat{a}_1) = \cos 90 = 0\)
• \(\hat{b}_3 \cdot \hat{a}_2 = \cos( \hat{b}_3, \hat{a}_2) = \cos 90 = 0\)
• \(\hat{b}_3 \cdot \hat{a}_3 = \cos( \hat{b}_3, \hat{a}_3) = \cos 0 = 1\)

So, the DCM matrix becomes:

The rotation matrix \({}^{A}_{B}R\) that maps vectors from frame B, rotated reference frame, to frame A, static reference frame, can be found by transposing the previous matrix

APPROACH 2 => Finding the Direction Cosine Matrix \(\text{DCM}_{BA}\) by projecting basis vectors.

STEP 1 => Write each basis vector of frame B in terms of basis vectors of frame A, in other words write the frame A coordinates of each basis vector of B by projecting each of them onto frame A basis vectors.

STEP 2 => Build the DCM matrix using the previous vectors.

The rotation matrix that maps vectors from frame B to frame A or \({}^{A}_{B}R\) can be found by tranposing the previous matrix.

The rotation matrix \({}^{A}_{B}R\) could also be found in this way.

References:

• [PDF] - Lecture Notes - Modelling of Mechanical Systems (Thomas Bak, 2002)
• DCM Tutorial – An Introduction to Orientation Kinematics – Starlino Electronics
• Rotation Matrix (Direction Cosine Matrix (DCM)) | Academic Flight
• Direction Cosine Matrix (DCM) — Introduction to Symbolic Computational Dynamics
• Attitude representation in inertial navigation · VectorNav
• Direction cosine matrix · Sea of Tranquility
• Kinematics: On Direction Cosine Matrices | IntechOpen
• Convert rotation angles to direction cosine matrix - MATLAB angle2dcm
• Attitude Representations – Understanding Direct Cosine Matrices, Euler Angles and Quaternions – Steven Dumble | PhD
• Convert rotation angles to direction cosine matrix - MATLAB angle2dcm
• Create rotation angles from direction cosine matrix - MATLAB dcm2angle

In order to properly use a graphics application programming interface, one must know the default coordinate system and conventions defined by the graphics API. This convention includes: the coordinate system origin, which may be on the upper left corner of drawing window or lower left corner of drawing window; the default direction of axis and whether the coordinate system is RHS (Right-Handed System) or LHS (Left-Handed System).

OpenGL Default Coordinate System

OpenGL API uses a default coordinate system with the origin place at the screen center, X axis as the horizontal axis, Y axis as the vertical axis and Z axis going from the screen center towards the user. This default coordinate system NDC is normalized, each dimension X, Y and Z has the range from -1.0 to +1.0. Any vertex with any dimension out of this range is not displayed on the screen.

In this coordinate system, the coordinate of points: point A at the OpenGL window top right corner is (x = 1, y = 1); the coordinate of point B is (x = 1,y= -1); the coordinate of the point C is (x = -1, y = -1).

OpenGL Coordinate System X Math-textbook Coordinate System

OpenGL default coordinate system should not be confused with the coordinate system commonly found in math and physics textbooks, where Z axis is the vertical axis, Y axis is the horizontal axis and the X axis is going from the screen twoards the screen observer (out of screen). The advantage of OpenGL coordinate system over the coordinate system used by textbooks is that 2D drawing requires just ignoring the Z axis, setting any Z coordinate to zero. If one does not like the OpenGL default coordinate system, it is possible to obtain the math textbook coordinate system just by applying affine rotation transform the OpenGL default coordinate system.

The transform matrix that maps coordinates from the math and physics style coordinate system (M) to the OpenGL coordinate (C) can be obtained by:

Where:

• \(Rot_x\) - Is the X axis rotation matrix.
• \(Rot_z\) - Is the Z axis rotation matrix.
• \(T_{M -> C}\) - Is the affine matrix transform that maps coordinates from math-and-physics coordinate system to OpenGL coordinate system.
• M - designates the coordinate system similar to the math-and-physics textbooks as in the right side of the previous picture.
• C - designates the openGL coordinate system.

Coordinates from math textbook style coordinate systems can be mapped to OpenGL by using the previous affine transform matrix.

With the previous equation it is possible to verify that the Y axis, vector (0, 1, 0) from the M coordinate system corresponds to the axis X (1, 0, 0) in the openGL (C) coordinate system.

The same procedure can be applied to X axis in (M) coordinate system, which is mapped to Z axis in OpenGL coordinate system and to Z axis from (M) coordinate system which is mapped to Y axis in OpenGL coordinate system.

Alternative Approach - using basis vectors

The rotation matrix that maps from math and physics reference frame M to Opengl reference frame, designated as C, is a matrix whose columns are the basis vectors (unit vectors) of each axis of reference frame M expressed in coordinates of reference frame C.

Where:

• \(^C\hat{X}_M\) - is the X axis basis vector of frame M expressed in coordinates of reference frame C (OpenGL reference frame).
• \(^C\hat{Y}_M\) - is the Y axis basis vector of frame M expressed in coordinates of reference frame C (OpenGL reference frame).
• \(^C\hat{Z}_M\) - is the Z axis basis vector of frame M expressed in coordinates of reference frame C (OpenGL reference frame).

So, the rotation matrix that maps from frame M to frame C becomes:

The homogenous 4 x 4 by matrix that maps from M frame M to frame C (OpengL frame) can be obtained as follow:

Coordinate System positions from other computer graphics APIS

The following diagram shows the possible default coordinate systems positions from many computer graphics APIS:

• (C) => Center of display screen
  • Default position of OpenGL NDC coordinate system (clip-space).
• (B) => Lower-left corner of display screen
  • Default position of OpenGL screen coordinate system
• (U) => Upper-left corner of display screen
  • => The Y axis is inverted, it is positive in the downward direction instead of upward direction.
  • Default position of most 2D graphics APIs including: Html5 canvas; SVG; Java AWT 2D.

Transform matrix, \(T_{B \rightarrow U}\), that maps coordinates from lower-left corner (B) coordinate system to upper-left-corner coordinate system.

Transform matrix, \(T_{C \rightarrow U}\), that maps coordinates from (C) window-center coordinate system to upper-left corner (U) coordinate system. This transform is useful for drawing with origin on screen center when using gaphics APIs which the default coordinate system is located at the upper-left corner.

Transform matrix, \(T_{B \rightarrow C}\), that maps coordinates from lower-left corner (B) to screen-center origin, OpenGL normalized coordinate system (C).

The matrix transform, \(T_{C \rightarrow B}\), that maps from coordinate system C to coordinate system B can be computed as the inverse of \(T_{B \rightarrow C}\).

OpenGL ES (embedded version of OpenGL) and WebGL APIs no longer support matrices in row-major order memory layout, where elements of rows are contiguous in memory. Instead, those APIs only support matrices in column-major memory layout, where elements of columns are contiguous in memory.

So, in modern OpenGL, the following API no longer works with boolean value 'true' that indicates that the matrix passed as a argument is in row-major oder format:

// WebGL API - (NO LONGER POSSIBLE =>> ERROR!) 
gl.uniformMatrix4fv(u_model, true, matrix_in_row_major_order); 

 // OpenGL ES - API (NO LONGER POSSIBLE!!!! =>> ERROR!!!)
 // Where u_model is a shader uniform variable 
 glUniformMatrix4fv(u_model, 1, GL_TRUE,  matrix_in_row_major_order );

Now, only this following usage works:

// WebGL API 
gl.uniformMatrix4fv(u_model, false, matrix_in_column_major_order); 

 // OpenGL ES - API 
 // Where u_model is a shader uniform variable 
 glUniformMatrix4fv(u_model, 1, GL_FALSE,  matrix_in_column_major_order);

Further reading:

• In-place matrix transposition
• Row- and column-major order - Wikipedia
• Row major vs. column major, row vectors vs. column vectors | The ryg blog
• Geometry (Row Major vs Column Major Vector)
• Matlab - Row-Major and Column-Major Array Layouts
• Row major versus Column major layout of matrices
• Matrices, Handedness, Pre and Post Multiplication, Row vs Column Major, and Notations
• Memory layout of multi-dimensional arrays - Eli Bendersky's website
• Writing Efficient Code

Row-major matrices

Let the array T of type float[16], containing 16 32 bits floating point numbers, be a comumn-major matrix represeting some affine transform or projection transform.

T: float[16] = [ T[0], T[1], T[2], T[3], T[4], T[5], T[6], T[7], T[8], T[9], T[11], T[12], T[13], T[14], T[15] ]

The array T float[16], represents the following matrix in row-major order:

     /                           \
     |  T[00] T[01] T[02]  T[03]  |   
     |  T[04] T[05] T[06]  T[07]  | 
 T = |  T[08] T[09] T[10]  T[11]  |
     |  T[12] T[13] T[14]  T[15]  |
     \                            /

     /                             \       
     |  A(1,1) A(1,2) A(1,3) A(1,4) |  
     |  A(2,1) A(2,2) A(2,3) A(2,4) |
 T = |  A(3,1) A(3,2) A(3,3) A(3,4) |
     |  A(4,1) A(4,2) A(4,3) A(4,4) |
     \                             /

It follows that: 

  T[00] = A(1,1) ; T[01] = A(1,2) ; T[02] = A(1,3) ; T[03] = A(1,4)
  T[04] = A(2,1) ; T[05] = A(2,2) ; T[06] = A(2,3) ; T[07] = A(2,4)
  T[08] = A(3,1) ; T[09] = A(3,2) ; T[10] = A(3,3) ; T[11] = A(3,4)
  T[12] = A(4,1) ; T[13] = A(4,2) ; T[14] = A(4,3) ; T[15] = A(4,4)


 T: float[16] = [  A(1,1), A(1,2), A(1,3), A(1,4)
                 , A(2,1), A(2,2), A(2,3), A(2,4)
                 , A(3,1), A(3,2), A(3,3), A(3,4)
                 , A(4,1), A(4,2), A(4,3), A(4,4)
               ];

Column-major matrices

Let the array T of type float[16], containing 16 32 bits floating point numbers, be a comumn-major matrix represeting some affine transform or projection transform.

T: float[16] = [ T[0], T[1], T[2], T[3], T[4], T[5], T[6], T[7], T[8], T[9], T[11], T[12], T[13], T[14], T[15] ]

The array T represents the following column-major matrix:

             Matrix in Column-Major Order  
       _                                        _ 
      /                                          \
     +                                            +
T  = | COLUMN[0], COLUMN[1], COLUMN[2], COLUMN[3] |
     +                                            +
      \_                                        _/  

    /                           \
    |  T[00] T[04] T[08]  T[12]  |
    |  T[01] T[05] T[09]  T[13]  |
T = |  T[02] T[06] T[10]  T[14]  |
    |  T[03] T[07] T[11]  T[15]  |
    \                            /

  Where:

                   | T[0] |               | T[4] |              | T[8]  |              | T[12] |
                   | T[1] |               | T[5] |              | T[9]  |              | T[13] |
       COLUMN[0] = | T[2] |   COLUMN[1] = | T[6] |  COLUMN[2] = | T[10] |  COLUMN[3] = | T[14] |
                   | T[3] |               | T[7] |              | T[11] |              | T[15] |

By applying the previous idea to this matrix, it is possible to find:

    /                             \       
    |  A(1,1) A(1,2) A(1,3) A(1,4) |  
    |  A(2,1) A(2,2) A(2,3) A(2,4) |
T = |  A(3,1) A(3,2) A(3,3) A(3,4) |
    |  A(4,1) A(4,2) A(4,3) A(4,4) |
    \                             /

T[00] = A(1,1) ; T[01] = A(2,1) ; T[02] = A(3,1) ; T[03] = A(4,1)
T[04] = A(1,2) ; T[05] = A(2,2) ; T[06] = A(3,2) ; T[07] = A(4,2)
T[08] = A(1,3) ; T[09] = A(2,3) ; T[10] = A(3,3) ; T[11] = A(4,3)
T[12] = A(1,4) ; T[13] = A(2,4) ; T[14] = A(3,4) ; T[15] = A(4,4)


T: float[16] = [  A(1,1), A(2,1), A(3,1), A(4,1)
                , A(1,2), A(2,2), A(3,2), A(4,2)
                , A(1,3), A(2,3), A(3,3), A(4,3)
                , A(1,3), A(2,4), A(3,4), A(4,4)
              ];

The array T can be converted to a matrix by using the next formula:

//----------------------------------//
// Matrix using 1 as initial index  //
//----------------------------------//
// Line 1 
A(1, 1) = T[0]  
A(1, 2) = T[4]  
A(1, 3) = T[8] 
A(1, 4) = T[12] 
// Line 2 
A(2, 1) = T[1]   
A(2, 2) = T[5] 
A(2, 3) = T[9] 
A(2, 4) = T[13]
// Line 3 
A(3, 1) = T[2]
A(3, 2) = T[6]
A(3, 3) = T[10]
A(3, 4) = T[14]
// Line 4 
A(4, 1) = T[3]
A(4, 2) = T[7]
A(4, 3) = T[11]
A(4, 4) = T[15]

//----------------------------------//
// Matrix using 0 as initial index  //
//----------------------------------//
// Line 1 
A(0, 0) = T[0]  
A(0, 1) = T[4]  
A(0, 2) = T[8] 
A(0, 3) = T[12] 
// Line 2 
A(1, 0) = T[1]   
A(1, 1) = T[5] 
A(1, 2) = T[9] 
A(1, 3) = T[13]
// Line 3 
A(2, 0) = T[2]
A(2, 1) = T[6]
A(2, 2) = T[10]
A(2, 3) = T[14]
// Line 4 
A(3, 0) = T[3]
A(3, 1) = T[7]
A(3, 2) = T[11]
A(3, 3) = T[15]

Example: Given the translation matrix T

     | 1  0  0  tx |
     | 0  1  0  ty |
 T = | 0  0  1  tz |
     | 0  0  0   1 |

 T[00] = A(1,1) = 1 ; T[01] = A(2,1) = 0 ; T[02] = A(3,1) = 0 ; T[03] = A(4,1) = 0 (column 1)
 T[04] = A(1,2) = 0 ; T[05] = A(2,2) = 1 ; T[06] = A(3,2) = 0 ; T[07] = A(4,2) = 0 (column 2)
 T[08] = A(1,3) = 0 ; T[09] = A(2,3) = 0 ; T[10] = A(3,3) = 1 ; T[11] = A(4,3) = 0 (column 3)
 T[12] = A(1,4) = tx; T[13] = A(2,4) = ty; T[14] = A(3,4) = tz; T[15] = A(4,4) = 1 (column 4)

Hence, 

 T: float[16] = {  1,  0,  0,  0
                 , 0,  1,  0,  0
                 , 0,  0,  1,  0
                 , tx, ty, tz, 1 }; 

 T: float[16] = { 1,  0,  0,  0, 0,  1,  0,  0, 0,  0,  1,  0, tx, ty, tz, 1 }

 gl.uniformMatrix4fv(u_model, false, T);

This matrix in column-major order is represented as:

                                | 1   0   0  0 |
                                | 0   1   0  0 |
T             =  TRANSPOSE(T) = | 0   0   1  0 |
 column-major                   | tx ty  tz  1 |

Quaternions are a generalization of complex numbers, that was proposed and formulated by the irish mathematician and physicist William Rowan Hamilton, in 1843. This mathematical construct provides a convenient notation for expressing rotation and has a broad range of applications, including computer graphics, computer vision, robotics, orbital mechanics, aerospace and flight dynamics.

A quaternion q is defined as:

• Where:
  • w, x, y, z are real numbers.
  • w => is the quaternion's real (scalar) component.
  • [x, y, z] => is the quaternion's vector componet
  • i, j, k are unit vectors (basis elements) that satisfies a set of properties.

The quaternion q can be represented in the following manner:

• Where:
  • \(\vec{v} = [x, y, z]\) is the quaternion vector component.

The quaternion unit vectors or basis elements satifies the following properties:

And also:

Types of Quaternion

• Real quaternion

Real quaternion, is a quaterion which all elements of vector components are zero.

• Pure quaternion

A pure quaternion, has the real part equal to zero. It can represent a vector in the space.

• Unit quaternion

A quaternion which the norm or magnitude is equal to one.

• Identity quaternion

An identity quaternion when multiplied by another quaternion results in no rotation.

In similar way to identity matrices, the following is valid:

Quaternion Conjugate

The conjugate of a quaternion, \(q^*\), is defined as:

Inverse Quaternion \(q^{-1}\)

The inverse quaternion satifies the following equality:

Quaternion Norm

The norm (magnitude) of a quaternion is given by:

Unit quaternion

An unit quaternion \(q_u\) has norm equal to one.

An unit quaternion is equal to its conjugate:

Quaternion Dot Product

A quaternion dot product is similar to vector dot product. The quaternion dot product between two quaternions q1 and q2 is defined as the sum of the product between the quaternions components.

Angle between two quaternion

The angle between two quaternions q1 and q2 is defined as:

Quaternion Sum

The sum of two quaternions q1 and q2 is:

Quaternion Product (a.k.a Hamilton Product)

The product p of two quaternons q1 and q2 can be computed as:

• Where:
  • p is the quaternion product between q1 and q2
  • \(\vec{v_1} \cdot \vec{v_2}\) is the scalar product between \(\vec{v_1}\) and \(\vec{v_2}\).
  • \(\vec{v1} \times \vec{v_2}\) is the cross product between \(\vec{v_1}\) and \(\vec{v_2}\).
  • \(p_v\) is the vector component of p
  • \(p_w\) is the scalar/real component of p.

The components of vector p are:

The quaternion product is non commuative, as a result, the following holds:

The quaternion product satifies the following properties:

Quaternion product between pure quaternions

The product between two pure quaternions yields the dot product and cross simultaneously, the dot product as the real part of the operation result and cross product as the vector part of the result.

Where:

• Note: A pure quaternion, is the same as space vector, the real part is zero, w = 0.
• \(p = 0 + x_p \cdot i + y_p \cdot j + z_p \cdot k\)
• \(q = 0 + x_q \cdot i + y_q \cdot j + z_q \cdot k\)
• \(p \cdot q\) is the dot product between the vector parts of p and q.
• \(p \times q\) is the cross product (vector product) between p and q.

Quaternion division (MATH432 - page 11)

The division between two quaternions q1 and q2 is defined as:

Quaternion product using matrices representation

The product of two quaternions q1 and q2 can be computed in a more convenient and faster way by using matrices and column vectors multiplication, instead of dot product and cross product.

Let \([q]_v\) be the vector column representation of a quaternion:

Let \([q]_m\) be the following matrix representation of a quaternion:

The quaternion product between q1 and q2, designated by \([q_r]_v\), can be expressed as a matrix-column vector multiplication, which is best suitable for a computer implementation.

By performing the previous matrix multiplication, the quaternion product becomes:

Rotation Quaternion

A rotation around an axis \(\hat{n}\) (unit vector) by angle \(\theta\) can be encoded by quaternion \(q_r\):

• Where:
  • \(\hat{n}\) is a unit vector.
  • \(\hat{n} = x_n \cdot i + y_n \cdot j + z_n \cdot k\)
  • \(|| \hat{n} || = 1\)

The norm of quaternion \(q_r = 1\), thus the following holds:

The coordinate transformation from a vector in reference frame A to a reference frame B, obtained by the rotation of coordinate frame A around a unit vector \(\hat{n}\) can be computed as:

Where:

• \(q_A\) is the coordinate of a vector in the reference frame A (aka coordinate system, aka coordinate frame).
• \(q_B\) is the coordinate of a the same vector in the reference frame B.
• \(q_A = 0 + x_A \cdot i + y_A \cdot j + z_A \cdot k\)
• \(q_B = 0 + x_B \cdot i + y_B \cdot j + z_B \cdot k\)

The previous rotation coordinate transformation can be determined in a more efficient way by applying the following algorithm. (The-Ryg-Blog)

Where:

• \([ Q ]_w\) - real part of a quaternion Q
• \([ Q ]_v\) - vector part of a quaterion Q
• \(w_r = \cos ( \theta / 2)\) - real part of rotation quaternion \(q_r\)
• \(v_r = \hat{n} \cdot \sin( \theta / 2 )\) - vector part of rotation quaternion \(q_r\)
• \(v_B\) - vector part of quaternion \(q_B\)
• \(v_A\) - Vector part of quaternion \(q_A\)
• \(v_r \times v_a\) - Cross product between those two vectors.

Determining rotation and axis from a normalized quaternion

Supposing that the components (w, x, y, z) of normalized quaternion q are known. The represented rotation angle, \(\theta\) and rotation axis (a.k.a orientation) vector, \(\hat{n} = [x_n \quad y_n \quad z_n]\), can be determined in the following way.

Rotation quaternion (known: x, y, z, w)

Rotation angle \(\theta\)

Rotation axis: \(\hat{n} = x_n \cdot i + y_n \cdot j + z_n \cdot k\)

Where:

• Atan2(Y, X) is the 2-argument tangent function which determines the angle quadrant.
• \(q = w + x \cdot i + y \cdot j + z \cdot k\) - is a unitary quaternion, which components x, y and z are assumed to be known. The quaternion q is presumed to be normalized.
• \(\theta\) is the rotation angle. Assumed to be unknown.
• \(\hat{n}\) - is the rotation axis vector. Also assumed to be unknown, which components are \(x_n\), \(y_n\) and \(z_n\).
  • \(\hat{n} = 0 + x_n \cdot i + y_n + \cdot j + z_n \cdot k\)

Quaternion to rotation matrix conversion

The rotation matrix that transforms coordinates from reference frame A to reference frame B can be computed as: (Moti Ben-Ari - Wizmann), (Jernej Barbic)

Where:

• a, b, c, d are the components of the unit quaternion \(q_r\)
• \(a = \cos( \theta / 2)\)
• \(b = x_n \cdot \sin (\theta / 2)\)
• \(c = y_n \cdot \sin (\theta / 2)\)
• \(d = z_n \cdot \sin (\theta / 2)\)

The matrix \(R_n\) can also be simplified to: (Moti Ben-Ari - Wizmann)

Coordinates from frame A can be turned into coordinates from frame B by performing the following matrix multiplication.

From the previous matrix \(R_n\) it is possible to determine the yaw-pitch-roll (Z-X-Y) Euler's angles. (chrobotics)

Where:

• \(\phi\) - is the yaw angle
• \(\theta\) - is the pitch angle
• \(\psi\) -is the roll angle

Quaternion to homogeneous coordinate rotation matrix

The matrix \(R_n\) can be converted to a homogenous coordinate form, \(T_n\), which is more convenient, since homogenous coordinate allows expressing translations in a similar fashion to rotation transformations, just as vector-matrix multiplication.

Homogeneous coordinates from frame A can be turned into coordinates of frame B by using the following expression:

Table with Useful Normalized Quaternions (Rotation Quaternions) - (Ogre3D Game Engine)

Note:

• \(\sqrt{1/2} = \sqrt{0.5} \approx 0.70710\)
• \(\pi = pi = \approx 3.14159265\)

Further Reading

Quaternions:

• Quaternions - Wikipedia
• Quaternions and spatial rotation
• History of Quaternions
• An algorithm for dividing quaternions [PAPER]
• Understanding Quaternions - ChRobotics
• Understanding Euler Angles - ChRobotics
• Quaternion Calculator - Energid / Actin Advanced Control For Robotics
  • Note: Online calculator for Quaternions.
• A Tutorial on Euler Angles and Quaternions - Weizmann Institute of Science
• Quaternion Rotation Primer / Ogre3D Game Engine
• https://web.cs.ucdavis.edu/~amenta/3dphoto/quaternion.pdf
• quatinterp() - Quaternion interpolation between two quaternions - Matlab Aerospace Toolbox
• Maths - Quaternion Interpolation (SLERP) - Euclidianspace.com
• Understanding Quaternions - 3DGEP - 3D Game Engine Programming
• Quaternion - in computer animation
• Quaternion interpolation
• [PDF] Spherical Skinning withDual-Quaternions and QTangents - Spherical Skinning withDual-Quaternions and QTangentsIvo Zoltan FreyCrytek R&D
• [PDF] Orientation & Quaternions - University of Texas
• [PDF] Quaternions: From Classical Mechanics to Computer Graphics, and Beyond
• [PDF] QUATERNION-BASED AUTOPILOT FOR DODECACOPTERS - PART I
• [PDF] Introduction into quaternions for spacecraft attitude representation
• https://github.com/mattlisle/quaternion-ukf (Kalman-Filter and Quaternion)
  • "This implementation of a UKF for tracking orientation of a drone with gyro and accelerometer data follows closely that described in the paper "A Quaternion-based Unscented Kalman Filter for Orientation Tracking" by Edgar Kraft."
• https://stackoverflow.com/questions/tagged/quaternions
• https://stackoverflow.com/questions/tagged/gyroscope

Dual Quaternion:

• Dual quaternions are a combination of quaternions and dual numbers. Unlike ordinary quaternions, which can only express rotation, dual quaternions can express both rotation and translation in a more concise and efficient way than transformation matrices. Some applications of this mathematical construct are: computer graphics; robotics; aerospace engineering and more.
• Dual Quaternions
• Feature guide » 2D and 3D transformations - Magnum Engine
• [PDF] - Paper: A Beginners Guide to Dual-Quaternions - What They Are, How They Work, and How to Use Them for 3D Character Hierarchies
• [PDF] - Paper: Modeling of Spacecraft-Mounted Robot Dynamics and Control Using Dual Quaternions
• [PDF] - Paper: Dual Quaternions for Rigid Transformation Blending
• [PDF] - Report: Dual Quaternions - Yan-Bin Jia
• Matalab File Exchange - Dual quaternion toolbox

Quaternion Computer Implementations:

• t::Quaternion - ROS Framework (C++)
• Qt5 C++ Framework - QQuaternion class (C++)
• Qt5 - Rotation Example
• MinGfx C++ - OpenGL Wrapper (C++)
• FQuat - Unreal Game Engine (C++)
  • Brief: "Floating point quaternion that can represent a rotation about an axis in 3-D space."
• FTransform - Unreal Game Engine (C++)
  • Brief: "Transform composed of Scale, Rotation (as a quaternion), and Translation."
  • Note: It has three member variables: Rotation (class FQaut), Scale3D (FVector) and Translation (FVector).
• FRotator - Unreal Gamer Engine (C++)
  • Brief: "Implements a container for rotation information. All rotation values are stored in degrees."
• Unity3D - Quaternion (C#, CSharp - Unity Game Engine)
  • Brief: "Quaternions are used to represent rotations. They are compact, don't suffer from gimbal lock and can easily be interpolated. Unity internally uses Quaternions to represent all rotations. They are based on complex numbers and are not easy to understand intuitively. You almost never access or modify individual Quaternion components (x,y,z,w); most often you would just take existing rotations (e.g. from the Transform) and use them to construct new rotations (e.g. to smoothly interpolate between two rotations). The Quaternion functions that you use 99% of the time are: Quaternion.LookRotation, Quaternion.Angle, Quaternion.Euler, Quaternion.Slerp, Quaternion.FromToRotation, and Quaternion.identity. (The other functions are only for exotic uses.)"
• Magnum Graphics Library - Quaternion (C++)
• Quaternion Struct - Dotnet (C#)
• Octave - Quaternion Library (OCtave Language, akin to Matlab)
• Blender - MathUtils module

Quaternions (non-dual quaternions) can only encode and represent rotations around some orientation, which is not enough to specify position and scaling. As a result, quaternions rotation matrices need to be concatenated with translation and scaling transform matrices for computing an object's model transform matrix, which is often passed to a shader program via uniform variables.

Assuming that a rendered object has a transform member variable of type RendererTransform containing a a quaternion, for representing rotation; a scale vector and a position vector. The method get_transform() compute the rendered object transform matrix, which is the concatenation of rotation (quaternion), scaling and translation transforms. The model transform determined by get_transform() method can be passed to a shader uniform variable.

Note: Dual-quaternions, which are a combination of dual-numbers and ordinary quaternions, can express both rotation and translation.

• Class RendererTransform:

struct RenderTransform
{
      // Rotation quaternion
      // Set to identity quaternion by default.
      // Members. The quaternion has the fields (w, x, y, z)
      quaternion rotation {1.0, 0.0, 0.0};

      // Scale vector containing scales (sx, sy, sz) of axis (X, Y, Z)
      vec3 scale(1.0, 1.0, 1.0);

      // Contains object's position
      vec3 position(tx, ty, tz);

      // ??? Get overall (end result) transform matrix that
      // will be passed to shader uniform variable.
      mat4 get_transform() const { ????? }

      void translate(float dx, float dy, float dz)
      {
          this->position.x += dx;
          this->position.y += dy;
          this->position.z += dz;
      }

      void set_position(float x, float y, float z)
      {
          this->position.x = x;
          this->position.y = y;
          this->position.z = z;
      }

      void set_scale(float sx float sy, float sz)
      {
         this->scale.x = sx;
         this->scale.y = sy;
         this->scale.z = sz;
      }
};

The transforms are often applied in the following sequence, first scaling, then rotation and finally translation.

Where:

• T is the total transform matrix, which is the product between intermediate transforms and later can be passed to a shader uniform variable.
• \(T_s\) - is the scale transform, determined from the scale vector.
• \(T_t\) - is the translation transform, computed from the position vector.
• \(T_r\) - is the quaternion rotation transform.

Intermediate transforms:

• Translation affine transform form position vector.

• Scaling affine transform from scale vector.

• Rotation affine transform obtainted from rotation quaternion.
  • Note: most linear algebrar or computer graphics math libraries provide subroutines (aka functions) or class methods for converting quaternions to 4x4 rotation matrices. Those libraries also provide subroutines for determining a quaternion's yaw, pitch and roll angles and constructing the quaternion from yaw, pitch and roll angles.

The end transform T is determined by concatenating all intermediate transforms:

Finally, the model transform becomes:

Then, the algorithm for determing the object model matrix, designated by the transform T, could be written as the following pseudo-code:

mat4 RenderTransform::get_transform() const
{
    // Convert quaternion to a 4x4 rotation matrix (affine transform)
    mat4 model_transform = quaternion_to_rotation_transform( this->quaternion );

    // Multiply all elements of first column by X axis scale
    model_transform.comlumn[0] = scale.x * model_transform.column[0];
    // Multiply all elemnts of second column by Y axis scale
    model_transform.comlumn[1] = scale.y * model_transform.column[1];
    // Multiply all elemnts of third column by Z axis scale
    model_transform.comlumn[2] = scale.z * model_transform.column[2];

    // Set all elements of forth column as the 4x1  column vector
    // column4 = [position.x position.y position.z 1.0 ]
    model_transform.column[4] = vec4( position.x, position.y, position.z, 1.0f );

    return model_transform;
}

To rotate the rendered object, it is only necessary to multiply the quaternion variable, called 'rotation' by another quaternion. So, a C++-like pseudo-code for setting or adding a new rotation could be written as:

RenderTransform object_transform;

// Set a new rotation - overriding the old one.
void RenderTransform::set_rotation(float angle_radians, vec3 axis_direction)
{
    // Normalize axis
    vec3 axis = axis_direction.normalize();
    // Make explicit assumption the assumption that the quaternion is normalized.
    assert( norm(axis) ~= 1.0  );
    this->object_transform.quaternion = make_quaternion_from_AngleAxis(angle, axis);
}

// Add a new rotation, rotate the previous quaternion.
void RenderTransform::rotate(float angle_radians, vec3 axis_direction)
{
    vec3 axis = axis_direction.normalize();
    assert( norm(axis) ~= 1.0  );
    quaternon q = make_quaternion_from_AngleAxis(angle, axis);
    this->object_transform.quaternion = this->object_transform.quaternion * q;
}


// Set a new Yaw-Pitch-Roll Euler's angles rotation.
void RenderTransform::set_rotation_YRP(float yaw, float pitch, float roll)
{
    vec3 axis = axis_direction.normalize();
    assert( norm(axis) ~= 1.0  );
    quaternon q = make_quaternion_from_YAW_PITCH_ROLL(yaw, pitch, roll);
    assert( norm(q) ~= 1.0 );
    this->object_transform.quaternion = q;
}

See also:

• Performance of quaternions in the GPU - Metail Tech
• c++ - Matrix multiply with position, quaternion and scale components - Stack Overflow
• Rotating Objects Using Quaternions - Noick Bobic / Gamasutra

It is worth prototyping quaternion math in a numerically-oriented programming language such as Julia Language is helpful for understanding and validating the inner works of quaternion computations. This implementation of quaternion math represents quaternion as an array of four elements and uses a matrix representation for implemeting quaternion multiplication.

# Struct for storing quaternion
# q = w + x·i + y·j + z·k
struct quatp
   arr::Array{Float64,1}
end

function quat(x, y, z, w)
  return quatp([x ; y; z; w])
end

# Turns a quaternion into 4x4 matrix representation
function quat_to_mat4x4(q::quatp)
    w, x, y, z  = q.arr
    mat = [  w -x -y -z
           ; x  w -z  y
           ; y  z  w -x
           ; z -y x w ]
    return mat
 end

function Base.getproperty(q::quatp, sym::Symbol)
    if sym == :arr
        return getfield(q, :arr)
    elseif sym == :w
        return getfield(q, :arr)[1]
    elseif sym == :x
        return getfield(q, :arr)[2]
    elseif sym == :y
        return getfield(q, :arr)[3]
    elseif sym == :z
        return getfield(q, :arr)[4]
    elseif sym == :v
        # Return only the quaternion (vector part)
        return getfield(q, :arr)[2:end]
    end
end

# Norm - quaternion magnitude, akin to vector norm
function norm(q::quatp)
   w, x, y, z = q.arr
   return sqrt(w * w + x * x + y * y + z * z)
end

# Conjugate of a quaternion
function conj(q::quatp)
  return quat(q.w, -q.x, -q.y, -q.z)
end

# Inverse of quaternion q^-1
function Base.:inv(q::quatp)
   w, x, y, z = q.arr
   norm_square =  w * w +  x * x + y * y + z * z
   return  quat(w, -x, -y, -z) / norm_square
end

# Quaternion sum
Base.:+(q1::quatp, q2::quatp) = quat(q1.w + q2.w, q1.x + q2.x, q1.y + q2.y, q1.z + q2.z)
Base.:-(q1::quatp, q2::quatp) = quat(q1.w - q2.w, q1.x - q2.x, q1.y - q2.y, q1.z - q2.z)

Base.:+(q::quatp,   w::Float64) = quat(w + q.w, q.x, q.y, q.z)
Base.:+(w::Float64, q::quatp  ) = quat(w + q.w, q.x, q.y, q.z)
Base.:-(q::quatp,   w::Float64) = quat(q.w - w, q.x, q.y, q.z)

# Mutliplication => scalar X quaternion
Base.:*(a::Float64, q::quatp) = quat( a * q.w, a * q.x, a * q.y, a * q.z)
Base.:*(q::quatp, a::Float64) = quat( a * q.w, a * q.x, a * q.y, a * q.z)

Base.:/(q::quatp, a::Float64) = quat( q.w / a, q.x / a, q.y / a, q.z / a)

# Quaternion multiplication / product
Base.:*(q1::quatp, q2::quatp) = begin
     A = quat_to_mat4x4(q1)
     return quatp( A * q2.arr )
end

# Unit quaternions i, j, k
const i = quat(0.0, 1.0, 0.0, 0.0)
const j = quat(0.0, 0.0, 1.0, 0.0)
const k = quat(0.0, 0.0, 0.0, 1.0)


# Obtain a rotation quaternion that represents a rotation
# by some angle around some direction (normalized vector 'vec')
# => The angle is given in degrees.
function rot_quat(angle::Float64, vec)
   # Normalize vector vec and store the result in v
   xv, yv, zv = vec
   x, y, z = vec / sqrt( xv * xv + yv * yv + zv * zv )
   S = cosd(angle / 2)
   C = sind(angle / 2)
   return quat(C, S * x, S * y, S * z)
end

## Vector cross product VA x VB of two column vectors (3 x 1)
function cross(va, vb)
     x, y, z = va
     T = [ 0 -z y ; z 0 -x; -y x 0]
     T * vb
 end

## Apply a rotation quaternion to a vector (column vector 3 x 1)
##  =>> Vector 3 rows and 1 column
## Algorithm from:
#     =>> https://fgiesen.wordpress.com/2019/02/09/rotating-a-single-vector-using-a-quaternion/
function apply_rot(qr::quatp, va)
   wr = qr.w # Scalar/real part of rotation quaternion
   vr = qr.v # Vector part of rotation quaternion
   vt = 2 * cross(vr, va)
   vb = va + wr * vt + cross(vr, vt)
   return vb
end

# Obtain a rotation matrix - that represents rotation
# around an arbitrary axis (vector)
# Note: The angle is should be in degrees.
#
function make_rotmat(angle::Float64, vec)
   # Normalize vector vec and store components in x, y, z
   x, y, z = vec
   x, y, z = vec / sqrt( x * x + y * y + z * z )
   # sin and cos caching
   C = cosd(angle / 2)
   S = sind(angle / 2)
   # Parameters a, b, c, d => Components of the rotation
   # quaternion
   a = C
   b = x * S
   c = y * S
   d = z * S
   # rotation matrix
   mat = [
        ( a^2 + b^2 - c^2 - d^2 )  ( 2 * (b * c - a * d )  ) ( 2 * (b * d + a * c)   )  0
     ;  ( 2 * (b * c + a * d)   )  ( a^2 - b^2 + c^2 - d^2 ) ( 2 * (c * d - a * b)   )  0
     ;  ( 2 * (b * d - a * c)   )  ( 2 * (c * d + a * b)   ) ( a^2 - b^2 - c^2 + d^2 )  0
     ;  0                          0                         0                          1
   ]
   return mat
end

Quaternion components:

# Construct a quaternion

julia> q = quat(10.5, -25.1, 4.5, 5.6)
quatp([10.5, -25.1, 4.5, 5.6])

julia> q.w # Real part
10.5

julia> q.x # Imaginary part - i axis component
-25.1

julia> q.y # Imaginary part - j axis component
4.5

julia> q.z # Imaginary part - k axis component (Z axis)
5.6

julia> q.arr # Internal array data representation
4-element Array{Float64,1}:
  10.5
 -25.1
   4.5
   5.6

julia> q.v # Vector component (i, j, k) components only
3-element Array{Float64,1}:
 -25.1
   4.5
   5.6

Quaternion conjugate and inverse:

julia> q = quat(1.0, 2.0, 3.0, 4.0)
quatp([1.0, 2.0, 3.0, 4.0])

julia> inv(q)
quatp([0.03333333333333333, -0.06666666666666667, -0.1, -0.13333333333333333])

# Expected (-1.0)
julia> q * inv(q)
quatp([1.0, 0.0, 0.0, 0.0])

# Expected (-1.0)
julia> inv(q) * q
quatp([1.0, -5.551115123125783e-17, -6.938893903907228e-18, 5.551115123125783e-17])

Normalize a quaternion:

# Expected normalized to be: (w = 0.18257, x = 0.36515, y = 0.54772, z = 0.7303)
# Test case from: Matlab Robotics Toolbox - https://www.mathworks.com/help/robotics/ref/quaternion.normalize.html

 julia> q = quat(1.0, 2.0, 3.0, 4.0)
 quatp([1.0, 2.0, 3.0, 4.0])

 julia> qn = normalize(q)
 quatp([0.18257418583505536, 0.3651483716701107, 0.5477225575051661, 0.7302967433402214])

 julia> norm(qn)
 0.9999999999999999

Test case 1 - for quaternion multiplication:

• Test quaternion laws:
• Expected: i * i = j * j = k * k = -1
• Expected: i * j = k
• Expected: j * k = i
• Expected: k * i = j
• Expected: j * i = -k

# ---- Show components i, j, k --------#

julia> i
quatp([0.0, 1.0, 0.0, 0.0])

julia> j
quatp([0.0, 0.0, 1.0, 0.0])

julia> k
quatp([0.0, 0.0, 0.0, 1.0])

# ----- Check multiplication of i, j, k by themselves -----------#
#

julia> i * i # Expected -1 or (-1) + 0·i + 0·j + 0·k
quatp([-1.0, 0.0, 0.0, 0.0])

julia> j * j # Expected -1
quatp([-1.0, 0.0, 0.0, 0.0])

julia> k * k # Expected -1
quatp([-1.0, 0.0, 0.0, 0.0])

#------- Check multiplication i * j, j * k, k * i and j * i
#
julia> i * j # Expected k
quatp([0.0, 0.0, 0.0, 1.0])

julia> j * i # Expected -k
quatp([0.0, 0.0, 0.0, -1.0])

julia> k * i # Expected j
quatp([0.0, 0.0, 1.0, 0.0])

julia> i * k # Expected -j
quatp([0.0, 0.0, -1.0, 0.0])

Test case 2 - for quaternion multiplication:

• From: Matlab Aerospace Toolbox
• q1 = [ 1.0 0.0 1.0 0.00 ]
• q2 = [ 1.0 0.5 0.5 0.75 ]
• Expected result =>> q1 * q2 = [ 0.5 1.25 1.5 0.25 ]

# -------- Test A ------------------------#
julia> q1 = quat(1.0, 0.0, 1.0, 0.0)
quatp([1.0, 0.0, 1.0, 0.0])

julia> q2 = quat(1.0, 0.5, 0.5, 0.75)
quatp([1.0, 0.5, 0.5, 0.75])

# ------- Test B ---------------------#

julia> q1 = 1.0 + 0.0i + 1.0j + 0.0k
quatp([1.0, 0.0, 1.0, 0.0])

julia> q2 = 1.0 + 0.5i + 0.5j + 0.75k
quatp([1.0, 0.5, 0.5, 0.75])

julia> q1 * q2
quatp([0.5, 1.25, 1.5, 0.25])

Test case 3 - for quaternion multiplication:

• From: Maplesoft - Quaternion
• q1 = ( -6) + (-8)i + (-5)j + 7k
• q2 = (-14) + (-12)i + (-19)j + (-17)k
• Expected: q1 * q2 = 12 + 402i + (-36)j + 96k
• Expected: q2 * q1 = 12 + (-34)i + 404j + (-88)k
• Expected: q2 + q1 = (-20) + (-20)i + (-24)j + (-10)k

julia> q1 = -6.0 + -8.0i + -5.0j + 7.0k
quatp([-6.0, -8.0, -5.0, 7.0])

julia> q2 = -14.0 + -12.0i + -19.0j + -17.0k
quatp([-14.0, -12.0, -19.0, -17.0])

julia> q1 * q2
quatp([12.0, 402.0, -36.0, 96.0])

julia> q2 * q1
quatp([12.0, -34.0, 404.0, -88.0])

julia> q1 + q2
quatp([-20.0, -20.0, -24.0, -10.0])

julia> q2 + q1
quatp([-20.0, -20.0, -24.0, -10.0])

Test algorithm - make_rotmat that computes rotation matrices around an arbitrary axis:

Define functions for computing rotation matrices:

# Rotation around X axis
function rot_x(angle::Float64)
   C = cosd(angle)
   S = sind(angle)
   return [  1   0   0  0
           ; 0   C  -S  0
           ; 0   S   C  0
           ; 0   0   0  1
          ]
end

# Rotation around X axis
function rot_y(angle::Float64)
   c = cosd(angle)
   s = sind(angle)
   return [  c   0   s  0
           ; 0   1   0  0
           ;-s   0   c  0
           ; 0   0   0  1
          ]
end

# Returns homogeneous rotation matrix around Z axis
# => The angle must be in degrees, not radians
function rot_z(angle::Float64)
   c = cosd(angle)
   s = sind(angle)
   return [  c  -s  0  0
           ; s   c  0  0
           ; 0   0  1  0
           ; 0   0  0  1
          ]
end

Test rotation matrices:

const axis_x = [1.0 0.0 0.0]';
const axis_y = [0.0 1.0 0.0]';
const axis_z = [0.0 0.0 1.0]';

# -------- Rotation around Z axis ---------- #
#

julia> rot_z(90.0)
4×4 Array{Float64,2}:
 0.0  -1.0  0.0  0.0
 1.0   0.0  0.0  0.0
 0.0   0.0  1.0  0.0
 0.0   0.0  0.0  1.0

julia> make_rotmat(90.0, axis_z)
4×4 Array{Float64,2}:
 0.0  -1.0  0.0  0.0
 1.0   0.0  0.0  0.0
 0.0   0.0  1.0  0.0
 0.0   0.0  0.0  1.0


julia> rot_z(45.0)
4×4 Array{Float64,2}:
 0.707107  -0.707107  0.0  0.0
 0.707107   0.707107  0.0  0.0
 0.0        0.0       1.0  0.0
 0.0        0.0       0.0  1.0

julia> make_rotmat(45.0, axis_z)
]4×4 Array{Float64,2}:
 0.707107  -0.707107  0.0  0.0
 0.707107   0.707107  0.0  0.0
 0.0        0.0       1.0  0.0
 0.0        0.0       0.0  1.0

julia> rot_z(125.0)
4×4 Array{Float64,2}:
 -0.573576  -0.819152  0.0  0.0
  0.819152  -0.573576  0.0  0.0
  0.0        0.0       1.0  0.0
  0.0        0.0       0.0  1.0

julia> make_rotmat(125.0, axis_z)
4×4 Array{Float64,2}:
 -0.573576  -0.819152  0.0  0.0
  0.819152  -0.573576  0.0  0.0
  0.0        0.0       1.0  0.0
  0.0        0.0       0.0  1.0


# ----- Rotation around X axis ------------#
#
julia> rot_x(0.0)
4×4 Array{Float64,2}:
 1.0  0.0   0.0  0.0
 0.0  1.0  -0.0  0.0
 0.0  0.0   1.0  0.0
 0.0  0.0   0.0  1.0

julia> make_rotmat(0.0, axis_x)
4×4 Array{Float64,2}:
 1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0


julia> rot_x(45.0)
4×4 Array{Float64,2}:
 1.0  0.0        0.0       0.0
 0.0  0.707107  -0.707107  0.0
 0.0  0.707107   0.707107  0.0
 0.0  0.0        0.0       1.0

julia> make_rotmat(45.0, axis_x)
4×4 Array{Float64,2}:
 1.0  0.0        0.0       0.0
 0.0  0.707107  -0.707107  0.0
 0.0  0.707107   0.707107  0.0
 0.0  0.0        0.0       1.0


julia> rot_x(90.0)
4×4 Array{Float64,2}:
 1.0  0.0   0.0  0.0
 0.0  0.0  -1.0  0.0
 0.0  1.0   0.0  0.0
 0.0  0.0   0.0  1.0

julia> make_rotmat(90.0, axis_x)
4×4 Array{Float64,2}:
 1.0  0.0   0.0  0.0
 0.0  0.0  -1.0  0.0
 0.0  1.0   0.0  0.0
 0.0  0.0   0.0  1.0

Test rotation quaternions:

# Rotation around X axis
function rot_xx(angle::Float64)
   C = cosd(angle)
   S = sind(angle)
   return [ 1  0  0  ; 0   C  -S   ; 0   S   C     ]
end

julia> va = [6.0 10.5 2.56]'
3×1 LinearAlgebra.Adjoint{Float64,Array{Float64,2}}:
  6.0
 10.5
  2.56

angle = 60.0 # 60 degrees

julia> vb_rotmat = rot_xx(angle) * va
3×1 Array{Float64,2}:
  6.0
  3.032974966311837
 10.373266739736605

# Rotation quaternion that represents rotation of 60 degrees around X axis.
julia> qr = rot_quat(60.0, [1.0 0.0 0.0]')
quatp([0.8660254037844386, 0.5, 0.0, 0.0])

# Apply transform to vector va
julia> vb_qr = apply_rot(qr, va)
3×1 Array{Float64,2}:
  6.0
  3.032974966311837
 10.373266739736607

The library GLM (OpenGL math library) contains many classes and subroutines for computer graphics computations, such as: homogeneous coordinates; quaternios; 1D, 2D, 3D and homogeneous coordinates vectors; vector-matrix operations and so on. Aside those facilities, the library also provides the subroutines glm::lookAt(), for computing view matrix transform, that turns world coordinates into camera coordinates; glm::perspective() - for computing the projection matrix, that turns camera coordinates into clip-space coordinates (NDC coordinates with range -1.0 to 1.0) and also glm::ortho for computing the orthogonal perspective matrix. The GLM library is designed with syntax based on GLSL (OpenGL shading language).

Web Site:

• OpenGL Mathematics

Repository:

• https://github.com/g-truc/glm

Other Documentations:

• https://openframeworks.cc/documentation/glm/

Type signature of most relevant GLM functions:

// Converts an input angle in degrees to radians.
//
//  angle_radians = angle_degrees * (PI / 180.0 )
//
float glm::radians(foat  degrees);

// ------ Basic Matrix Transformation =>> Useful for model matrix ---------------//
glm::mat4 glm::rotate   ( glm::mat4 const & m, float angle, glm::vec3 const & axis );
glm::mat4 glm::scale    ( glm::mat4 const & m, glm::vec3 const & factors           );
glm::mat4 glm::translate( glm::mat4 const & m, glm::vec3 const & translation       );


// ---- Camera View matrix and Camera's matrix -------------//
glm::mat4 glm::lookAt( glm::vec3 const & eye, glm::vec3 const & look, glm::vec3 const & up );
glm::mat4 glm::ortho( float left, float right, float bottom, float top, float near, float far );
glm::mat4 glm::ortho( float left, float right, float bottom, float top );
glm::mat4 glm::frustum( float left, float right, float bottom, float top, float near, float far );
glm::mat4 glm::perspective( float fovy, float aspect, float near, float far);

Subroutines for vector computations

• glm::normalize()
  • => Normalize a vector.
• glm::dot(vec1, vec2)
  • => Vector dot product.
• glm::length(vec)
  • => Vector norm, aka magnitude.
• glm::distance(vec1, vec2)
  • => Distance between two vectors thatrepresents position.
• glm::cross(vec1, vec2)
  • => Returns the cross product between two vectors.

Subroutine for Camera Implementation

Camera's View Matrix:

• glm::mat4 glm::lookAt (eye, look, up );
  • Brief: Return View matrix transform - which transforms world coordinates to camera's space coordinates.
  • eye => Vector containing the camera's current position in world coordinates.
  • look => Vector (Position) to where camera is looking at.
  • up => Camera orientation. Suitable fault value Y axis or (x = 0, y = 1, z = 0).

Camera's Projection Matrix:

• glm::mat4 glm::perspective ( fovy, float aspect, near, far);
  • Returns projection matrix that transforms camera-space coordinates to clip-space coordinates (NDC - Normalized Device Coordinates within range -1.0 to 1.0). Anything out of the near, far range is not visible.
  • fovy => Field view angle in radians.
  • near => Distance to near plane.
  • far => Distance to far plane.
• glm::mat4 glm::ortho ( left, right, bottom, top );
  • Returns a projection matrix for ortographic projection. This type projection is most suitable for 2D drawing; CAD - Computer Aided Design; Technical drawing view and on.

GLM functions - transform matrix multiplication order

Consider the subroutine, glm::translate() - which translate a coordinate system.

glm::mat4 SOURCE_MATRIX = something();
glm::mat4 OUTPUT_MATRIX =  glm::translate( SOURCE_MATRIX, glm::vec3(DX, DY, DZ));

GLM computes the OUTPUT_MATRIX by performing the matrix multiplication in the following way. The same logic is applicable to the subroutines: glm::rotate, glm::scale and so on.

OUTPUT_MATRIX = SOURCE_MATRIX * Translate_Transform(DX, DY, DZ)

                                | 1  0  0  DX |
OUTPUT_MATRIX = SOURCE_MATRIX * | 0  1  0  DY |
                                | 0  0  1  DZ |
                                | 0  0  0  1  |

Consider the following sequence of operations:

 const glm::vec4 axis_Z = glm::vec3(0, 0, 1);

// Reset model matrix to identity matrix
 glm::mat4 model(1.0);

 // Move to (X, Y) position
 // model = model * T_translate(_x, y, 0.0)
 // model = identity * T_translate(_x, y, 0.0) = T_translate(_x, y, 0.0)
 model = glm::translate( model, glm::vec3(_x, _y, 0.0)  );

 // Scale object (increase or decrease object size)
 // model = model * T_scale
 // model = T_translate * T_scale
 model =  glm::scale( model, glm::vec3(_scale, _scale, _scale) );

 // Rotate from a given angle around Z axis at current object X, Y  postion
 // model = model * T_rotation
 // model =  (T_translate * T_scale) * T_rot
 model = glm::rotate( model, glm::radians(_angle),  axis_Z);

The final value of the model matrix, is product between the following affine transforms:

model  =  Identity * T_translate * T_scale * T_rotation
model  =  T_translate * T_scale * T_rotation

Download Library Source

$ mkdir -p /tmp/temp && cd /temp

# Download source code archive
$ >> curl -o glm.zip -L https://github.com/g-truc/glm/archive/master.zip

# Extract code
$ >> unzip  glm.zip

# Enter the in the extracted directory
$ >> cd glm-master/

# List directory content
$ >> ls
cmake/  CMakeLists.txt  copying.txt  doc/  glm/  manual.md  readme.md  test/  util/

Load the library in CERN's ROOT repl

 $ >> ~/Applications/root/bin/root
ERROR in cling::CIFactory::createCI(): cannot extract standard library include paths!
Invoking:
  LC_ALL=C ccache  -O2 -DNDEBUG -xc++ -E -v /dev/null 2>&1 | sed -n -e '/^.include/,${' -e '/^ \/.*++/p' -e '}'
Results was:
With exit code 0
   ------------------------------------------------------------------
  | Welcome to ROOT 6.22/02                        https://root.cern |
  | (c) 1995-2020, The ROOT Team; conception: R. Brun, F. Rademakers |
  | Built for linuxx8664gcc on Aug 17 2020, 12:46:52                 |
  | From tags/v6-22-02@v6-22-02                                      |
  | Try '.help', '.demo', '.license', '.credits', '.quit'/'.q'       |
   ------------------------------------------------------------------

root [0]

root [0] .I .

root [0] .I .
root [1] #include <glm/glm.hpp>
root [2] #include <glm/gtc/matrix_transform.hpp>
root [3] #include <glm/gtc/type_ptr.hpp>
root [4] #include <glm/gtx/string_cast.hpp>

// Note: GLM matrices are stored in Column-major order
void show_matrix(const char* label, glm::mat4 const& m){
    std::cout << "\n [MATRIX] " << label << " = " << '\n';
    std::cout << std::fixed << std::setprecision(3);
    for(size_t i = 0; i < 4; i++)
    {
        for(size_t j = 0; j < 4; j++)
        {
            std::cout << std::setw(8) << m[j][i];
        }
        std::cout << '\n';
    }
}

Construct a matrix from column vectors:

auto col0 = glm::vec4(25, 40, 7, 0)
auto col1 = glm::vec4(8, 9, 10, 0)
auto col2 = glm::vec4(-5, 6, 9, 0)
auto col3 = glm::vec4(100, 50, 20, 1)

root [9] auto M = glm::mat4(col0, col1, col2, col3)

root [24] show_matrix("M", M)

 [MATRIX] M =
  25.000   8.000  -5.000 100.000
  40.000   9.000   6.000  50.000
   7.000  10.000   9.000  20.000
   0.000   0.000   0.000   1.000

This matrix could also be constructed in the following mode:

 auto MM = glm::mat4(25, 40, 7, 0  ,8, 9, 10, 0  ,-5, 6, 9, 0 ,100, 50, 20, 1)

root [27] show_matrix("MM", MM)

 [MATRIX] MM =
  25.000   8.000  -5.000 100.000
  40.000   9.000   6.000  50.000
   7.000  10.000   9.000  20.000
   0.000   0.000   0.000   1.000

// Note: It shows the matrix columns, not the matix rows.
root [39] glm::to_string(M)
(std::string) "mat4x4((25.000000, 40.000000, 7.000000, 0.000000), (8.000000, 9.000000, 10.000000, 0.000000), (-5.000000, 6.000000, 9.000000, 0.000000), (100.000000, 50.000000, 20.000000, 1.000000))"
root [40]

Extract columns from matrix MM:

#include <glm/gtc/matrix_access.hpp>

auto col0 = glm::column(MM, 0)
auto col1 = glm::column(MM, 1)
auto col2 = glm::column(MM, 2)
auto col3 = glm::column(MM, 3)

root [58] glm::to_string(col0)
(std::string) "vec4(25.000000, 40.000000, 7.000000, 0.000000)"

root [60] glm::to_string(col1)
(std::string) "vec4(8.000000, 9.000000, 10.000000, 0.000000)"

root [62] glm::to_string(col2)
(std::string) "vec4(-5.000000, 6.000000, 9.000000, 0.000000)"

root [64] glm::to_string(col3)
(std::string) "vec4(100.000000, 50.000000, 20.000000, 1.000000)"

// Get elements from column 0 (col0)
root [65] col0[0]
(float) 25.0000f

root [66] col0[1]
(float) 40.0000f

root [67] col0[2]
(float) 7.00000f

root [68] col0[3]
(float) 0.00000f

Extract rows from matrix MM:

auto row0 = glm::row(MM, 0)
auto row1 = glm::row(MM, 1)
auto row2 = glm::row(MM, 2)
auto row3 = glm::row(MM, 3)

root [73] glm::to_string(row0)
(std::string) "vec4(25.000000, 8.000000, -5.000000, 100.000000)"

root [74] glm::to_string(row1)
(std::string) "vec4(40.000000, 9.000000, 6.000000, 50.000000)"

root [75] glm::to_string(row2)
(std::string) "vec4(7.000000, 10.000000, 9.000000, 20.000000)"

root [76] glm::to_string(row3)
(std::string) "vec4(0.000000, 0.000000, 0.000000, 1.000000)"

Modify matrix columns:

root [77] show_matrix("MM", MM)

 [MATRIX] MM =
  25.000   8.000  -5.000 100.000
  40.000   9.000   6.000  50.000
   7.000  10.000   9.000  20.000
   0.000   0.000   0.000   1.000

root [80] MM[0] = glm::vec4(-90.f, 50.0f, 32.5f, 26.0f)


root [81] show_matrix("MM", MM)

 [MATRIX] MM =
 -90.000   8.000  -5.000 100.000
  50.000   9.000   6.000  50.000
  32.500  10.000   9.000  20.000
  26.000   0.000   0.000   1.000

Obtain transpose matrix:

 root [82] auto tmat = glm::transpose(MM)
 (glm::mat<4, 4, float, glm::qualifier::packed_highp> &) @0x7f7b71f0c20c


 root [83] show_matrix("tmat", tmat)

  [MATRIX] tmat =
  -90.000  50.000  32.500  26.000
    8.000   9.000  10.000   0.000
   -5.000   6.000   9.000   0.000
  100.000  50.000  20.000   1.000

// ----- Access first column elements (column 0) ----//
//  =>> The matrix is column-major order (Fortran-like data layout),
//      not C or C++ data-layout, which is row-major order
//
// Tmat[j column-number][i row-number]

// ----- Column 0 elements ----------//

root [84] tmat[0][0]
(float) -90.0000f

root [85] tmat[1][0]
(float) 50.0000f

root [86] tmat[0][0]
(float) -90.0000f

root [87] tmat[0][1]
(float) 8.00000f

root [88] tmat[0][2]
(float) -5.00000f

root [89] tmat[0][3]
(float) 100.000f

Get pointer to first element:

root [90] float* ptr = glm::value_ptr(tmat)
(float *) 0x7f7b71f0c20c

root [91] ptr[0]
(float) -90.0000f

root [92] ptr[1]
(float) 8.00000f

root [93] ptr[15]
(float) 1.00000f

root [94] ptr[10]
(float) 9.00000f

root [95] *(ptr + 1)
(float) 8.00000f

root [96] *(ptr + 10)
(float) 9.00000f

root [97] *(ptr + 15)
(float) 1.00000f

Null 4x4 matrix:

root [53] auto zero_4x4 = glm::mat4()
(glm::mat<4, 4, float, glm::qualifier::packed_highp> &) @0x7f904e6c20e0

root [54] show_matrix("zero_4x4", zero_4x4)

 [MATRIX] zero_4x4 =
   0.000   0.000   0.000   0.000
   0.000   0.000   0.000   0.000
   0.000   0.000   0.000   0.000
   0.000   0.000   0.000   0.000

Identity matrix:

root [55] auto id_4x4 = glm::mat4(1.0)
(glm::mat<4, 4, float, glm::qualifier::packed_highp> &) @0x7f904e6c2120
root [56]
root [56] show_matrix("id_4x4", id_4x4)

 [MATRIX] id_4x4 =
   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

Matrix translation coordinate transform:

root [58] auto t1 = glm::translate(id_4x4, glm::vec3(2.0, 10.0, 30.0))
(glm::mat<4, 4, float, glm::qualifier::packed_highp> &) @0x7f904e6c2160

root [59] show_matrix("t1", t1)

 [MATRIX] t1 =
   1.000   0.000   0.000   2.000
   0.000   1.000   0.000  10.000
   0.000   0.000   1.000  30.000
   0.000   0.000   0.000   1.000


root [62] auto t2 = glm::translate(t1, glm::vec3(-10, 100.0, 200.0))
(glm::mat<4, 4, float, glm::qualifier::packed_highp> &) @0x7f904e6c21a0

root [63] show_matrix("t2", t2)

 [MATRIX] t2 =
   1.000   0.000   0.000  -8.000
   0.000   1.000   0.000 110.000
   0.000   0.000   1.000 230.000
   0.000   0.000   0.000   1.000

Rotation around Z axis of 90 degrees:

• glm::mat4 glm::rotate(glm::mat4 const& matrix, float angle_radians, glm::vec3 const& axis)
• Rotate around a given axis. The angle is given in radians.

root [65] const auto axis_z = glm::vec3(0.0f, 0.0f, 1.0f);
root [66] const auto axis_y = glm::vec3(0.0f, 1.0f, 0.0f);
root [67] const auto axis_x = glm::vec3(1.0f, 0.0f, 0.0f);


/*          | cos(t)  -sin(t)   0   0  |
 *          | sin(t)   cos(t)   0   0  |
 *  Rz(t) = |  0        0       1   0  |
 *          |  0        0       0   1  |
 *
 *
 *  t_rotZ = Rz(90) x id_4x4 = Rz(90)
 */
root [71] auto t_rotZ = glm::rotate(id_4x4, glm::radians(90.0f), axis_z)
(glm::mat<4, 4, float, glm::qualifier::packed_highp> &) @0x7f904e6c2204

root [72] show_matrix("t_rotZ", t_rotZ)

 [MATRIX] t_rotZ =
  -0.000  -1.000   0.000   0.000
   1.000  -0.000   0.000   0.000
   0.000   0.000   1.000   0.000
   0.000   0.000   0.000   1.000
root [73]

 root [96] glm::to_string(t_rotZ)
 (std::string) "mat4x4((-0.000000, 1.000000, 0.000000, 0.000000), ... "
 root [97]

Scaling transformation:

root [75] auto s1 = glm::scale(id_4x4, glm::vec3(2.0, 2.0, 2.0))
(glm::mat<4, 4, float, glm::qualifier::packed_highp> &) @0x7f904e6c2244

root [76] show_matrix("s1", s1)

 [MATRIX] s1 =
   2.000   0.000   0.000   0.000
   0.000   2.000   0.000   0.000
   0.000   0.000   2.000   0.000
   0.000   0.000   0.000   1.000

root [78] s1 = glm::scale(s1, glm::vec3(5.0, 5.0, 5.0))
(glm::mat &) @0x7f904e6c2244

root [79] show_matrix("s1", s1)

 [MATRIX] s1 =
  10.000   0.000   0.000   0.000
   0.000  10.000   0.000   0.000
   0.000   0.000  10.000   0.000
   0.000   0.000   0.000   1.000
root [80]

// Apply transform to vector:

root [82] auto res = s1 * glm::vec4(2.0, 5.0, 10.0, 1.0)
(glm::vec<4, float, glm::qualifier::packed_highp> &) @0x7f904e6c2284

root [91] glm::to_string(res)
(std::string) "vec4(20.000000, 50.000000, 100.000000, 1.000000)"

root [92] res[0]
(float) 20.0000f

root [93] res[1]
(float) 50.0000f

root [94] res[2]
(float) 100.000f

root [95] res[3]
(float) 1.00000f

Scaling transformation:

root [97] auto s = glm::scale(id_4x4, glm::vec3(4.0, 5.0, 6.0))
(glm::mat<4, 4, float, glm::qualifier::packed_highp> &) @0x7f904e6c22a4

root [98] show_matrix("s", s)

 [MATRIX] s =
   4.000   0.000   0.000   0.000
   0.000   5.000   0.000   0.000
   0.000   0.000   6.000   0.000
   0.000   0.000   0.000   1.000
root [99]
root [99] s = glm::translate(s, glm::vec3(4, -5, 9))
(glm::mat &) @0x7f904e6c22a4


root [100] show_matrix("s", s)

 [MATRIX] s =
   4.000   0.000   0.000  16.000
   0.000   5.000   0.000 -25.000
   0.000   0.000   6.000  54.000
   0.000   0.000   0.000   1.000
root [101]


root [105] s = glm::scale(s, glm::vec3(5.0, 2.0, 3.0))
(glm::mat &) @0x7f904e6c22a4


root [106] show_matrix("s", s)

 [MATRIX] s =
  25.000   0.000   0.000  20.000
   0.000   4.000   0.000 -10.000
   0.000   0.000   9.000  27.000
   0.000   0.000   0.000   1.000

Inverse matrix:

root [111] show_matrix("s", s)

 [MATRIX] s =
  25.000   0.000   0.000  20.000
   0.000   4.000   0.000 -10.000
   0.000   0.000   9.000  27.000
   0.000   0.000   0.000   1.000

root [112] inv_s = glm::inverse(s)
(glm::mat<4, 4, float, glm::qualifier::packed_highp> &) @0x7f904e6c2324

root [113] show_matrix("inv_s", inv_s)

 [MATRIX] inv_s =
   0.040  -0.000   0.000  -0.800
  -0.000   0.250  -0.000   2.500
   0.000  -0.000   0.111  -3.000
  -0.000   0.000  -0.000   1.000
root [114]
root [114]

GLM Quaternions

Create an identity quaternion, which represents no rotation.

root [26] auto q = glm::quat(1.0, 0.0, 0.0, 0.0)
(glm::qua<float, glm::qualifier::packed_highp> &) @0x7fc8b47e6020

root [27] q.x
(float) 0.00000f
root [28]
root [28] q.y
(float) 0.00000f
root [29]
root [29] q.z
(float) 0.00000f
root [30]
root [30] q.w
(float) 1.00000f
root [31]
root [31]

Create a quaternion that represents a rotation of 90 degrees around Z axis:

root [39] auto qrot_z = glm::angleAxis(glm::radians(90.0f), glm::vec3{0.0f, 0.0f, 1.0f} )
(glm::qua<float, glm::qualifier::packed_highp> &) @0x7fc8b47e6050

root [40] glm::to_string(qrot_z)
(std::string) "quat(0.707107, {0.000000, 0.000000, 0.707107})"

root [42] qrot_z.w
(float) 0.707107f

root [43] qrot_z.x
(float) 0.00000f

root [44] qrot_z.y
(float) 0.00000f

root [45] qrot_z.z
(float) 0.707107f

Get quaternion angle:

// Angle in radians.
root [61] glm::angle(qrot_z)
(float) 1.57080f

// Angle in degrees
root [62] glm::angle(qrot_z) * 180.0 / M_PI
(double) 90.000003

Get quaternion rotation axis vector:

root [63] auto axis = glm::axis(qrot_z)
(glm::vec<3, float, glm::qualifier::packed_highp> &) @0x7fc8b47e6150

root [64] glm::to_string(axis)
(std::string) "vec3(0.000000, 0.000000, 1.000000)"

Determine quaternion conjugate:

// Conjugate
root [68] auto res = glm::conjugate(qrot_z)
(glm::qua<float, glm::qualifier::packed_highp> &) @0x7fc8b47e616c

root [69] glm::to_string(res)
(std::string) "quat(0.707107, {-0.000000, -0.000000, -0.707107})"

root [71] glm::to_string( qrot_z * res )
(std::string) "quat(1.000000, {0.000000, 0.000000, 0.000000})"

Determine quaternion inverse:

root [72] glm::to_string( glm::inverse(qrot_z) )
(std::string) "quat(0.707107, {-0.000000, -0.000000, -0.707107})"

Computer quaternion norm (aka magnitude):

root [74] glm::length(qrot_z)
(float) 1.00000f