Merge pull request #267 from ornlneutronimaging/541_rotation_center

   git submodule init

   git submodule update

Useful Functions

----------------

The [cache](https://docs.python.org/3/library/functools.html#functools.cache) decorator can be used to hold onto returns of functions.

The data is stored in a ``dict`` with the function parameters as the key, and the return as the value.

The cache remains in memory for the lifetime of the pytest run and a [pytest fixture](https://docs.pytest.org/en/7.1.x/how-to/fixtures.html) with more limited scope may be more appropriate.

Access Development Version on Analysis Cluster

----------------------------------------------

        whether angles are in degrees or radians, default is True (degrees)

    atol_deg: float = 1e-3

        tolerance for the search of 180 deg paris, default is 0.1 degrees

    num_pairs: int = 1

        Number of pairs to look for. Specifying -1 means as many pairs as possible.

    max_workers: int = 0

        number of cores to use for parallel median filtering, default is 0, which means using all available cores.

    tqdm_class: imars3d.ui.widgets.TqdmType

        default=1e-3,

        doc="tolerance for the search of 180 deg paris, default is 0.1 degrees",

    )

    num_pairs = param.Integer(

        default=1, bounds=(-1, None), doc="Number of pairs to look for. Specifying -1 means as many pairs as possible."

    )

    max_workers = param.Integer(

        default=0,

        bounds=(0, None),

        self.max_workers = clamp_max_workers(params.max_workers)

        val = self._find_rotation_center(

            params.arrays, params.angles, params.in_degrees, params.atol_deg, self.max_workers, params.tqdm_class

            arrays=params.arrays,

            angles=params.angles,

            in_degrees=params.in_degrees,

            atol_deg=params.atol_deg,

            num_pairs=params.num_pairs,

            max_workers=self.max_workers,

            tqdm_class=params.tqdm_class,

        )

        logger.info("FINISHED Executing Filter: Find Rotation Center")

        return val

        angles: np.ndarray,

        in_degrees: bool = True,

        atol_deg: float = 1e-3,

        num_pairs: int = 1,

        max_workers: int = -1,

        tqdm_class=None,

    ) -> float:

        # locate 180 degree pairs

        atol = atol_deg if in_degrees else np.radians(atol_deg)

        idx_low, idx_hgh = find_180_deg_pairs_idx(angles, atol=atol, in_degrees=in_degrees)

        if num_pairs <= 0 or num_pairs >= idx_low.size:

            logger.info("Using all pairs of angles")

        elif num_pairs == 1:

            idx_low = [idx_low[0]]

            idx_hgh = [idx_hgh[0]]

            logger.info("Using one pair of angles")

        else:

            # integer division to get correct size if possible

            span = idx_low.size // num_pairs

            # get equally spaced items if possible

            if span > 1:

                idx_low = idx_low[::span]

                idx_hgh = idx_hgh[::span]

            # trim down to the requested number

            # the selected angels are not equally spaced

            if idx_low.size > num_pairs:

                idx_low = idx_low[:num_pairs]

                idx_hgh = idx_hgh[:num_pairs]

            logger.info(f"Using {idx_low.size} pairs of angles")

        # process

        max_workers = clamp_max_workers(max_workers)

        # use shared memory model and tqdm wrapper for multiprocessing to reduce

# standard imports

from datetime import datetime

import logging

import multiprocessing

import resource

from typing import Union

logger = logging.getLogger(__name__)

def clamp_max_workers(max_workers: Union[int, None]) -> int:

    """Calculate the number of max workers.

    """

    if max_workers is None:

        max_workers = 0

    return min(resource.RLIMIT_NPROC, max(1, multiprocessing.cpu_count() - 2)) if max_workers <= 0 else max_workers

    result = min(resource.RLIMIT_NPROC, max(1, multiprocessing.cpu_count() - 2)) if max_workers <= 0 else max_workers

    # log if the value was different than what was asked for

    if max_workers == 0 and result != max_workers:

        logger.info(f"Due to system load, setting maximum workers to {result}")

    return result

def to_time_str(value: datetime = datetime.now()) -> str:

#!/usr/bin/env python3

from functools import cache

import numpy as np

import pytest

import tomopy

from imars3d.backend.diagnostics.rotation import find_rotation_center

# all tests share a consistent set of omeaga angles

# these are in radians

OMEGAS = np.linspace(0, np.pi * 2, 181)

# decorator creates a dict of previous parameter/return pairs

# the cache does not empty

# each call takes ~1s

# Time difference on 2022-12-09 was 17s vs 11s

@cache

def get_synthetic_stack(

    center: float,

    omegas: np.ndarray,

) -> np.ndarray:

    """

    Generate synthetic radiograph stack with given rotation center and rotation angles.

    ---------

    @param center:

        Rotation center

    @param omegas:

        Rotation angles in radians

    Return

    ------

    shepp3d = tomopy.misc.phantom.shepp3d(size=129)

    projs = tomopy.sim.project.project(

        shepp3d,

        omegas,

        OMEGAS,

        emission=False,

        center=center,

    )

    return projs

@pytest.mark.parametrize(

    "num_pairs",

    [-1, 0, 1, 44, 45, 46, 200],

)

def test_pairs(num_pairs):

    CENTER_REF = 80.5

    projs = get_synthetic_stack(CENTER_REF)

    center_calc = find_rotation_center(arrays=projs, angles=OMEGAS, in_degrees=False, num_pairs=num_pairs)

    # answer within the same pixel should be sufficient for most cases

    np.testing.assert_almost_equal(center_calc, CENTER_REF, decimal=1)

@pytest.mark.parametrize(

    "center_ref",

    [

        119.5,

    ],

)

def test_find_rotation_center(center_ref):

    #

    # case 0: valid input

def test_differrent_centers(center_ref):

    # NOTE: unit of omegas is handled by find_180_deg_pairs_idx

    omegas = np.linspace(0, np.pi * 2, 181)

    projs = get_synthetic_stack(center_ref, omegas)

    center_calc = find_rotation_center(arrays=projs, angles=omegas, in_degrees=False)

    projs = get_synthetic_stack(center_ref)

    # this is using default number of pairs (1)

    center_calc = find_rotation_center(arrays=projs, angles=OMEGAS, in_degrees=False)

    # verify

    # NOTE:

    # answer within the same pixel should be sufficient for most cases

    np.testing.assert_almost_equal(center_calc, center_ref, decimal=1)

    #

    # case 1: wrong dimension

    np.testing.assert_allclose(center_calc, center_ref, atol=0.2)

def test_wrong_dimension():

    projs = np.random.random(100).reshape(10, 10)

    with pytest.raises(ValueError):

        find_rotation_center(arrays=projs, angles=omegas)

        find_rotation_center(arrays=projs, angles=OMEGAS)

if __name__ == "__main__":