Introduction

dictk (Digital Image Correlation Toolkit) is a Python library for digital image correlation (DIC) — comparing images of a specimen before and after deformation to measure displacement and strain fields.

Installation

pip install dictk

Overview

DIC compares a reference image (often of an undeformed specimen) to a subject image (often of a deformed specimen) as a means to optically quantify displacement and strain fields.

Before we discuss the image portion of DIC, it is useful to introduce basic concepts of Continuum Mechanics and the Finite Element Method. This introduction will formalize definitions of displacement and strain and lay the groundwork for motion of discrete points locatable in reference/subject image pairs that are mapped to nodes of a finite element mesh (i.e., nodes).

The basic workflow is as follows:

• Collect a pair of before and after images that capture an area of interest.
  • The before image, often called the reference image, captures some (possibly undeformed or deformed) configuration of the specimen.
  • The after image, often called the current image or subject image, captures the specimen in a newly deformed state, different from the before state.
• Within each pair of images, identify points of interest.
• Use the $(X,Y)$ coordinates of the image points as coordinates of nodes that compose a finite element mesh.
• Use DIC to quantify the displacement field of each (nodal) point.
• Use the mesh (nodal) deformation to calculate a discrete strain field at Gauss points.

The following sections explicate this workflow in detail.

Literature Overview

Sutton, Orteu, and Schreier provide the field's most comprehensive reference text: basic image-correlation and computer-vision theory, camera calibration and optics, two-dimensional, stereo, and volumetric correlation methods, and the continuum-mechanics, statistics, and optimization background needed to apply them all.¹ For a shorter, more recent introduction alongside an annotated bibliography of the field, see Brown.² The development of modern Digital Image Correlation (DIC) spans from specialized hardware acquisition to advanced mathematical frameworks and open-source implementations. For high-resolution SEM imaging, Lenthe et al. describe advanced detector signal acquisition and scanning methods,³ while Black et al. demonstrate high-throughput measurements using multi-beam SEM imaging.⁴ When dealing with plastic localization, the Heaviside-DIC approach offers a robust method for capturing discontinuities.⁵

Mathematical foundations for the DICe engine are established through the work of Turner, covering gradient-based local formulations,⁶ stereo correlation and triangulation,⁷ and the implementation of virtual strain gauges.⁸ Alternatively, the Augmented Lagrangian DIC (ALDIC) framework provides a global approach to displacement fields, available as both a theoretical framework⁹ and a MATLAB implementation.¹⁰ For complex geometries and large deformations, the SpatioTemporally Adaptive Quadtree (STAQ) mesh offers significant improvements in resolution.¹¹

Software accessibility is primarily supported by the µDIC toolkit, an open-source Python library for DIC tasks.¹² The project's evolution is tracked through specific version forks,¹³ comprehensive online documentation, and its primary repository.¹⁴ Finally, looking toward performance optimization, the Rust GPU project provides the infrastructure for hardware-accelerated computations that could benefit future DIC processing pipelines.¹⁵

References

1. Sutton MA, Orteu JJ, Schreier HW. Image Correlation for Shape, Motion and Deformation Measurements: Basic Concepts, Theory and Applications. Springer; 2009. ISBN 978-0-387-78746-6. link ↩

2. Brown C. Introduction to Digital Image Correlation (DIC) with annotated bibliography. Lawrence Livermore National Laboratory. 2025 Dec. LLNL-TR-2013494. link ↩

3. Lenthe WC, Stinville JC, Echlin MP, Chen Z, Daly S, Pollock TM. Advanced detector signal acquisition and electron beam scanning for high resolution SEM imaging. Ultramicroscopy. 2018 Dec 1;195:93-100. link ↩

4. Black RL, Garbowski T, Bean C, Eberle AL, Nickell S, Texier D, Valle V, Stinville JC. High-throughput high-resolution digital image correlation measurements by multi-beam SEM imaging. Experimental Mechanics. 2023 Jun;63(5):939-53. link ↩

5. Bourdin F, Stinville JC, Echlin MP, Callahan PG, Lenthe WC, Torbet CJ, Texier D, Bridier F, Cormier J, Villechaise P, Pollock TM. Measurements of plastic localization by heaviside-digital image correlation. Acta Materialia. 2018 Sep 15;157:307-25. link ↩

6. Turner DZ. An overview of the gradient-based local DIC formulation for motion estimation in DICe. 2016 Aug 19. SAND2016-7360R. link ↩

7. Turner DZ. An overview of the stereo correlation and triangulation formulations used in DICe. 2017 Mar 10. SAND2017-1876R. link ↩

8. Turner DZ. An overview of the virtual strain gauge formulation in DICe. 2018 May 21. SAND2018-5463R. link ↩

9. Yang J, Bhattacharya K. Augmented Lagrangian digital image correlation. Experimental Mechanics. 2019 Feb 15;59:187-205. link ↩

10. Augmented Lagrangian Digital Image Correlation (2D_ALDIC) MATLAB code. link ↩

11. Yang J, Rubino V, Ma Z, Tao J, Yin Y, McGhee A, Pan W, Franck C. SpatioTemporally adaptive quadtree mesh (STAQ) digital image correlation for resolving large deformations around complex geometries and discontinuities. Experimental Mechanics. 2022 Sep;62(7):1191-215. link ↩

12. Olufsen SN, Andersen ME, Fagerholt E. μDIC: An open-source toolkit for digital image correlation. SoftwareX. 2020 Jan 1;11:100391. link ↩

13. Olufsen SN. Fork: µDIC: A Python toolkit for Digital Image Correlation (DIC), GitHub: https://github.com/ElsevierSoftwareX/SOFTX_2019_193 link fork of https://github.com/PolymerGuy/muDIC for the 2020 paper. ↩

14. Olufsen SN. docs for µDIC: A toolkit for digital image correlation, docs and repo ↩

15. Rust GPU https://rust-gpu.github.io/ and https://github.com/rust-gpu/rust-gpu ↩

Continuum Mechanics

This section summarizes the kinematics of general (finite) motion, the motion map, the deformation gradient and its Jacobian, the family of finite-strain measures and their linearizations, and the polar and spectral decompositions.

• Upper case indices $(I,J,K)$ denote vector components in the reference configuration and
• Lower case indices $(i,j,k)$ denote vector components in the current configuration.

All vectors are referenced to origin $O$ and dextral, orthonormal basis vectors ${\mathbf{E}}_1$, ${\mathbf{E}}_2$, ${\mathbf{E}}_3$ (equivalently, ${\mathbf{e}}_1$, ${\mathbf{e}}_2$, ${\mathbf{e}}_3$).

Motion

Let the arbitrary time interval be defined as $\mathbb{T}:=[t_0,t_f]$, from initial to final time, inclusive.¹ Let the motion, a one-parameter family of configurations, $\mathbf{\phi }:\mathbb{B}\times \mathbb{T}\mapsto {\mathbb{R}}^3$, map the material particle $\mathbf{X}$ (the reference configuration) into the current configuration $\mathbf{x}$,

$$\mathbf{x}=\mathbf{\phi }(\mathbf{X},t)={\mathbf{\phi }}_t(\mathbf{X}).$$

A motion $\mathbf{\phi }$ evaluated at a particular time $t\in \mathbb{T}$ is referred to as a current configuration or placement. For any placement at time $t$, there is a displacement field $\mathbf{u}:\mathbb{B}\times \mathbb{T}\mapsto {\mathbb{R}}^3$,

$$\mathbf{u}(\mathbf{X},t):=\mathbf{\phi }(\mathbf{X},t)-\mathbf{X}⟺u_i:={\phi }_i-{\delta }_{iJ}{X}_J.$$

Thus, the current configuration $\mathbf{\phi }$ is simply a function of the original placement $\mathbf{X}$, plus a displacement $\mathbf{u}$, which is a function of reference configuration $\mathbf{X}$ and time $t$,

$$\mathbf{\phi }(\mathbf{X},t)=\mathbf{X}+\mathbf{u}(\mathbf{X},t)⟺{\phi }_i={\delta }_{iJ}{X}_J+u_i.$$

The initial condition ${\mathbf{u}}_0$ is found from the initial placement ${\mathbf{\phi }}_0$ and the reference configuration $\mathbf{X}$,

$${\mathbf{u}}_0(\mathbf{X})=\mathbf{u}(\mathbf{X},0)=\mathbf{\phi }(\mathbf{X},0)-\mathbf{X}.$$

Deformation Gradient

To each configuration $\mathbf{\phi }$, we define a deformation gradient $\mathbf{F}:\mathbb{B}\times \mathbb{T}\mapsto {\mathbb{M}}_{+}^3$,

$$\mathbf{F}:={\nabla }_0\mathbf{\phi }(\mathbf{X},t)⟺F_{iJ}:=\frac{\partial {\phi }_i}{\partial X_J}.$$

Real, square matrices of dimension three with positive determinant are denoted ${\mathbb{M}}_{+}^3$. Gradient operations with and without a subscript "$0$" are gradients taken in the reference and current configurations, respectively:

$$\begin{array}{rl}{\nabla }_0(\bullet )& =\frac{\partial (\bullet )}{\partial X_J}\text{(reference configuration gradient)},\\ \nabla (\cdot )& =\frac{\partial (\cdot )}{\partial x_j}\text{(current configuration gradient)}.\end{array}$$

Alternative notations are $\mathrm{GRAD}(\bullet )$ and $\mathrm{grad}(\cdot )$, respectively.

Jacobian of the Deformation Gradient

The Jacobian of the deformation gradient,

$$J:=\det \mathbf{F},$$

describes the (generally non-uniform) volumetric expansion or contraction of the motion from the reference configuration $\mathbf{X}$. All configurations must be admissible in the sense that the Jacobian of the deformation must be positive $(J>0)$. This requirement keeps the deformations from mapping the body to a single, infinitesimally small point $(J=0)$ or turning the body inside-out $(J<0)$.

Isochoric motions preserve the body's total volume. A Jacobian of unity $(J=1)$ describes an isochoric motion. The table below describes the categories of motions (expansion, volume-preserving, contraction, and inadmissible) by Jacobian measure.

Four important isochoric deformations are (1) pure translation, (2) pure rotation, (3) isochoric stretch, and (4) isochoric shear.

Displacement Gradient

From the displacement field defined above, the relationship between the displacement gradient and the deformation gradient is given by

$${\nabla }_0\mathbf{u}={\nabla }_0\mathbf{\phi }-\mathbf{I}=\mathbf{F}-\mathbf{I}⟺\frac{\partial u_i}{\partial X_J}=\frac{\partial {\phi }_i}{\partial X_J}-{\delta }_{iJ}.$$

Right Cauchy-Green Deformation

The right Cauchy-Green deformation arises from the inner product of two differential fiber elements in the reference configuration, $d{\mathbf{X}}_1$ and $d{\mathbf{X}}_2$, mapped by the deformation gradient $\mathbf{F}$ to obtain the inner product of the same differential fibers in the current configuration, $d{\mathbf{x}}_1$ and $d{\mathbf{x}}_2$,

$$d{\mathbf{x}}_1\cdot d{\mathbf{x}}_2=d{\mathbf{x}}_1^{\top }d{\mathbf{x}}_2=(\mathbf{F}d{\mathbf{X}}_1{)}^{\top }\mathbf{F}d{\mathbf{X}}_2=d{\mathbf{X}}_1^{\top }{\mathbf{F}}^{\top }\mathbf{F}d{\mathbf{X}}_2=d{\mathbf{X}}_1^{\top }\mathbf{C}d{\mathbf{X}}_2,$$

where

$$\mathbf{C}:={\mathbf{F}}^{\top }\mathbf{F}⟺C_{IJ}:=F_{Ii}{F}_{iJ}.$$

The right Cauchy-Green deformation tensor $\mathbf{C}$:

1. is defined in the reference configuration,
2. is symmetric and positive-definite,
3. gets its name from the location of the deformation gradient $\mathbf{F}$ in the definition, which is to the right,
4. is a metric that maps fiber lengths from the reference configuration to the current configuration, and
5. is second-order in reference displacement gradients, as shown below:

$$\begin{array}{rl}\mathbf{C}=({\nabla }_0\mathbf{\phi }{)}^{\top }{\nabla }_0\mathbf{\phi }& =(\mathbf{I}+{\nabla }_0\mathbf{u}{)}^{\top }(\mathbf{I}+{\nabla }_0\mathbf{u}),\\ & =\mathbf{I}+\underset{\text{1st order}}{\underbrace{{\nabla }_0\mathbf{u}+({\nabla }_0\mathbf{u}{)}^{\top }}}+\underset{\text{2nd order}}{\underbrace{({\nabla }_0\mathbf{u}{)}^{\top }{\nabla }_0\mathbf{u}}}.\end{array}$$

This result can be expected since, by definition, $\mathbf{C}$ is second-order in the deformation gradient ${\nabla }_0\mathbf{\phi }$, and the relationship between the deformation gradient ${\nabla }_0\mathbf{\phi }$ and the displacement gradient ${\nabla }_0\mathbf{u}$ is linear.

Left Cauchy-Green Deformation

The left Cauchy-Green deformation arises from similar multiplication as with the right Cauchy-Green deformation, but with the stretching going in reverse, from the current configuration back to the reference configuration,

$$d{\mathbf{X}}_1\cdot d{\mathbf{X}}_2=d{\mathbf{X}}_1^{\top }d{\mathbf{X}}_2=({\mathbf{F}}^{-1}d{\mathbf{x}}_1{)}^{\top }{\mathbf{F}}^{-1}d{\mathbf{x}}_2=d{\mathbf{x}}_1^{\top }{\mathbf{F}}^{-\top }{\mathbf{F}}^{-1}d{\mathbf{x}}_2=d{\mathbf{x}}_1^{\top }{\mathbf{b}}^{-1}d{\mathbf{x}}_2,$$

where

$$\mathbf{b}:=\mathbf{F}{\mathbf{F}}^{\top }⟺b_{ij}:=F_{iJ}{F}_{Jj}.$$

The left Cauchy-Green deformation tensor $\mathbf{b}$:

1. is defined in the current configuration,
2. is symmetric and positive-definite,
3. gets its name from the location of the deformation gradient $\mathbf{F}$ in the definition, which is to the left,
4. is a metric whose inverse maps fiber lengths from the current configuration to the reference configuration, and
5. is second-order in current displacement gradients.

Green-Lagrange Strain

The Green-Lagrange strain tensor,

$$\mathbf{E}:=\frac{1}{2}(\mathbf{C}-\mathbf{I})⟺E_{IJ}:=\frac{1}{2}(C_{IJ}-{\delta }_{IJ}),$$

is closely related to the right Cauchy-Green deformation tensor and is often used in defining constitutive law relationships because the measure, when linearized about the reference configuration, coincides with the small strain tensor of linear deformation elasticity, denoted $\mathbf{ϵ}$ and defined in the Infinitesimal Strain section. This relationship can be seen as follows:

$$\begin{array}{rl}2\mathbf{E}& =\mathbf{I}+{\nabla }_0\mathbf{u}+({\nabla }_0\mathbf{u}{)}^{\top }+({\nabla }_0\mathbf{u}{)}^{\top }{\nabla }_0\mathbf{u}-\mathbf{I},\\ 2{\mathbf{E}}_{\text{LIN}}& ={\nabla }_0\mathbf{u}+({\nabla }_0\mathbf{u}{)}^{\top },\end{array}$$

where the higher-order (quadratic) term in the first line is set to zero to achieve the linearized second line.

Almansi-Euler Strain

The Almansi-Euler strain tensor,

$$\mathbf{e}:=\frac{1}{2}(\mathbf{I}-{\mathbf{b}}^{-1})⟺e_{ij}:=\frac{1}{2}\left({\delta }_{ij}-(b_{ij}{)}^{-1}\right),$$

can likewise be used to approximate the small strain tensor $\mathbf{ϵ}$ by combining the definitions of the left Cauchy-Green deformation and the deformation gradient as follows:

$$\begin{array}{rl}2\mathbf{e}& =\mathbf{I}-{\left[{\nabla }_0\mathbf{\phi }({\nabla }_0\mathbf{\phi }{)}^{\top }\right]}^{-1}\\ & ={\delta }_{ij}-{\left[\frac{\partial {\phi }_i}{\partial X_J}\frac{\partial {\phi }_j}{\partial X_J}\right]}^{-1}={\delta }_{ij}-\left[\frac{\partial X_J}{\partial {\phi }_i}\frac{\partial X_J}{\partial {\phi }_j}\right]\\ & =\mathbf{I}-\left[(\mathbf{I}-\nabla \mathbf{u})(\mathbf{I}-\nabla \mathbf{u}{)}^{\top }\right]\\ & =\nabla \mathbf{u}+\nabla {\mathbf{u}}^{\top }-\nabla \mathbf{u}\nabla {\mathbf{u}}^{\top },\\ 2{\mathbf{e}}_{\text{LIN}}& =\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^{\top },\end{array}$$

where the higher-order (quadratic) term is set to zero to achieve the linearized final line.

Small Strain

When displacement gradients are small in the reference configuration,

$${\nabla }_0\mathbf{u}\ll 1,$$

or in the current configuration,

$$\nabla \mathbf{u}\ll 1,$$

respectively, the nonlinear gradient terms are negligible and the finite strain theory simplifies to small strain theory, which occurs when finite strain measures are linearized to obtain ${\mathbf{E}}_{\text{LIN}}$ and ${\mathbf{e}}_{\text{LIN}}$ in the previous sections.

Note that we have restricted the gradients of displacement, and not the displacement $\mathbf{u}$ itself. Thus, displacements between the reference and current configurations can be large (finite), but the gradients of the displacement, either in the reference or current configuration, are small.

The Strain Tensors and Finite Rotations section will demonstrate that the small strain tensors are not suitable to describe motion that contains finite rotation. This makes sense because, in finite rotation, gradients of displacement are large, not small. To adequately describe motion that includes finite rotation, a fully nonlinear strain measure, such as the Seth-Hill strain family, must be used.

Infinitesimal Strain

If we further restrict the small strain theory such that the displacement is small compared to unity,

$$\mathbf{u}\ll 1,$$

the infinitesimal strain theory is obtained, which has no distinction between Lagrangian and Eulerian strain tensors.

In this case, the two small strain tensors, ${\mathbf{E}}_{\text{LIN}}$ and ${\mathbf{e}}_{\text{LIN}}$, converge to a single definition of strain, called the infinitesimal strain tensor $\mathbf{ϵ}$, defined as

$$\mathbf{ϵ}:=\frac{1}{2}\left(\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^{\top }\right).$$

Note that the ${\nabla }_0(\bullet )$ notation has been dropped since the distinction between the reference and current configurations is nonexistent. Also, note that the factor of $\frac{1}{2}$ appears because it then follows that the infinitesimal strain is simply the symmetric part of the displacement gradient,

$$\mathbf{ϵ}=\mathrm{sym}(\nabla \mathbf{u}).$$

Finally, note that the finite Lagrangian and Eulerian strain tensors were defined with the factor of $\frac{1}{2}$ so that their expressions, once linearized and subject to a small displacement assumption, simplify to exactly the infinitesimal strain tensor $\mathbf{ϵ}$.

Seth-Hill Strain Family

We now return to finite strain definitions. Seth and Hill showed that the Green-Lagrange strain tensor $\mathbf{E}$ and the Almansi-Euler strain tensor $\mathbf{e}$ are special cases of the so-called Seth-Hill family of strain measures, defined as

$${\mathbf{E}}^{(m)}:=\{\begin{array}{ll}\frac{1}{m}\left({\mathbf{U}}^m-\mathbf{I}\right)& \text{for }m\ne 0,\\ \ln \mathbf{U}& \text{for }m=0;\text{ and,}\end{array}$$

$${\mathbf{e}}^{(m)}:=\{\begin{array}{ll}\frac{1}{m}\left(\mathbf{I}-{\mathbf{v}}^{-m}\right)& \text{for }m\ne 0,\\ \ln \mathbf{v}& \text{for }m=0.\end{array}$$

The principal stretches ${\lambda }_{\alpha }=l_{\alpha }/L_{\alpha }$, $0<{\lambda }_{\alpha }<\infty$, allow the strain measure ${\mathbf{E}}^{(m)}$ to be written as principal strains, as a function of principal stretch, $f({\lambda }_{\alpha })$,

$$\begin{array}{rl}{\mathbf{E}}^{(m)}& =\sum _{\alpha =1}^{3}f({\lambda }_{\alpha }){\mathbf{N}}_{\alpha }\otimes {\mathbf{N}}_{\alpha }⟺\\ E_{IJ}^{(m)}& =f({\lambda }_1)\left[\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]+f({\lambda }_2)\left[\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right]+f({\lambda }_3)\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 1\end{array}\right],\end{array}$$

where the stretch function

$$f({\lambda }_{\alpha }):=\{\begin{array}{ll}\frac{1}{m}\left({\lambda }_{\alpha }^m-1\right)& \text{for }m\ne 0,\\ \ln {\lambda }_{\alpha }& \text{for }m=0.\end{array}$$

For integer values² of $m\in [-2,2]$, five common strain measures result, listed in the table below, in their three-dimensional and one-dimensional forms. Similar relationships can be constructed for the spatial tensors using

$${\mathbf{e}}^{(m)}=\sum _{\alpha =1}^{3}g({\lambda }_{\alpha }){\mathbf{n}}_{\alpha }\otimes {\mathbf{n}}_{\alpha },g({\lambda }_{\alpha }):=\{\begin{array}{ll}\frac{1}{m}\left(1-{\lambda }_{\alpha }^{-m}\right)& \text{for }m\ne 0,\\ \ln {\lambda }_{\alpha }& \text{for }m=0.\end{array}$$

The one-dimensional strains are illustrated as a function of stretch ratio $\lambda =\ell /L$ in the figure below.

Source: Terpsma et al.³

The figure illustrates several results:

• For small stretches, $\ell \approx L$, (a) the stretch ratio is near unity, $\lambda \approx 1$, (b) the strain values are small, $f(\lambda )\approx 0$, and (c) the tangent of the strains with respect to the stretch ratio is near unity, $df/d\lambda \approx 1$.
• For elongations, $\lambda >1$, the strain monotonically increases since $df/d\lambda >0$ when $\lambda >0$.
• For extreme compressions, $\lambda \approx 0$, (a) the Green-Lagrange strain goes to a value of $-\frac{1}{2}$, (b) the engineering (Biot, nominal) strain tensor goes to a value of $-1$, and (c) the log, true, and Eulerian strains tend to $-\infty$.
• The engineering (Biot, nominal) strain is a linear function of stretch $\lambda$; all other measures are nonlinear functions of stretch $\lambda$.

Neff (2013)⁴ suggested "reasonable requirements" on $f(\lambda ):{\mathbb{R}}^{+}\mapsto \mathbb{R}$, summarized in the table below, wherein a "+" indicates the requirement is satisfied and a "−" indicates the requirement is not satisfied.

The results above illustrate that the log strain ${\mathbf{E}}^{(0)}$ retains more of the desired qualities than any other strain tensor, in the context of finite compression and extension.⁵

• For infinitesimal deformation, all tensors converge to the infinitesimal strain tensor $\mathbf{ϵ}=\mathrm{sym}(\nabla \mathbf{u})$.
• For finite deformation, the Seth-Hill strain measures given by the $f(\lambda )$ function diverge quickly for both large compression and large tension.

The figure below illustrates the one-dimensional strains subtracted from the natural logarithmic strain, $\ln \lambda$, as a function of stretch ratio $\lambda =\ell /L$. The log strain is considered as the finite deformation baseline.

The results show, for example, that in compression at $\lambda =0.8$, the Green-Lagrange strain tensor underreports the log strain by nearly 5%. Such a result illustrates that for finite deformation:

1. Strain measures are not interchangeable, and
2. It is ambiguous and incomplete to simply refer to "strain" as a deformation measure.

For finite deformation, both the strain value and the strain tensor must be specified.

Source: Terpsma et al.³

Strain Tensors and Finite Rotations

Because it takes on nonzero values under finite rotation, the linearized strain tensor should not be used for geometrically nonlinear analysis. These nonzero values are completely artificial and strictly a result of using a linear strain definition with geometrically nonlinear motions. This result is shown as follows.

Let $\mathbf{R}=\mathbf{R}(t)$ be a two-dimensional, rigid body rotation parameterized by time $t$ and scaled by constant $\omega =2\pi$ radians per second. Then, the motion $\mathbf{\phi }$ of a body can be written as

$$\mathbf{\phi }(\mathbf{X},t)=\mathbf{R}(t)\mathbf{X}⟺\left\{\begin{array}{c}{x}_1\\ x_2\end{array}\right\}=\left[\begin{array}{cc}\cos \omega t& \sin \omega t\\ -\sin \omega t& \cos \omega t\end{array}\right]\left\{\begin{array}{c}{X}_1\\ X_2\end{array}\right\}.$$

Then the deformation gradient $\mathbf{F}=\mathbf{F}(t)$ is a function of time alone,

$$\mathbf{F}(t)⟺F_{iJ}=\left[\begin{array}{cc}{F}_{11}& F_{12}\\ F_{21}& F_{22}\end{array}\right]=\left[\begin{array}{cc}\cos \omega t& \sin \omega t\\ -\sin \omega t& \cos \omega t\end{array}\right].$$

The linearized strain tensor is found to be

$$\begin{array}{rl}2{\mathbf{E}}_{\text{LIN}}& ={\nabla }_0\mathbf{u}+({\nabla }_0\mathbf{u}{)}^{\top },\\ & =\mathbf{F}-\mathbf{I}+{\left(\mathbf{F}-\mathbf{I}\right)}^{\top },\\ & =F_{iJ}-{\delta }_{iJ}+F_{Ji}-{\delta }_{Ji},\\ & =\left[\begin{array}{cc}\cos \omega t& \sin \omega t\\ -\sin \omega t& \cos \omega t\end{array}\right]+\left[\begin{array}{cc}\cos \omega t& -\sin \omega t\\ \sin \omega t& \cos \omega t\end{array}\right]-\left[\begin{array}{cc}2& 0\\ 0& 2\end{array}\right],\\ & =2\left(\cos \omega t-1\right)\mathbf{I}.\end{array}$$

Now, for small angles, $\omega t\approx 0$, which is for small deviations $x_1\approx X_1$, $x_2\approx X_2$, then ${\mathbf{E}}_{\text{LIN}}=0$ for rigid body rotations. However, for arbitrary finite angles, $\omega t\ne 0$, and the linearized strain tensor ${\mathbf{E}}_{\text{LIN}}$ reports nonzero strain for rigid body rotations, which is nonsensical.

A correct strain tensor will report zero strain for rigid body rotations. One such strain tensor is the fully nonlinear Green-Lagrange strain tensor. This result is shown as follows:

$$\begin{array}{rl}2\mathbf{E}& ={\mathbf{F}}^{\top }\mathbf{F}-\mathbf{I}\\ & =\left[\begin{array}{cc}\cos \omega t& -\sin \omega t\\ \sin \omega t& \cos \omega t\end{array}\right]\left[\begin{array}{cc}\cos \omega t& \sin \omega t\\ -\sin \omega t& \cos \omega t\end{array}\right]-\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]=0.\end{array}$$

Polar Decomposition

Given the rotation tensor $\mathbf{R}$, the material stretch tensor $\mathbf{U}$, and the spatial stretch tensor $\mathbf{v}$, the deformation gradient $\mathbf{F}$ has the multiplicative decomposition,

$$\mathbf{F}=\mathbf{R}\mathbf{U}=\mathbf{v}\mathbf{R}⟺F_{iJ}=R_{iK}{U}_{KJ}=v_{ik}{R}_{kJ}.$$

Here we have a slight abuse of notation, where intermediate configurations that have stretched but not yet rotated are denoted with capital letter indices. Thus the "$K$" subscript in $R_{iK}{U}_{KJ}$ is an intermediate stretched but non-rotated configuration.

The stretch tensors $\mathbf{U}$ and $\mathbf{v}$ are both symmetric and positive definite. The rotation tensor $\mathbf{R}$ is non-symmetric and orthogonal. The figure below shows the polar decomposition about a material point $\mathbf{X}$ and fibers $d\mathbf{X}$ in its vicinity mapped to the spatial point $\mathbf{x}$ with the same fibers mapped to $d\mathbf{x}$.

Principal Stretches and Axes

The stretch tensors $\mathbf{U}$ and $\mathbf{v}$ have the same eigenvalues, $({\lambda }_1,{\lambda }_2,{\lambda }_3)$, called principal stretches. For non-trivial rotations, i.e. $\mathbf{R}\ne \mathbf{I}$, $\mathbf{U}$ and $\mathbf{v}$ have unique eigenvectors, called principal stretch directions. The principal stretch directions of $\mathbf{U}$ are $\left\{{\mathbf{N}}_1,{\mathbf{N}}_2,{\mathbf{N}}_3\right\}$. The principal stretch directions of $\mathbf{v}$ are $\left\{{\mathbf{n}}_1,{\mathbf{n}}_2,{\mathbf{n}}_3\right\}$. The two sets of eigenvectors are related through rotation $\mathbf{R}$,

$${\mathbf{n}}_1=\mathbf{R}{\mathbf{N}}_1,{\mathbf{n}}_2=\mathbf{R}{\mathbf{N}}_2,{\mathbf{n}}_3=\mathbf{R}{\mathbf{N}}_3;$$

or generally,

$${\mathbf{n}}_{\alpha }=\mathbf{R}{\mathbf{N}}_{\alpha }\text{for }\alpha =1,2,3.$$

Spectral Representation

The deformation gradient, its polar decomposition, and the Cauchy-Green deformations have spectral decompositions in terms of the principal stretches and stretch directions,

$$\begin{array}{rlrl}\mathbf{F}& =\sum _{\alpha =1}^3{\lambda }_{\alpha }{\mathbf{n}}_{\alpha }\otimes {\mathbf{N}}_{\alpha },& \mathbf{U}& =\sum _{\alpha =1}^3{\lambda }_{\alpha }{\mathbf{N}}_{\alpha }\otimes {\mathbf{N}}_{\alpha },\\ \mathbf{R}& =\sum _{\alpha =1}^3{\mathbf{n}}_{\alpha }\otimes {\mathbf{N}}_{\alpha },& \mathbf{v}& =\sum _{\alpha =1}^3{\lambda }_{\alpha }{\mathbf{n}}_{\alpha }\otimes {\mathbf{n}}_{\alpha },\\ \mathbf{C}& =\sum _{\alpha =1}^3{\lambda }_{\alpha }^2{\mathbf{N}}_{\alpha }\otimes {\mathbf{N}}_{\alpha },& \mathbf{b}& =\sum _{\alpha =1}^3{\lambda }_{\alpha }^2{\mathbf{n}}_{\alpha }\otimes {\mathbf{n}}_{\alpha }.\end{array}$$

The Green-Lagrange strain tensor $\mathbf{E}$ and the Almansi-Euler strain tensor $\mathbf{e}$, in principal stretches and stretch directions, are

$$\mathbf{E}=\sum _{\alpha =1}^3\frac{1}{2}\left({\lambda }_{\alpha }^2-1\right){\mathbf{N}}_{\alpha }\otimes {\mathbf{N}}_{\alpha },\mathbf{e}=\sum _{\alpha =1}^3\frac{1}{2}\left(1-{\lambda }_{\alpha }^{-2}\right){\mathbf{n}}_{\alpha }\otimes {\mathbf{n}}_{\alpha }.$$

The generalization of the Seth-Hill material strain tensor ${\mathbf{E}}^{(m)}$ and spatial strain tensor ${\mathbf{e}}^{(m)}$, in principal stretches and stretch directions, are

$${\mathbf{E}}^{(m)}=\sum _{\alpha =1}^3\frac{1}{m}\left({\lambda }_{\alpha }^m-1\right){\mathbf{N}}_{\alpha }\otimes {\mathbf{N}}_{\alpha },{\mathbf{e}}^{(m)}=\sum _{\alpha =1}^3\frac{1}{m}\left(1-{\lambda }_{\alpha }^{-m}\right){\mathbf{n}}_{\alpha }\otimes {\mathbf{n}}_{\alpha },$$

and the relationship between the two strain tensors is given through a rotation $\mathbf{R}$ transformation,

$${\mathbf{e}}^{(-m)}=\mathbf{R}{\mathbf{E}}^{(m)}{\mathbf{R}}^{\top }.$$

In the case when $m\to \infty$, the material and spatial logarithmic strain tensors, also known as the Hencky material and spatial strain tensors, $\mathbf{H}$ and $\mathbf{h}$, are obtained as⁶

$$\begin{array}{rl}\mathbf{H}=\ln \mathbf{U}={\mathbf{E}}^{(0)}& =\sum _{\alpha =1}^3\ln {\lambda }_{\alpha }{\mathbf{N}}_{\alpha }\otimes {\mathbf{N}}_{\alpha },\text{ and}\\ \mathbf{h}=\ln \mathbf{v}={\mathbf{e}}^{(0)}& =\sum _{\alpha =1}^3\ln {\lambda }_{\alpha }{\mathbf{n}}_{\alpha }\otimes {\mathbf{n}}_{\alpha }.\end{array}$$

Two concrete illustrations follow: Rigid Body Motion works through pure translation as the simplest possible deformation, and Simple Shear works through an isochoric shear in closed form, computing $\mathbf{F}$, $\mathbf{C}$, $\mathbf{E}$, $\mathbf{b}$, and $\mathbf{e}$ explicitly.

1. Note that $t_0$, while typically zero, may be any real number less than $t_f$. ↩

2. Technically, $m$ can be any real number, not just an integer. ↩

3. Terpsma RJ, Hovey CB. Blunt impact brain injury using cellular injury criterion. Sandia National Lab. (SNL-NM), Albuquerque, NM (United States); 2020 Oct 1. link ↩ ↩2

4. Neff, P. (2013). The Hencky strain measure is the geodesic distance to SO($n$), at 6. ↩

5. The Bažant strain, $f(\lambda )=\frac{1}{2}\left(\lambda -\frac{1}{\lambda }\right)$, not considered here, also satisfies $-f(\lambda )=f\left(\frac{1}{\lambda }\right)$. ↩

6. See Xiao H, Bruhns OT, Meyers A. Hypo-elasticity model based upon the logarithmic stress rate. Journal of Elasticity. 1997 Apr 1;47(1):51-68, at page 54, Eq. (2.2). ↩

Rigid Body Motion

Translation

A translation is a homogeneous deformation of the form

$$\mathbf{\phi }(\mathbf{X})=\bar{\mathbf{u}}+\mathbf{X}.$$

This deformation occurs when the displacement $\mathbf{u}(\mathbf{X},t)$ is a constant $\bar{\mathbf{u}}$, and thus not a function of reference position $\mathbf{X}$ or time $t$.

For the concrete example in the figure below, let $\bar{\mathbf{u}}=⟨{\bar{u}}_1,{\bar{u}}_2{⟩}^{\top }=⟨3,2{⟩}^{\top }$. In this case, we see the placement $\mathbf{\phi }$ moves right and up on the page, relative to the reference configuration $\mathbf{X}$, by an amount of 3 and 2, respectively. The reference configuration is shown in blue. The current configuration is shown in red.

Simple Shear

Note: The source code for this section is listed below.

In this section, we cover simple shear, a deformation that illustrates several concepts without being too complicated. That said, we will find that simple shear isn't quite as trivial as its name would suggest. In fact, a paper titled Simple shear is not so simple¹ reviews some details not discussed here.

The figure below illustrates simple shear, with relative motion of planes normal to the vertical axis. For all configurations (reference and current):

• Horizontal fibers remain horizontal.
• The length of horizontal fibers remains constant.
• The vertical distance between the horizontal fibers remains constant.
• The body's volume is constant. The deformation is isochoric since $\det (\mathbf{F})=1$.

The relative motion is characterized by the non-dimensional ratio of length $a$ per unit height, where $a\in \mathbb{R}\subset [0,\infty )$. The angle, $\gamma \in \mathbb{R}\subset [0,\pi /2)$, subtended by simple shear is $\gamma =\arctan (a)$. In the limit as $a\mapsto \infty$, the shear angle $\gamma \mapsto \pi /2$. For small values of $a$, the small-angle assumption is valid, with $\tan (\gamma )\approx \gamma \approx a$.

Source: simple_shear.py

Rate-Independent Form

The simple shear deformation $\mathbf{\phi }$ in rate-independent form is

$$\mathbf{\phi }(\mathbf{X})=(X_1+a{X}_2){\mathbf{e}}_1+X_2{\mathbf{e}}_2+X_3{\mathbf{e}}_3.$$

The deformation gradient $\mathbf{F}$ is

$$\begin{array}{rl}\mathbf{F}& :={\mathbf{e}}_i\frac{\partial {\phi }_i}{\partial X_J}{\mathbf{E}}_J={\mathbf{e}}_1\otimes {\mathbf{E}}_1+a{\mathbf{e}}_1\otimes {\mathbf{E}}_2+{\mathbf{e}}_2\otimes {\mathbf{E}}_2+{\mathbf{e}}_3\otimes {\mathbf{E}}_3\\ & =\{{\mathbf{e}}_1{\mathbf{e}}_2{\mathbf{e}}_3\}\left[\begin{array}{ccc}1& a& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\left\{\begin{array}{c}{\mathbf{E}}_1\\ {\mathbf{E}}_2\\ {\mathbf{E}}_3\end{array}\right\}⟺{\mathbf{F}}_{iJ}=\left[\begin{array}{ccc}1& a& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right].\end{array}$$

Note that the volume remains constant for all deformations since $\det (\mathbf{F})=1$. The deformation thus belongs to the group of isochoric motions. The right Cauchy-Green strain, $\mathbf{C}:={\mathbf{F}}^T\mathbf{F}$, in simple shear, is

$$\begin{array}{rl}\mathbf{C}=& {\mathbf{E}}_1\otimes {\mathbf{E}}_1+a({\mathbf{E}}_1\otimes {\mathbf{E}}_2+{\mathbf{E}}_2\otimes {\mathbf{E}}_1)+\\ & (a^2+1){\mathbf{E}}_2\otimes {\mathbf{E}}_2+{\mathbf{E}}_3\otimes {\mathbf{E}}_3⟺{\mathbf{C}}_{IJ}=\left[\begin{array}{ccc}1& a& 0\\ a& a^2+1& 0\\ 0& 0& 1\end{array}\right].\end{array}$$

The principal directions of $\mathbf{C}$ are given with the three eigenvectors $\{{\mathbf{N}}_1,{\mathbf{N}}_2,{\mathbf{N}}_3\}$ and their respective eigenvalues ${\lambda }_1^2,{\lambda }_2^2,{\lambda }_3^2$,

$$\mathbf{C}=\sum _{\alpha =1}^3{\lambda }_{\alpha }^2{\mathbf{N}}_{\alpha }\otimes {\mathbf{N}}_{\alpha }.$$

For simple shear, the eigenvalues of $\mathbf{C}$ and their corresponding eigenvectors, satisfying $[\mathbf{C}]\{\mathbf{N}\}={\lambda }^2\{\mathbf{N}\}$, are

$$\begin{array}{rl}{\lambda }_1^2=\frac{1}{2}\left(a^2+a\sqrt{a^2+4}+2\right),& \text{and}{\mathbf{N}}_1={\left\{\frac{1}{2}\left(-a+\sqrt{a^2+4}\right),1,0\right\}}^T;\\ {\lambda }_2^2=1,& \text{and}{\mathbf{N}}_2=\{0,0,1{\}}^T;\\ {\lambda }_3^2=\frac{1}{2}\left(a^2-a\sqrt{a^2+4}+2\right),& \text{and}{\mathbf{N}}_3={\left\{\frac{1}{2}\left(-a-\sqrt{a^2+4}\right),1,0\right\}}^T.\end{array}$$

The Green-Lagrange strain, $\mathbf{E}:=\frac{1}{2}(\mathbf{C}-\mathbf{I})$, in simple shear is

$${\mathbf{E}}_{IJ}=\frac{1}{2}\left[\begin{array}{ccc}0& a& 0\\ a& a^2& 0\\ 0& 0& 0\end{array}\right].$$

The left Cauchy-Green strain, $\mathbf{b}:=\mathbf{F}{\mathbf{F}}^T$, (and its inverse), in simple shear, are

$${\mathbf{b}}_{ij}=\left[\begin{array}{ccc}{a}^2+1& a& 0\\ a& 1& 0\\ 0& 0& 1\end{array}\right],\text{and}{\mathbf{b}}_{ji}^{-1}=\left[\begin{array}{ccc}1& -a& 0\\ -a& a^2+1& 0\\ 0& 0& 1\end{array}\right].$$

The Almansi-Euler strain, $\mathbf{e}:=\frac{1}{2}(\mathbf{I}-{\mathbf{b}}^{-1})$, for simple shear is

$${\mathbf{e}}_{ij}=\frac{1}{2}\left[\begin{array}{ccc}0& a& 0\\ a& -a^2& 0\\ 0& 0& 0\end{array}\right].$$

Source

simple_shear.py

"""Demonstrate simple shear deformation of a square body."""

from typing import Tuple
import os
import numpy as np
from numpy.typing import NDArray

from matplotlib import rc
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator

rc("text", usetex=True)
# rc('font', **{'family': 'serif', 'serif': ['Computer Modern Roman']})
rc("font", family="serif")


def rotate(x0: NDArray, y0: NDArray, rotation: float) -> Tuple[NDArray, NDArray]:
    """Given list of reference points (x0, y0), rotate them about the
    z-axis by rotation angle (radians) to the current points (x1, y1).
    """
    x1 = np.cos(rotation) * x0 - np.sin(rotation) * y0
    y1 = np.sin(rotation) * x0 + np.cos(rotation) * y0
    return x1, y1


def simple_shear(x0: NDArray, y0: NDArray, shear: float) -> Tuple[NDArray, NDArray]:
    """Given a list of reference points (x0, y0), simple shear them in
    the x-axis by distance shear_x (Length) to the current points (x1, y1).
    """
    x1 = x0 + shear * y0
    y1 = y0
    return x1, y1


def draw_shape(
    axis,
    xs: np.ndarray,
    ys: np.ndarray,
    color: str = "dimgray",
    linestyle: str = "-",
    alpha: float = 1.0,
):
    """Draw the body and points."""

    # Draw body outline
    axis.plot(
        xs,
        ys,
        color=color,
        linestyle=linestyle,
        alpha=alpha,
    )  # body outline

    # Plot points on the body
    axis.plot(xs, ys, "o", color=color, markersize=5, alpha=alpha)


SHOW = False
SAVE = True

fig = plt.figure(figsize=(6, 3))  # inches, (wide, tall)
ax1 = fig.add_subplot(1, 2, 1)
ax2 = fig.add_subplot(1, 2, 2)

# Subfigure 1

OFFSET_X = 0
OFFSET_Y = 0
RADTODEG = 180.0 / np.pi
DEGTORAD = 1.0 / RADTODEG

# Defined reference configuration path, in a closed circle
xs_0 = np.array([0, 1, 1, 0, 0])
ys_0 = np.array([0, 0, 1, 1, 0])

# Draw reference shape
draw_shape(axis=ax1, xs=xs_0, ys=ys_0, color="dimgray", linestyle="-", alpha=0.5)

# Draw first sheared shape
SHEAR_1 = 0.5
(xs_1, ys_1) = simple_shear(xs_0, ys_0, shear=SHEAR_1)
draw_shape(axis=ax1, xs=xs_1, ys=ys_1, color="green", linestyle=":", alpha=0.9)

# Draw second sheared shape
SHEAR_2 = 1.0
(xs_2, ys_2) = simple_shear(xs_0, ys_0, shear=SHEAR_2)
draw_shape(axis=ax1, xs=xs_2, ys=ys_2, color="red", linestyle="--", alpha=0.9)

# Draw origin point
ax1.plot(0, 0, "o", color="black", label="origin = (0, 0, 0)")
# Draw origin label
ax1.text(0.25, -0.25, r"$O, o$", ha="center", va="center")

# Draw 1:a rise over run lines
SCALE = 0.90  # scale
HAIRLINE_OFFSET_Y = 0.1
epsx, epsy = 0.125, 0.25 + HAIRLINE_OFFSET_Y
slope_x, slope_y = np.array([0, 0, 0.25]) * SCALE, np.array([0, 0.5, 0.5])
ax1.plot(slope_x + epsx, slope_y + epsy, lw=0.5, color="green")
ax1.text(0.125, 0.5, "1", color="green", ha="right", va="center")
ax1.text(0.25, 0.88, r"$a$", color="green", ha="center")


# SHEAR_12 = 0.5  # Length units, shear in the X_1 direction
# draw(ax1, ux=OFFSET_Y, uy=OFFSET_Y, ur=dr, shear=SHEAR_12, t0=False, c="green", ls=":")

# SHEAR_12 = 1.0  # Length units, shear in the X_1 direction, larger shear
# draw(ax1, ux=OFFSET_X, uy=OFFSET_X, ur=dr, shear=SHEAR_12, t0=False, c="red", ls="--")


ax1.axis("equal")
# ax2.axis('equal')
# major axes
ax1.xaxis.set_major_locator(MultipleLocator(1.0))
ax1.yaxis.set_major_locator(MultipleLocator(1.0))
ax2.xaxis.set_major_locator(MultipleLocator(1.0))
ax2.yaxis.set_major_locator(MultipleLocator(1.0))
# minor axes
# ax1.xaxis.set_minor_locator(MultipleLocator(0.5))
# ax1.yaxis.set_minor_locator(MultipleLocator(0.5))

ax1.grid(
    visible=True, which="major", linestyle="solid", linewidth=0.5, color="lightgray"
)  # FIX APPLIED HERE
ax2.grid(
    visible=True, which="major", linestyle="solid", linewidth=0.5, color="lightgray"
)  # FIX APPLIED HERE
#
# ax.grid(b=True, which='minor', linestyle=':')
ax1.set_xlabel(r"configuration $X_1, x_1$")
ax1.set_ylabel(r"configuration $X_2, x_2$")

# Subfigure 2

X_MIN = 0
X_MAX = 10
epsx, epsy = 0.4, np.pi / 16
x = np.linspace(X_MIN, X_MAX)
y = np.arctan(x)
ax2.plot(x, y, linewidth=2, color="blue")
ax2.text(
    X_MAX - epsx,
    np.pi / 2 + epsy / 2,
    r"$\gamma \mapsto \frac{\pi}{2}$",
    ha="right",
    backgroundcolor="white",
)
ax2.plot(
    [X_MIN, X_MAX],
    np.pi / 2 * np.array([1, 1]),
    lw=2,
    alpha=0.5,
    color="black",
    linestyle="--",
    zorder=4,
)

ax2.plot(0, 0, "o", color="dimgray", alpha=0.5, zorder=4)
ax2.text(0 + epsx, 0 - epsy, r"$(0, 0)$", backgroundcolor="white")

ax2.plot(0.5, 0.46, "o", color="green", alpha=0.9, zorder=4)
ax2.text(0.50 + epsx, 0.46 - epsy, r"$(0.50, 0.46)$", backgroundcolor="white")

ax2.plot(1, np.pi / 4, "o", color="red", alpha=0.9, zorder=4)
ax2.text(1 + epsx, np.pi / 4 - epsy, r"$(1, \frac{\pi}{4})$", backgroundcolor="white")

ax2.set_xlabel(r"non-dimensional distance $a\;[l/L]$")
ax2.set_ylabel(r"$\gamma = \arctan(a)$ [rad]")
# https://matplotlib.org/3.1.1/gallery/ticks_and_spines/tick-locators.html
# ax2.xaxis.set_major_locator(ticker.FixedLocator([0, 5, 10]))
ax2.set_xticks([0, 5, 10])
ax2.set_yticks([0, np.pi / 4, np.pi / 2])
ax2.set_yticklabels(["0", r"$\frac{\pi}{4}$", r"$\frac{\pi}{2}$"])

ax1.set_xlim(-epsx, 2 + epsx)
ax2.set_xlim(0 - epsx, X_MAX + epsx)
eps = np.pi / 8
ax2.set_ylim(0 * np.pi / 4 - eps, np.pi / 2 + eps)

if SHOW:
    plt.show()

if SAVE:
    script_name = os.path.basename(__file__)
    figure_name = os.path.splitext(script_name)[0]
    print(f"Saving figure as {figure_name}.pdf")
    fig.savefig(figure_name + ".pdf", bbox_inches="tight")

References

1. Destrade M, Murphy JG, Saccomandi G. Simple shear is not so simple. International Journal of Non-Linear Mechanics. 2012 Mar 1;47(2):210-4. download ↩

Finite Element Method

A finite element mesh is a collection of nodes (points) connected into elements — small regions used to interpolate a quantity of interest (e.g. displacement) across the whole domain. Multi-Point Motion (covered later in this book) builds exactly this kind of nodal point collection and, in Tracking the Grid, tracks every one of its 12 points to its current position — exactly the per-node displacement data a mesh needs to represent how an object deformed.

Kernel Size First also covers the kernel-size-versus-point-spacing tradeoff involved in getting that per-node data reliably — the same considerations apply whether the points come from a toy grid or a real mesh.

Once every node's current position is known, an actual finite element mesh still needs one more thing this page doesn't provide: element connectivity — which nodes join together into which elements. Building that connectivity is future work, not implemented here; what this page covers below is the element formulation that connectivity would plug into — shape functions, strain, and deformation gradient, for the four-node quadrilateral element.

Four-Noded Quadrilateral Finite Element (Q4)

The four-node quadrilateral element is one of the most commonly used elements in 2D FEA. It has four corner nodes, with each node having two degrees of freedom (DOFs): displacements in the $X$ and $Y$ directions.

Image credit: James et al.¹

Note: Since we are using a finite deformation continuum mechanics framework, we will use $X$ and $Y$ (instead of $x$ and $y$ in the James et al. figure above).

Shape Functions

For the element in local coordinates $(\xi ,\eta )\in [-1,1]\times [-1,1]$, the bilinear shape functions are defined:

$$\begin{array}{rl}{N}_1(\xi ,\eta )& :=\frac{1}{4}(1-\xi )(1-\eta )\\ N_2(\xi ,\eta )& :=\frac{1}{4}(1+\xi )(1-\eta )\\ N_3(\xi ,\eta )& :=\frac{1}{4}(1+\xi )(1+\eta )\\ N_4(\xi ,\eta )& :=\frac{1}{4}(1-\xi )(1+\eta )\end{array}$$

The shape functions satisfy the following properties:

• Kronecker delta property: $N_a({\xi }_b,{\eta }_b)={\delta }_{ab}$ (equals 1 at node $a$, 0 at other nodes)
• Partition of unity: ${\sum }_{a=1}^4{N}_a(\xi ,\eta )=1$ for all $(\xi ,\eta )$

Local Coordinates

The key concept in finite element analysis is the isoparametric mapping between the local coordinate system and the global coordinate system.

This mapping allows:

• Integration to be performed on the local domain (parent element)
• Handling of arbitrarily shaped quadrilaterals
• Use of the same shape functions for geometry and displacement (isoparametric concept)

The isoparametric coordinates $(\xi ,\eta )$ range from $-1$ to $+1$ in both the $X$ and $Y$ directions.

The mapping between global coordinates $(X,Y)$ and local coordinates is introduced as a linear combination of local shape functions $N_{\mathrm{node}}(\mathbf{X})$:

$$\mathbf{X}(\mathbf{\xi })=\left\{\begin{array}{c}X(\xi ,\eta )\\ Y(\xi ,\eta )\end{array}\right\}=\sum _{a=1}^4{N}_a(\mathbf{\xi }){\mathbf{X}}_a=\sum _{a=1}^4{N}_a(\xi ,\eta )\left\{\begin{array}{c}{X}_a\\ Y_a\end{array}\right\}$$

where $(X_a,Y_a)$ is the position of node $a$, and $a=1\dots 4$.

Shape Function Derivatives in Local Coordinates

The derivatives with respect to the local coordinate system are

$$\frac{\partial N_a(\xi ,\eta )}{\partial \xi }=\frac{1}{4}\left[\begin{array}{c}-(1-\eta )\\ (1-\eta )\\ (1+\eta )\\ -(1+\eta )\end{array}\right],\frac{\partial N_a(\xi ,\eta )}{\partial \eta }=\frac{1}{4}\left[\begin{array}{c}-(1-\xi )\\ -(1+\xi )\\ (1+\xi )\\ (1-\xi )\end{array}\right]$$

These are assembled into a $2\times 4$ matrix:

$$\frac{\partial \mathbf{N}(\mathbf{\xi })}{\partial \mathbf{\xi }}={\left[\begin{array}{cccc}\frac{\partial N_1}{\partial \xi }& \frac{\partial N_2}{\partial \xi }& \frac{\partial N_3}{\partial \xi }& \frac{\partial N_4}{\partial \xi }\\ \frac{\partial N_1}{\partial \eta }& \frac{\partial N_2}{\partial \eta }& \frac{\partial N_3}{\partial \eta }& \frac{\partial N_4}{\partial \eta }\end{array}\right]}_{2\times 4}$$

Jacobian Matrix

The Jacobian matrix ${\mathbf{j}}_0$ relates derivatives in local coordinates to derivatives in global coordinates. It is important to include the "matrix" part of "Jacobian matrix". It is distinct from the Jacobian of the Deformation Gradient $J$, which is a scalar value (not a matrix). For nodal coordinates organized as:

$$\mathbf{X}={\left[\begin{array}{cc}{X}_1& Y_1\\ X_2& Y_2\\ X_3& Y_3\\ X_4& Y_4\end{array}\right]}_{4\times 2}$$

the Jacobian matrix ${\mathbf{j}}_0$ is computed as:

$${\mathbf{j}}_0(\mathbf{\xi }):=\frac{\partial \mathbf{X}(\mathbf{\xi })}{\partial \mathbf{\xi }}={\left[\begin{array}{cc}\frac{\partial X}{\partial \xi }& \frac{\partial Y}{\partial \xi }\\ \frac{\partial X}{\partial \eta }& \frac{\partial Y}{\partial \eta }\end{array}\right]}_{2\times 2}=\frac{\partial \mathbf{N}(\mathbf{\xi })}{\partial \mathbf{\xi }}\cdot \mathbf{X}$$

The individual components (dropping the ${\bullet }_0$ reference configuration notation to avoid subscript confusion) are:

$$\begin{array}{rl}{j}_{11}& =\frac{\partial X}{\partial \xi }=\sum _{a=1}^4\frac{\partial N_a}{\partial \xi }{X}_a\\ j_{12}& =\frac{\partial Y}{\partial \xi }=\sum _{a=1}^4\frac{\partial N_a}{\partial \xi }{Y}_a\\ j_{21}& =\frac{\partial X}{\partial \eta }=\sum _{a=1}^4\frac{\partial N_a}{\partial \eta }{X}_a\\ j_{22}& =\frac{\partial Y}{\partial \eta }=\sum _{a=1}^4\frac{\partial N_a}{\partial \eta }{Y}_a\end{array}$$

Shape Function Derivatives in Global Coordinates

The transformation from local to global coordinate derivatives requires the inverse Jacobian matrix through the chain rule. Since

$$\frac{\partial N_a}{\partial X}=\frac{\partial N_a}{\partial \xi }\frac{\partial \xi }{\partial X}+\frac{\partial N_a}{\partial \eta }\frac{\partial \eta }{\partial X}$$

$$\frac{\partial N_a}{\partial Y}=\frac{\partial N_a}{\partial \xi }\frac{\partial \xi }{\partial Y}+\frac{\partial N_a}{\partial \eta }\frac{\partial \eta }{\partial Y}$$

then

$$\left[\begin{array}{c}\frac{\partial N_a}{\partial X}\\ \frac{\partial N_a}{\partial Y}\end{array}\right]={\left[\begin{array}{cc}\frac{\partial \xi }{\partial X}& \frac{\partial \eta }{\partial X}\\ \frac{\partial \xi }{\partial Y}& \frac{\partial \eta }{\partial Y}\end{array}\right]}_{2\times 2}\left[\begin{array}{c}\frac{\partial N_a}{\partial \xi }\\ \frac{\partial N_a}{\partial \eta }\end{array}\right]={\mathbf{j}}_0^{-1}(\mathbf{X})\left[\begin{array}{c}\frac{\partial N_a}{\partial \xi }\\ \frac{\partial N_a}{\partial \eta }\end{array}\right]$$

In matrix form for all shape functions:

$$\frac{\partial \mathbf{N}}{\partial \mathbf{X}}={\mathbf{j}}_0^{-1}\cdot \frac{\partial \mathbf{N}}{\partial \mathbf{\xi }}$$

where:

$$\frac{\partial \mathbf{N}}{\partial \mathbf{X}}={\left[\begin{array}{cccc}\frac{\partial N_1}{\partial X}& \frac{\partial N_2}{\partial X}& \frac{\partial N_3}{\partial X}& \frac{\partial N_4}{\partial X}\\ \frac{\partial N_1}{\partial Y}& \frac{\partial N_2}{\partial Y}& \frac{\partial N_3}{\partial Y}& \frac{\partial N_4}{\partial Y}\end{array}\right]}_{2\times 4}$$

Displacement Field

The displacement $\mathbf{u}(\mathbf{X})$ is defined as the difference between the current configuration $\mathbf{\phi }(\mathbf{X})$ and the reference configuration $\mathbf{X}$,

$$\mathbf{u}(\mathbf{X}):=\mathbf{\phi }(\mathbf{X})-\mathbf{X}$$

The displacement field within the element is interpolated using shape functions:

$$\mathbf{u}(\xi ,\eta )=\left\{\begin{array}{c}u(\xi ,\eta )\\ v(\xi ,\eta )\end{array}\right\}=\sum _{a=1}^4{N}_a(\mathbf{\xi }){\mathbf{u}}_a=\sum _{a=1}^4{N}_a(\xi ,\eta )\left\{\begin{array}{c}{u}_a\\ v_a\end{array}\right\}$$

where $(u_a,v_a)$ is the respective $(X,Y)$ displacement of node $a$, and $a=1\dots 4$.

Displacement Gradient

$${\nabla }_0\mathbf{u}(\mathbf{X})={\left[\begin{array}{cc}\frac{\partial u}{\partial X}& \frac{\partial u}{\partial Y}\\ \frac{\partial v}{\partial X}& \frac{\partial v}{\partial Y}\end{array}\right]}_{2\times 2}$$

Each component is computed using the chain rule:

$$\begin{array}{rl}\frac{\partial u}{\partial X}& =\sum _{a=1}^4\frac{\partial N_a}{\partial X}{u}_a\\ \frac{\partial u}{\partial Y}& =\sum _{a=1}^4\frac{\partial N_a}{\partial Y}{u}_a\\ \frac{\partial v}{\partial X}& =\sum _{a=1}^4\frac{\partial N_a}{\partial X}{v}_a\\ \frac{\partial v}{\partial Y}& =\sum _{a=1}^4\frac{\partial N_a}{\partial Y}{v}_a\end{array}$$

In compact matrix notation:

$${\nabla }_0\mathbf{u}={\left(\frac{\partial \mathbf{N}}{\partial \mathbf{X}}\cdot \mathbf{u}\right)}^T$$

where $\mathbf{u}$ is the nodal displacement matrix:

$$\mathbf{u}={\left[\begin{array}{cc}{u}_1& v_1\\ u_2& v_2\\ u_3& v_3\\ u_4& v_4\end{array}\right]}_{4\times 2}$$

See Displacement Gradient for more information.

Deformation Gradient

The deformation gradient tensor $\mathbf{F}(\mathbf{X})$ maps material points in the reference configuration $\mathbf{X}$ to their positions in the current (deformed) configuration $\mathbf{\phi }$:

$$\mathbf{F}(\mathbf{X}):=\frac{\partial \mathbf{\phi }(\mathbf{X})}{\partial \mathbf{X}}$$

Because $\mathbf{u}(\mathbf{X})=\mathbf{\phi }(\mathbf{X})-\mathbf{X}$,

$${\nabla }_0\mathbf{u}=\mathbf{F}-\mathbf{I}⟹\mathbf{F}={\nabla }_0\mathbf{u}+\mathbf{I}$$

Explicitly:

$$\mathbf{F}=\left[\begin{array}{cc}{F}_{11}& F_{12}\\ F_{21}& F_{22}\end{array}\right]=\left[\begin{array}{cc}\frac{\partial u}{\partial X}+1& \frac{\partial u}{\partial Y}\\ \frac{\partial v}{\partial X}& \frac{\partial v}{\partial Y}+1\end{array}\right]$$

The determinant $J:=\det (\mathbf{F})$ represents the local volume ratio and must be positive for physically admissible deformations. See Deformation Gradient and Jacobian of the Deformation Gradient for more information.

Gauss Points

To evaluate quantities that depend on the displacement field and its gradient, such as strain, we use Gaussian Quadrature. We don't typically calculate quantities at the nodes. Rather, we quantify strain at specific integration points (also known as Gauss points) where mathematical precision is the highest.

For a 2D quadrilateral element, we typically use a $2\times 2$ Gauss rule. The integration points are located in the local coordinate system $(\xi ,\eta )$ at

$$\xi ,\eta \in \left\{\pm \frac{1}{\sqrt{3}}\right\}\approx \pm 0.57735$$

References

1. James KA, Lee E, Martins JR. Stress-based topology optimization using an isoparametric level set method. Finite Elements in Analysis and Design. 2012 Oct 1;58:20-30. link ↩

Image Generation

The source for the commands on this page is dictk's own rosta, checkerboard, and astronaut subcommands — see dictk --help.

Note: the images embedded on this page are rendered as PNG (--format png), not dictk's default TIFF. Browsers don't natively render TIFF in <img> tags, so a TIFF embedded here simply wouldn't display.

Among the alternatives, PNG also wins on its own merits: it is lossless, whereas JPG's compression tends to smear hard edges and speckle-pattern detail (for the 200x200 checkerboard on this page: TIFF 40,256 bytes, JPG 6,760 bytes, PNG only 418 bytes — JPG is actually larger than PNG here, because its block-based compression is a poor fit for hard-edged content like a checkerboard). SVG doesn't help either: since there's no vector structure to trace, dictk's SVG output just wraps that same PNG in a base64-encoded XML container, which comes out to 809 bytes here — roughly double the raw PNG for no rendering benefit.

TIFF remains dictk's command-line default, since it's the lossless, uncompressed format conventionally used for DIC and other scientific-imaging workflows.

CLI vs. API: the Command Line Interface (CLI) subcommands on this page (dictk rosta, dictk checkerboard, dictk astronaut) write an image file to disk — that's their whole job. The corresponding Python functions, dictk.rosta, dictk.checkerboard, and dictk.astronaut, take the same parameters but perform no file I/O: they return a NumPy array only. That keeps the Python API composable in a functional style — arrays can be piped through further functions (e.g. combine below) before anything touches disk — and callers who do want a file call dictk.image.write explicitly, as a separate step. See each function's docstring (rendered in the API reference) for details.

Rosta

We create a synthetic example speckle pattern with the built-in rosta image generator. It implements the Rosta algorithm described by Olufsen (Olufsen SN, Andersen ME, Fagerholt E. muDIC: An open-source toolkit for digital image correlation. SoftwareX. 2020 Jan 1;11:100391, Algorithm 1, page 6, repository).

The help text for rosta:

dictk rosta --help

returns

usage: dictk rosta [-h] [--dot-size DOT_SIZE] [--density DENSITY]
                   [--smoothness SMOOTHNESS] [--random-seed RANDOM_SEED]
                   [--output OUTPUT] [--format {tiff,png,jpg,svg}]
                   [width] [height]

positional arguments:
  width                 Image width in pixels (int), default: 200.
  height                Image height in pixels (int), default: 200.

options:
  -h, --help            show this help message and exit
  --dot-size DOT_SIZE, -s DOT_SIZE
                        Dot pattern size factor, 0.0 to 100.0 (float),
                        default: 4.0.
  --density DENSITY, -d DENSITY
                        Dot pattern density, 0.0 to 1.0 (float), default:
                        0.32.
  --smoothness SMOOTHNESS, -m SMOOTHNESS
                        Smoothness factor, 0.0 to 100.0 (float), default: 2.0.
  --random-seed RANDOM_SEED, -r RANDOM_SEED
                        Seed for reproducible pattern generation (int),
                        default: 42.
  --output OUTPUT, -o OUTPUT
                        Output directory (path), default: current directory.
  --format {tiff,png,jpg,svg}, -f {tiff,png,jpg,svg}
                        Output image format (str), default: tiff.

Create a synthetic image, 200 by 200 pixels, 50% dot density:

dictk rosta 200 200 --density 0.5 --format png -o .

Saved image: rosta_200w_by_200h_dot_4.0_den_0.5_smo_2.0.png

Note that the file name is automatically chosen based on the input parameters.

The result:

The Python equivalent returns the same pixel data as a NumPy array, with no file written:

import dictk

pattern = dictk.rosta(width=200, height=200, density=0.5)

shape=(200, 200), dtype=uint8

Checkerboard

To make it easier to manually identify discrete points in the speckle pattern, dictk can also generate a checkerboard test image.

The help text for checkerboard:

dictk checkerboard --help

returns

usage: dictk checkerboard [-h] [--count-x COUNT_X] [--count-y COUNT_Y]
                          [--output OUTPUT] [--format {tiff,png,jpg,svg}]
                          [width] [height]

positional arguments:
  width                 Image width in pixels (int), default: 200.
  height                Image height in pixels (int), default: 200.

options:
  -h, --help            show this help message and exit
  --count-x COUNT_X, -x COUNT_X
                        Number of rectangles along the width (int), default:
                        8.
  --count-y COUNT_Y, -y COUNT_Y
                        Number of rectangles along the height (int), default:
                        8.
  --output OUTPUT, -o OUTPUT
                        Output directory (path), default: current directory.
  --format {tiff,png,jpg,svg}, -f {tiff,png,jpg,svg}
                        Output image format (str), default: tiff.

Create a synthetic image, 200 by 200 pixels:

dictk checkerboard 200 200 --format png -o .

Saved image: checkerboard_200w_by_200h_8x8.png

The Python equivalent, again returning an array with no file written:

import dictk

board = dictk.checkerboard(width=200, height=200)

shape=(200, 200), dtype=uint8

Astronaut

Unlike rosta and checkerboard, which procedurally generate a fresh synthetic pattern from parameters, astronaut loads a bundled real-world photograph and converts it to grayscale — useful for exercising dictk's imaging utilities against something other than a synthetic pattern. The source is a NASA portrait of astronaut Eileen Collins, from the NASA Great Images database ("No known copyright restrictions, released into the public domain."). Its native resolution is 512x512; passing width/ height other than that resizes the source image rather than generating a new one at that size.

The help text for astronaut:

dictk astronaut --help

returns

usage: dictk astronaut [-h] [--output OUTPUT] [--format {tiff,png,jpg,svg}]
                       [width] [height]

positional arguments:
  width                 Image width in pixels (int), default: 512.
  height                Image height in pixels (int), default: 512.

options:
  -h, --help            show this help message and exit
  --output OUTPUT, -o OUTPUT
                        Output directory (path), default: current directory.
  --format {tiff,png,jpg,svg}, -f {tiff,png,jpg,svg}
                        Output image format (str), default: tiff.

Save it at 300 by 300 pixels — smaller downscales from the native 512x512 start to lose too much detail:

dictk astronaut 300 300 --format png -o .

Saved image: astronaut_300w_by_300h.png

The Python equivalent, again returning an array with no file written:

import dictk

photo = dictk.astronaut(width=300, height=300)

shape=(300, 300), dtype=uint8

Combining into a reference image

combine works on any two grayscale images of the same shape, so it isn't limited to combining the two synthetic images below — Speckle + Astronaut further down combines rosta with a real photograph instead.

Speckle + Checkerboard

We combine the rosta speckle pattern with the checkerboard into a reference image checkerboard0 by averaging their pixel values and normalizing back to uint8:

from dictk.image import combine, read, write

speckle = read(path="rosta_200w_by_200h_dot_4.0_den_0.5_smo_2.0.png")
checker = read(path="checkerboard_200w_by_200h_8x8.png")
checkerboard0 = combine(a=speckle, b=checker)
write(arr=checkerboard0, path="checkerboard0.png")

Saved image: checkerboard0.png

Because both inputs are averaged and rescaled together, the checkerboard's squares stay clearly black or white while the speckle pattern shows up as gray texture within them:

• Where the checkerboard is black, speckle white maps to gray and speckle black stays black.
• Where the checkerboard is white, speckle black maps to gray and speckle white stays white.

That trimodal structure is visible in the pixel-intensity histograms below: speckle and checkerboard are both roughly bimodal (dark/light), while checkerboard0 picks up a distinct middle hump from the black/white-speckle-on-opposite checkerboard combinations.

from dictk.image import read
from dictk.plot import histogram_save

speckle = read(path="rosta_200w_by_200h_dot_4.0_den_0.5_smo_2.0.png")
checker = read(path="checkerboard_200w_by_200h_8x8.png")
checkerboard0 = read(path="checkerboard0.png")

histogram_save(arr=speckle, path="rosta_histogram.png")
histogram_save(arr=checker, path="checkerboard_histogram.png")
histogram_save(arr=checkerboard0, path="checkerboard0_histogram.png")

Saved histograms: rosta_histogram.png, checkerboard_histogram.png, checkerboard0_histogram.png

rosta	checkerboard	checkerboard0

Speckle + Astronaut

The checkerboard above is a stand-in for an actual specimen — in a real DIC setup, the speckle pattern is applied directly to the surface being measured, not swapped in from another generator. Combining rosta with the astronaut photo instead of the checkerboard is closer to that: a speckle pattern overlaid on a realistic, non-uniform grayscale image.

This time the two source images are never written to disk at all — both dictk.rosta and dictk.astronaut return arrays directly, which combine accepts as-is, so only the combined result astronaut0 is saved:

import dictk
from dictk.image import combine, write

speckle = dictk.rosta(width=300, height=300, density=0.5)
photo = dictk.astronaut(width=300, height=300)
astronaut0 = combine(a=speckle, b=photo)
write(arr=astronaut0, path="astronaut0.png")

Saved image: astronaut0.png

Both checkerboard0.png and astronaut0.png are also bundled in src/dictk/data/, alongside the source astronaut.png, so later examples can reuse them without regenerating from scratch each time.

Subimage Generation

Now we consider extracting a subimage from a subject image. A subimage is a useful precursor to image computation: narrowing the data down to a region of interest makes the computation more efficient than considering the full image.

Reference Frames

When we have just a single image, we have a single, trivial reference frame $\mathcal{F}$: origin $O$ (blue dot) at the top-left corner, with the $x$-axis (red) running left-to-right and the $y$-axis (green) running top-to-bottom. origin, and every pixel coordinate on this page, is expressed in this frame — one that's always implicitly present, even in the left panel below where nothing is drawn to show it. astronaut0 here is the same reference image created in Image Generation:

from dictk.image import read
from dictk.plot import reference_frame_plot

astronaut0 = read(path="astronaut0.png")
reference_frame_plot(image=astronaut0, path="reference_frame.png")

Saved: reference_frame.png

When we extract a subimage from an image, it is useful to be explicit about reference frames: the subimage has its own frame $\mathcal{G}$, located within the image's frame $\mathcal{F}$. The Python API section below demonstrates this concept.

Python API

dictk.image.subimage extracts a rectangular crop from a source image: a width x height region whose top-left corner sits at origin. origin may place the requested region partially or completely outside the source image — rather than raising an error, subimage fills whatever doesn't overlap with black (zero) pixels, so the result is always a well-formed height x width array. This is the building block later tutorials use to pull a kernel or search area out of a larger reference/current image pair around a point of interest.

dictk.image.PixelCoordinate is a simple (x, y) NamedTuple used for origin. dictk.image.subimage itself returns the cropped array directly, with no file written.

The examples below use subimage_comparison_plot, which saves a two-panel figure: the left panel shows where the region falls relative to the source image (blue/red boxes), and the right panel shows the extracted result on its own, in its own local frame $\mathcal{G}$ — sharing the same axis limits as the left panel so the two red boxes render at matching scale. It's built from two smaller single-panel functions, also available individually: subimage_bounds_plot (the left panel alone) and subimage_plot (the right panel alone, but zoomed to the subimage's own size rather than sharing the source image's scale).

Square, fully inside

An 80x80 square region entirely within astronaut0's 300x300 bounds. subimage_comparison_plot draws both panels side by side, sharing the same axis limits, so the red box in the right panel renders at identical scale to the one on the left.

from dictk.image import PixelCoordinate
from dictk.plot import subimage_comparison_plot

origin = PixelCoordinate(x=100, y=40)
subimage_comparison_plot(image=astronaut0, origin=origin, width=80, height=80, path="subimage_comparison_80w_by_80h_at_100_40.png")

Saved: subimage_comparison_80w_by_80h_at_100_40.png

Rectangle, fully inside

A 180x70 region — wider than it is tall — also entirely within the source image bounds:

from dictk.image import PixelCoordinate
from dictk.plot import subimage_comparison_plot

origin = PixelCoordinate(x=50, y=200)
subimage_comparison_plot(image=astronaut0, origin=origin, width=180, height=70, path="subimage_comparison_180w_by_70h_at_50_200.png")

Saved: subimage_comparison_180w_by_70h_at_50_200.png

Partially outside

A 120x120 region with a negative origin, straddling the source image's top-left corner. subimage fills the part of the region above and to the left of the source with black:

from dictk.image import PixelCoordinate
from dictk.plot import subimage_comparison_plot

origin = PixelCoordinate(x=-20, y=-40)
subimage_comparison_plot(image=astronaut0, origin=origin, width=120, height=120, path="subimage_comparison_120w_by_120h_at_-20_-40.png")

Saved: subimage_comparison_120w_by_120h_at_-20_-40.png

Completely outside

A 40x100 region entirely beyond the source image's bounds — its x-range (310 to 350) shares no pixels with the source's (0 to 300), so there is no overlap at all and the result is entirely black:

from dictk.image import PixelCoordinate
from dictk.plot import subimage_comparison_plot

origin = PixelCoordinate(x=310, y=250)
subimage_comparison_plot(image=astronaut0, origin=origin, width=40, height=100, path="subimage_comparison_40w_by_100h_at_310_250.png")

Saved: subimage_comparison_40w_by_100h_at_310_250.png

Image Preprocessing

Certain preprocessing steps can make digital image correlation more robust to differences between a reference and current image — some in pixel intensity, others in how well an image's content suits an FFT-based technique such as phase correlation.

• Brightness and Contrast shift and stretch the pixel-intensity histogram, the two most basic differences a reference and current image can have.
• Windowing tapers an image's edges toward zero before a Fourier transform, reducing spectral leakage.

Brightness and Contrast

This page covers two preprocessing steps that can make digital image correlation more robust to pixel-intensity differences between a reference and deformed image, using the astronaut reference image from Image Generation as an example.

Brightness

Brightness shifts the entire pixel-intensity histogram up or down by a constant amount — the whole image gets lighter or darker together, dark areas included. Pushed too far, dark regions wash out to a flat gray and highlights clip at pure white (255), permanently losing detail.

import dictk
from dictk.image import brightness, write
from dictk.plot import histogram_save

photo = dictk.astronaut(width=300, height=300)
write(arr=photo, path="astronaut_original.png")
histogram_save(arr=photo, path="astronaut_original_histogram.png")

bright_1_5 = brightness(arr=photo, factor=1.5)
write(arr=bright_1_5, path="astronaut_brightness_1.5.png")
histogram_save(arr=bright_1_5, path="astronaut_brightness_1.5_histogram.png")

bright_2_0 = brightness(arr=photo, factor=2.0)
write(arr=bright_2_0, path="astronaut_brightness_2.0.png")
histogram_save(arr=bright_2_0, path="astronaut_brightness_2.0_histogram.png")

Saved: astronaut_original.png, astronaut_brightness_1.5.png, astronaut_brightness_2.0.png

factor=1.0 (original)	factor=1.5	factor=2.0

factor=1.0 (original)	factor=1.5	factor=2.0

At factor=1.5 the histogram shifts right as a whole — midtones move into the brighter half and the mean climbs, with a few highlights starting to clip at 255. At factor=2.0 the shift is large enough that a big share of pixels pile up at that 255 ceiling, visible as a tall spike at the histogram's right edge: real detail that's been clipped away and can't be recovered.

Contrast

Contrast is the spread between an image's darkest and lightest pixels. Increasing contrast stretches the histogram outward from its own mean — darks get darker, lights get lighter — while the mean itself stays roughly where it was.

import dictk
from dictk.image import contrast, write
from dictk.plot import histogram_save

photo = dictk.astronaut(width=300, height=300)

contrast_1_5 = contrast(arr=photo, factor=1.5)
write(arr=contrast_1_5, path="astronaut_contrast_1.5.png")
histogram_save(arr=contrast_1_5, path="astronaut_contrast_1.5_histogram.png")

contrast_2_0 = contrast(arr=photo, factor=2.0)
write(arr=contrast_2_0, path="astronaut_contrast_2.0.png")
histogram_save(arr=contrast_2_0, path="astronaut_contrast_2.0_histogram.png")

Saved: astronaut_contrast_1.5.png, astronaut_contrast_2.0.png

factor=1.0 (original)	factor=1.5	factor=2.0

factor=1.0 (original)	factor=1.5	factor=2.0

At factor=1.5 the histogram spreads outward from the mean rather than shifting — the astronaut's silhouette and helmet edges get sharper, while the mean barely moves. At factor=2.0 the spread is wide enough that more pixels pile up at both the 0 and 255 ends, crushing fine midtone detail even as high-contrast edges sharpen further.

Key Insight: Contrast stretches the histogram, while brightness translates it.

Windowing

The FFT implicitly treats an image as one period of an infinitely-repeating signal. If the content doesn't tile seamlessly, which is the general case since nothing arranges an image's edges to match up, that discontinuity leaks energy across many frequencies rather than the few the underlying content actually has, an effect called spectral leakage. In a correlation surface, leakage broadens and can shift the peak, hurting the precision of any technique that searches that surface for a match.

Windowing counters this by tapering an image's edges toward zero before transforming it, so the (still discontinuous, but now near-zero) seam contributes far less energy. Two standard 1D windows, applied to an image by taking the outer product of a window with itself along each axis:

$$w_{\mathrm{Hann}}(n)=0.5\left(1-\cos \left(\frac{2\pi n}{N-1}\right)\right)$$

$$w_{\mathrm{Hamming}}(n)=0.54-0.46\cos \left(\frac{2\pi n}{N-1}\right)$$

for $n=0,\dots ,N-1$ across a window of length $N$. Hann tapers all the way to exactly zero at both ends; Hamming stops short (around $0.08$), trading a little residual discontinuity for a narrower main lobe in the transformed signal.

See Harris FJ. "On the use of windows for harmonic analysis with the discrete Fourier transform." Proceedings of the IEEE 1978;66(1):51-83. A U.S. government work, not protected by U.S. copyright.

window()

dictk.correlation.window applies either taper to a 2D array. This reuses kernel from Cross Correlation (CC) and the Fourier Domain section of Correlation Criteria — the same checkerboard0, p0, and kernel_margin — to show what tapering actually does to an image before it's passed to an FFT:

import numpy as np
import matplotlib.pyplot as plt
from dictk.image import read, PixelCoordinate, subimage, write
from dictk.correlation import window, WindowingMethod

reference_image = read(path="checkerboard0.png")
p0 = PixelCoordinate(x=100, y=75)
kernel_margin = 25
kernel = subimage(
    image=reference_image,
    origin=PixelCoordinate(x=p0.x - kernel_margin, y=p0.y - kernel_margin),
    width=2 * kernel_margin,
    height=2 * kernel_margin,
)
write(arr=kernel, path="windowing_kernel_original.png")

kernel_hann = window(arr=kernel, method=WindowingMethod.HANN)
write(arr=kernel_hann.astype(np.uint8), path="windowing_kernel_hann.png")

kernel_hamming = window(arr=kernel, method=WindowingMethod.HAMMING)
write(arr=kernel_hamming.astype(np.uint8), path="windowing_kernel_hamming.png")

# window()'s own weights, isolated from kernel's content: windowing an
# all-ones array leaves exactly the 2D weight array behind. A single row
# at the kernel's mid-height cuts through the row axis's own peak (~1.0),
# so what's left is each method's column-axis taper alone.
mid_row = kernel.shape[0] // 2
ones = np.ones_like(kernel, dtype=np.float64)
weight_profiles = {
    "none": np.ones(kernel.shape[1]),
    "hann": window(arr=ones, method=WindowingMethod.HANN)[mid_row, :],
    "hamming": window(arr=ones, method=WindowingMethod.HAMMING)[mid_row, :],
}
for name, profile in weight_profiles.items():
    fig, ax = plt.subplots(figsize=(4, 2.5), constrained_layout=True)
    ax.plot(profile, color="black")
    ax.set_ylim(-0.05, 1.05)  # shared across all three, for a fair comparison
    ax.set_xlabel("x (pixels)")
    ax.set_ylabel("window weight")
    fig.savefig(f"windowing_kernel_cut_{name}.png", dpi=300)
    plt.close(fig)

Saved: windowing_kernel_original.png, windowing_kernel_hann.png, windowing_kernel_hamming.png, windowing_kernel_cut_none.png, windowing_kernel_cut_hann.png, windowing_kernel_cut_hamming.png

none	Hann	Hamming

Every edge fades toward black; Hann's corners go fully black (tapers to exactly 0), while Hamming's stay a faint gray (tapers to $0.08\times 0.08\approx 0.006$ of the original corner pixel, the product of both axes' own $0.08$ edge value).

The bottom row makes each method's own taper precise, independent of checkerboard0's content: a horizontal cut through the window's weight array at the kernel's mid-height, all three sharing the same $y$-axis. none is flat at $1.0$ everywhere -- no taper at all. Hann and Hamming both peak at $1.0$ at that same mid-height (the row axis's own window is near its own peak there), so this cut isolates the column axis's taper alone: Hann reaches exactly $0$ at both edges, Hamming levels off at $0.08$ -- not the smaller $0.006$ corner value above, since a corner is where both axes are simultaneously at their own edge, and a mid-height cut only ever passes through one axis's edge at a time.

See Correlation Visualization for windowing shown in action, tapering a real kernel and search area before they're compared.

Image Transformation

Image deformations, also called transformations in the computer vision literature (see Szeliski¹), fall into the categories shown below:

Each category preserves a different, nested set of geometric properties — every property a category preserves is also preserved by every more restrictive category to its left:

Pure Translation (Rigid Body Motion)

As the simplest of the categories above — no change in shape or size — dictk.image.translate shifts every pixel by a fixed displacement. This example shifts the image by dx=-60 pixels in x and dy=+80 pixels in y, representing rigid-body motion where the material moves without deforming.

import dictk
from dictk.image import translate, write

photo = dictk.astronaut(width=300, height=300)
write(arr=photo, path="astronaut_translate_original.png")

translated = translate(arr=photo, dx=-60, dy=80)
write(arr=translated, path="astronaut_translate_rigid_body.png")

Saved: astronaut_translate_original.png, astronaut_translate_rigid_body.png

Translation	Image
Original
dx=-60, dy=+80

Pure Rotation

A 30° counterclockwise rotation, another rigid-body motion that preserves distances and angles. dictk.image.rotate pivots on the image's top-left corner (0, 0), consistent with stretch and translate's pivot choice in this codebase — unlike the more typical "object spins in place" rotation about the center, most content swings away from that fixed corner, similar to a door on a hinge.

import dictk
from dictk.image import rotate, write

photo = dictk.astronaut(width=300, height=300)
write(arr=photo, path="astronaut_rotate_original.png")

rotated = rotate(arr=photo, angle=30.0)
write(arr=rotated, path="astronaut_rotate_30deg.png")

Saved: astronaut_rotate_original.png, astronaut_rotate_30deg.png

Rotation	Image
Original
30° (origin-pivoted)

X-Axis Stretch (Extension)

As a concrete example of the similarity category above, dictk.image.stretch applies a uniaxial stretch along the x-axis: the image's top-left corner (x=0, y=0) stays fixed, and content grows away from it, using backward mapping with bilinear interpolation so the result has no gaps (unlike naively moving each source pixel forward, which can leave holes). The two stretches below range from a small, realistic deformation (5%, similar in magnitude to a modest tensile strain in a materials test) up to a much larger one (50%).

import dictk
from dictk.image import stretch, write

photo = dictk.astronaut(width=300, height=300)
write(arr=photo, path="astronaut_stretch_original.png")

stretch_5pct = stretch(arr=photo, factor_x=1.05)
write(arr=stretch_5pct, path="astronaut_stretch_x_5pct.png")

stretch_50pct = stretch(arr=photo, factor_x=1.50)
write(arr=stretch_50pct, path="astronaut_stretch_x_50pct.png")

Saved: astronaut_stretch_original.png, astronaut_stretch_x_5pct.png, astronaut_stretch_x_50pct.png

Stretch	Image
Original
5% (factor_x=1.05)
50% (factor_x=1.50)