Enhancement: Average affines using Log-Euclidean Mean from lie group algebra?

[gdevenyi]

Is your feature request related to a problem? Please describe. Currently, as far as I can tell, transforms are averaged by simply averaging the individual serialized parameters.

I discovered in the MINC family of scripts a command which is also for averaging affine transforms: https://github.com/BIC-MNI/minc-widgets/blob/master/xfmavg/xfmavg#L57-L66

Which includes this note about the implementation details:

| This is done via Matrix logs and exponents (currently forks octave)
|
|    As affine transformations are in the 
|    Lie group of matrices, the premise is:
|       average(t1,..,tn) = mexp[ [mlog(t1) + .. + mlog(tn)] / n ]
|
|    In MATLAB/Octave speak this translates to:
|       AVG = expm((logm(t1) + ... + logm(tn)) / n)

Describe the solution you'd like Use lie group algebra concepts for the average

Describe alternatives you've considered Don't

Additional context This paper appears to argue that the Log-Euclidean Mean does not satisfy all the requirements of the mean in the lie group https://hal.inria.fr/hal-00938320/document (Section 3.1)

Despite its intuitive formulation, the Log-Euclidean mean is not admissible: it is not consistent with the left and right translations when the Lie group is not abelian

But I'm not sure of the subset of the lie group that is affine transformations avoids these restrictions as I don't have enough linear algebra background.

[gdevenyi]

I wrote a implementation in SimpleITK

#!/usr/bin/env python

# Based on https://github.com/BIC-MNI/pyezminc/blob/develop/examples/xfmavg_scipy.py

import SimpleITK as sitk
import numpy as np
import argparse

# needed for matrix log and exp
import scipy.linalg

def itk_to_homogeneous_matrix(itk_transform):
    # Dump the matrix 4x3 matrix
    input_itk_transform_parameters = itk_transform.GetParameters()

    M = np.zeros((4, 4))
    M[0:3, 0:3] = np.asarray(input_itk_transform_parameters[0:9]).reshape((3, 3))
    M[0:3, 3] = input_itk_transform_parameters[9:]
    M[3, 3] = 1

    return M

if __name__ == "__main__":
    parser = argparse.ArgumentParser(
        formatter_class=argparse.ArgumentDefaultsHelpFormatter
    )
    parser.add_argument(
        "-o",
        "--output",
        type=str,
        required=True,
        help="""
            Name of output average transform.
            """,
    )
    parser.add_argument(
        "--file_list",
        type=str,
        nargs="*",  # 0 or more values expected => creates a list
        required=True,
        help="""
            Specify a list of input files, space-separated (i.e. file1 file2 ...).
            """,
    )
    opts = parser.parse_args()

    average_matrix = np.zeros((4, 4), dtype=complex)

    for file in opts.file_list:
        # Input transform
        input_itk_transform = sitk.ReadTransform(file)
        average_matrix += scipy.linalg.logm(
            itk_to_homogeneous_matrix(input_itk_transform)
        )

    average_matrix /= len(opts.file_list)
    average_matrix = scipy.linalg.expm(average_matrix).real

    output_average_transform = sitk.AffineTransform(average_matrix[0:3, 0:3].flatten(), average_matrix[0:3, 3].flatten(), (0,0,0))

    print(output_average_transform)

    sitk.WriteTransform(output_average_transform, opts.output)