diff --git a/.gitignore b/.gitignore index 2fc8ce3..a81c8ee 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,138 @@ -build -dist -*.egg-info +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py,cover +.hypothesis/ +.pytest_cache/ +cover/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 +db.sqlite3-journal + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +.pybuilder/ +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py + +# pyenv +# For a library or package, you might want to ignore these files since the code is +# intended to run in multiple environments; otherwise, check them in: +# .python-version + +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# Pyre type checker +.pyre/ + +# pytype static type analyzer +.pytype/ + +# Cython debug symbols +cython_debug/ diff --git a/.pylintrc b/.pylintrc new file mode 100644 index 0000000..aa94ea9 --- /dev/null +++ b/.pylintrc @@ -0,0 +1,595 @@ +[MASTER] + +# A comma-separated list of package or module names from where C extensions may +# be loaded. Extensions are loading into the active Python interpreter and may +# run arbitrary code. +extension-pkg-whitelist=fire,numpy,pandas,cv2,skimage,scipy,sklearn,SimpleITK,matplotlib,torch, + +# Specify a score threshold to be exceeded before program exits with error. +fail-under=10 + +# Add files or directories to the blacklist. They should be base names, not +# paths. +ignore=CVS + +# Add files or directories matching the regex patterns to the blacklist. The +# regex matches against base names, not paths. +ignore-patterns= + +# Python code to execute, usually for sys.path manipulation such as +# pygtk.require(). +#init-hook= + +# Use multiple processes to speed up Pylint. Specifying 0 will auto-detect the +# number of processors available to use. +jobs=1 + +# Control the amount of potential inferred values when inferring a single +# object. This can help the performance when dealing with large functions or +# complex, nested conditions. +limit-inference-results=100 + +# List of plugins (as comma separated values of python module names) to load, +# usually to register additional checkers. +load-plugins= + +# Pickle collected data for later comparisons. +persistent=yes + +# When enabled, pylint would attempt to guess common misconfiguration and emit +# user-friendly hints instead of false-positive error messages. +suggestion-mode=yes + +# Allow loading of arbitrary C extensions. Extensions are imported into the +# active Python interpreter and may run arbitrary code. +unsafe-load-any-extension=no + + +[MESSAGES CONTROL] + +# Only show warnings with the listed confidence levels. Leave empty to show +# all. Valid levels: HIGH, INFERENCE, INFERENCE_FAILURE, UNDEFINED. +confidence= + +# Disable the message, report, category or checker with the given id(s). You +# can either give multiple identifiers separated by comma (,) or put this +# option multiple times (only on the command line, not in the configuration +# file where it should appear only once). You can also use "--disable=all" to +# disable everything first and then reenable specific checks. For example, if +# you want to run only the similarities checker, you can use "--disable=all +# --enable=similarities". If you want to run only the classes checker, but have +# no Warning level messages displayed, use "--disable=all --enable=classes +# --disable=W". +disable=print-statement, + parameter-unpacking, + unpacking-in-except, + old-raise-syntax, + backtick, + long-suffix, + old-ne-operator, + old-octal-literal, + import-star-module-level, + non-ascii-bytes-literal, + raw-checker-failed, + bad-inline-option, + locally-disabled, + file-ignored, + suppressed-message, + useless-suppression, + deprecated-pragma, + use-symbolic-message-instead, + apply-builtin, + basestring-builtin, + buffer-builtin, + cmp-builtin, + coerce-builtin, + execfile-builtin, + file-builtin, + long-builtin, + raw_input-builtin, + reduce-builtin, + standarderror-builtin, + unicode-builtin, + xrange-builtin, + coerce-method, + delslice-method, + getslice-method, + setslice-method, + no-absolute-import, + old-division, + dict-iter-method, + dict-view-method, + next-method-called, + metaclass-assignment, + indexing-exception, + raising-string, + reload-builtin, + oct-method, + hex-method, + nonzero-method, + cmp-method, + input-builtin, + round-builtin, + intern-builtin, + unichr-builtin, + map-builtin-not-iterating, + zip-builtin-not-iterating, + range-builtin-not-iterating, + filter-builtin-not-iterating, + using-cmp-argument, + eq-without-hash, + div-method, + idiv-method, + rdiv-method, + exception-message-attribute, + invalid-str-codec, + sys-max-int, + bad-python3-import, + deprecated-string-function, + deprecated-str-translate-call, + deprecated-itertools-function, + deprecated-types-field, + next-method-defined, + dict-items-not-iterating, + dict-keys-not-iterating, + dict-values-not-iterating, + deprecated-operator-function, + deprecated-urllib-function, + xreadlines-attribute, + deprecated-sys-function, + exception-escape, + comprehension-escape + +# Enable the message, report, category or checker with the given id(s). You can +# either give multiple identifier separated by comma (,) or put this option +# multiple time (only on the command line, not in the configuration file where +# it should appear only once). See also the "--disable" option for examples. +enable=c-extension-no-member + + +[REPORTS] + +# Python expression which should return a score less than or equal to 10. You +# have access to the variables 'error', 'warning', 'refactor', and 'convention' +# which contain the number of messages in each category, as well as 'statement' +# which is the total number of statements analyzed. This score is used by the +# global evaluation report (RP0004). +evaluation=10.0 - ((float(5 * error + warning + refactor + convention) / statement) * 10) + +# Template used to display messages. This is a python new-style format string +# used to format the message information. See doc for all details. +#msg-template= + +# Set the output format. Available formats are text, parseable, colorized, json +# and msvs (visual studio). You can also give a reporter class, e.g. +# mypackage.mymodule.MyReporterClass. +output-format=text + +# Tells whether to display a full report or only the messages. +reports=no + +# Activate the evaluation score. +score=yes + + +[REFACTORING] + +# Maximum number of nested blocks for function / method body +max-nested-blocks=5 + +# Complete name of functions that never returns. When checking for +# inconsistent-return-statements if a never returning function is called then +# it will be considered as an explicit return statement and no message will be +# printed. +never-returning-functions=sys.exit + + +[STRING] + +# This flag controls whether inconsistent-quotes generates a warning when the +# character used as a quote delimiter is used inconsistently within a module. +check-quote-consistency=no + +# This flag controls whether the implicit-str-concat should generate a warning +# on implicit string concatenation in sequences defined over several lines. +check-str-concat-over-line-jumps=no + + +[SPELLING] + +# Limits count of emitted suggestions for spelling mistakes. +max-spelling-suggestions=4 + +# Spelling dictionary name. Available dictionaries: none. To make it work, +# install the python-enchant package. +spelling-dict= + +# List of comma separated words that should not be checked. +spelling-ignore-words= + +# A path to a file that contains the private dictionary; one word per line. +spelling-private-dict-file= + +# Tells whether to store unknown words to the private dictionary (see the +# --spelling-private-dict-file option) instead of raising a message. +spelling-store-unknown-words=no + + +[TYPECHECK] + +# List of decorators that produce context managers, such as +# contextlib.contextmanager. Add to this list to register other decorators that +# produce valid context managers. +contextmanager-decorators=contextlib.contextmanager + +# List of members which are set dynamically and missed by pylint inference +# system, and so shouldn't trigger E1101 when accessed. Python regular +# expressions are accepted. +generated-members= + +# Tells whether missing members accessed in mixin class should be ignored. A +# mixin class is detected if its name ends with "mixin" (case insensitive). +ignore-mixin-members=yes + +# Tells whether to warn about missing members when the owner of the attribute +# is inferred to be None. +ignore-none=yes + +# This flag controls whether pylint should warn about no-member and similar +# checks whenever an opaque object is returned when inferring. The inference +# can return multiple potential results while evaluating a Python object, but +# some branches might not be evaluated, which results in partial inference. In +# that case, it might be useful to still emit no-member and other checks for +# the rest of the inferred objects. +ignore-on-opaque-inference=yes + +# List of class names for which member attributes should not be checked (useful +# for classes with dynamically set attributes). This supports the use of +# qualified names. +ignored-classes=optparse.Values,thread._local,_thread._local + +# List of module names for which member attributes should not be checked +# (useful for modules/projects where namespaces are manipulated during runtime +# and thus existing member attributes cannot be deduced by static analysis). It +# supports qualified module names, as well as Unix pattern matching. +ignored-modules= + +# Show a hint with possible names when a member name was not found. The aspect +# of finding the hint is based on edit distance. +missing-member-hint=yes + +# The minimum edit distance a name should have in order to be considered a +# similar match for a missing member name. +missing-member-hint-distance=1 + +# The total number of similar names that should be taken in consideration when +# showing a hint for a missing member. +missing-member-max-choices=1 + +# List of decorators that change the signature of a decorated function. +signature-mutators= + + +[BASIC] + +# Naming style matching correct argument names. +argument-naming-style=snake_case + +# Regular expression matching correct argument names. Overrides argument- +# naming-style. +#argument-rgx= + +# Naming style matching correct attribute names. +attr-naming-style=snake_case + +# Regular expression matching correct attribute names. Overrides attr-naming- +# style. +#attr-rgx= + +# Bad variable names which should always be refused, separated by a comma. +bad-names=foo, + bar, + baz, + toto, + tutu, + tata + +# Bad variable names regexes, separated by a comma. If names match any regex, +# they will always be refused +bad-names-rgxs= + +# Naming style matching correct class attribute names. +class-attribute-naming-style=any + +# Regular expression matching correct class attribute names. Overrides class- +# attribute-naming-style. +#class-attribute-rgx= + +# Naming style matching correct class names. +class-naming-style=PascalCase + +# Regular expression matching correct class names. Overrides class-naming- +# style. +#class-rgx= + +# Naming style matching correct constant names. +const-naming-style=UPPER_CASE + +# Regular expression matching correct constant names. Overrides const-naming- +# style. +#const-rgx= + +# Minimum line length for functions/classes that require docstrings, shorter +# ones are exempt. +docstring-min-length=-1 + +# Naming style matching correct function names. +function-naming-style=snake_case + +# Regular expression matching correct function names. Overrides function- +# naming-style. +#function-rgx= + +# Good variable names which should always be accepted, separated by a comma. +good-names=i, + j, + k, + ex, + Run, + _ + +# Good variable names regexes, separated by a comma. If names match any regex, +# they will always be accepted +good-names-rgxs= + +# Include a hint for the correct naming format with invalid-name. +include-naming-hint=no + +# Naming style matching correct inline iteration names. +inlinevar-naming-style=any + +# Regular expression matching correct inline iteration names. Overrides +# inlinevar-naming-style. +#inlinevar-rgx= + +# Naming style matching correct method names. +method-naming-style=snake_case + +# Regular expression matching correct method names. Overrides method-naming- +# style. +#method-rgx= + +# Naming style matching correct module names. +module-naming-style=snake_case + +# Regular expression matching correct module names. Overrides module-naming- +# style. +#module-rgx= + +# Colon-delimited sets of names that determine each other's naming style when +# the name regexes allow several styles. +name-group= + +# Regular expression which should only match function or class names that do +# not require a docstring. +no-docstring-rgx=^_ + +# List of decorators that produce properties, such as abc.abstractproperty. Add +# to this list to register other decorators that produce valid properties. +# These decorators are taken in consideration only for invalid-name. +property-classes=abc.abstractproperty + +# Naming style matching correct variable names. +variable-naming-style=snake_case + +# Regular expression matching correct variable names. Overrides variable- +# naming-style. +#variable-rgx= + + +[SIMILARITIES] + +# Ignore comments when computing similarities. +ignore-comments=yes + +# Ignore docstrings when computing similarities. +ignore-docstrings=yes + +# Ignore imports when computing similarities. +ignore-imports=no + +# Minimum lines number of a similarity. +min-similarity-lines=4 + + +[MISCELLANEOUS] + +# List of note tags to take in consideration, separated by a comma. +notes=FIXME, + XXX, + TODO + +# Regular expression of note tags to take in consideration. +#notes-rgx= + + +[LOGGING] + +# The type of string formatting that logging methods do. `old` means using % +# formatting, `new` is for `{}` formatting. +logging-format-style=old + +# Logging modules to check that the string format arguments are in logging +# function parameter format. +logging-modules=logging + + +[VARIABLES] + +# List of additional names supposed to be defined in builtins. Remember that +# you should avoid defining new builtins when possible. +additional-builtins= + +# Tells whether unused global variables should be treated as a violation. +allow-global-unused-variables=yes + +# List of strings which can identify a callback function by name. A callback +# name must start or end with one of those strings. +callbacks=cb_, + _cb + +# A regular expression matching the name of dummy variables (i.e. expected to +# not be used). +dummy-variables-rgx=_+$|(_[a-zA-Z0-9_]*[a-zA-Z0-9]+?$)|dummy|^ignored_|^unused_ + +# Argument names that match this expression will be ignored. Default to name +# with leading underscore. +ignored-argument-names=_.*|^ignored_|^unused_ + +# Tells whether we should check for unused import in __init__ files. +init-import=no + +# List of qualified module names which can have objects that can redefine +# builtins. +redefining-builtins-modules=six.moves,past.builtins,future.builtins,builtins,io + + +[FORMAT] + +# Expected format of line ending, e.g. empty (any line ending), LF or CRLF. +expected-line-ending-format= + +# Regexp for a line that is allowed to be longer than the limit. +ignore-long-lines=^\s*(# )??$ + +# Number of spaces of indent required inside a hanging or continued line. +indent-after-paren=4 + +# String used as indentation unit. This is usually " " (4 spaces) or "\t" (1 +# tab). +indent-string=' ' + +# Maximum number of characters on a single line. +max-line-length=100 + +# Maximum number of lines in a module. +max-module-lines=1000 + +# List of optional constructs for which whitespace checking is disabled. `dict- +# separator` is used to allow tabulation in dicts, etc.: {1 : 1,\n222: 2}. +# `trailing-comma` allows a space between comma and closing bracket: (a, ). +# `empty-line` allows space-only lines. +no-space-check=trailing-comma, + dict-separator + +# Allow the body of a class to be on the same line as the declaration if body +# contains single statement. +single-line-class-stmt=no + +# Allow the body of an if to be on the same line as the test if there is no +# else. +single-line-if-stmt=no + + +[IMPORTS] + +# List of modules that can be imported at any level, not just the top level +# one. +allow-any-import-level= + +# Allow wildcard imports from modules that define __all__. +allow-wildcard-with-all=no + +# Analyse import fallback blocks. This can be used to support both Python 2 and +# 3 compatible code, which means that the block might have code that exists +# only in one or another interpreter, leading to false positives when analysed. +analyse-fallback-blocks=no + +# Deprecated modules which should not be used, separated by a comma. +deprecated-modules=optparse,tkinter.tix + +# Create a graph of external dependencies in the given file (report RP0402 must +# not be disabled). +ext-import-graph= + +# Create a graph of every (i.e. internal and external) dependencies in the +# given file (report RP0402 must not be disabled). +import-graph= + +# Create a graph of internal dependencies in the given file (report RP0402 must +# not be disabled). +int-import-graph= + +# Force import order to recognize a module as part of the standard +# compatibility libraries. +known-standard-library= + +# Force import order to recognize a module as part of a third party library. +known-third-party=enchant + +# Couples of modules and preferred modules, separated by a comma. +preferred-modules= + + +[DESIGN] + +# Maximum number of arguments for function / method. +max-args=5 + +# Maximum number of attributes for a class (see R0902). +max-attributes=7 + +# Maximum number of boolean expressions in an if statement (see R0916). +max-bool-expr=5 + +# Maximum number of branch for function / method body. +max-branches=12 + +# Maximum number of locals for function / method body. +max-locals=15 + +# Maximum number of parents for a class (see R0901). +max-parents=7 + +# Maximum number of public methods for a class (see R0904). +max-public-methods=20 + +# Maximum number of return / yield for function / method body. +max-returns=6 + +# Maximum number of statements in function / method body. +max-statements=50 + +# Minimum number of public methods for a class (see R0903). +min-public-methods=2 + + +[CLASSES] + +# List of method names used to declare (i.e. assign) instance attributes. +defining-attr-methods=__init__, + __new__, + setUp, + __post_init__ + +# List of member names, which should be excluded from the protected access +# warning. +exclude-protected=_asdict, + _fields, + _replace, + _source, + _make + +# List of valid names for the first argument in a class method. +valid-classmethod-first-arg=cls + +# List of valid names for the first argument in a metaclass class method. +valid-metaclass-classmethod-first-arg=cls + + +[EXCEPTIONS] + +# Exceptions that will emit a warning when being caught. Defaults to +# "BaseException, Exception". +overgeneral-exceptions=BaseException, + Exception diff --git a/bin/install_mixmatch_dependencies b/bin/install_mixmatch_dependencies new file mode 100644 index 0000000..b977c0c --- /dev/null +++ b/bin/install_mixmatch_dependencies @@ -0,0 +1,3 @@ +#!/bin/bash +pip install git+https://github.com/jlevy44/airlab.git@no_mixed_precision --upgrade +pip install git+https://github.com/jlevy44/PathFlow-MixMatch@fix_nonlinear_memory_issues --upgrade diff --git a/pathflow_mixmatch/cli.py b/pathflow_mixmatch/cli.py index 9943594..300178a 100644 --- a/pathflow_mixmatch/cli.py +++ b/pathflow_mixmatch/cli.py @@ -1,5 +1,5 @@ import fire -import numpy as np, cv2 +import numpy as np, pandas as pd, cv2 import cv2 from skimage.morphology import watershed from skimage.feature import peak_local_max @@ -28,46 +28,86 @@ import matplotlib.pyplot as plt from contextlib import contextmanager import sys, os +import pickle +from skimage.exposure import rescale_intensity +from skimage.color import rgb2hed +from airlab.utils.image import Image + + + +# from apex import amp + @contextmanager def suppress_stdout(): - with open(os.devnull, "w") as devnull: - old_stdout = sys.stdout - sys.stdout = devnull - try: - yield - finally: - sys.stdout = old_stdout + with open(os.devnull, "w") as devnull: + old_stdout = sys.stdout + sys.stdout = devnull + try: + yield + finally: + sys.stdout = old_stdout def label_objects(I, min_object_size, threshold=220, connectivity=8, kernel=8, apply_watershed=False): - #try: - BW = (cv2.cvtColor(I, cv2.COLOR_BGR2GRAY)0.95 + x=PCA(n_components=4,random_state=42).fit_transform(pd.concat((props,props2)))#n_components=2 StandardScaler().fit_transform() + return enumerate(pd.DataFrame(sklearn.metrics.pairwise.pairwise_distances(x[:props.shape[0]],x[props.shape[0]:],'euclidean')).values.argmin(1))#pd.DataFrame(sklearn.metrics.pairwise.cosine_similarity(pd.DataFrame(x).T)).iloc[:props.shape[0],props.shape[0]:]#pd.DataFrame(x).T.corr()>0.95 def displace_image(img, displacement, gpu_device, dtype=th.float32): - channels=[] - for i in range(3): - im=sitk.GetImageFromArray(img[...,i]) - im=al.utils.image.create_tensor_image_from_itk_image(im, dtype=dtype, device=('cuda:{}'.format(gpu_device) if gpu_device>=0 else 'cpu')) - channels.append(al.transformation.utils.warp_image(im, displacement).numpy()) - return np.uint8(np.stack(channels).transpose((1,2,0))) + # channels=[] + image=al.image_from_numpy(np.squeeze(img),(),(), dtype=dtype, device=('cuda:{}'.format(gpu_device) if gpu_device>=0 else 'cpu'))#[...,i] + image_size = image.size + print(image_size) + grid = al.transformation.utils.compute_grid(image_size[:2], dtype=image.dtype, device=image.device) + + displacement, grid, _, _ = match_image_size(np.squeeze(displacement.cpu().numpy()), np.squeeze(grid.cpu().numpy())) + if displacement.ndim == 2: + displacement = displacement[:, :, np.newaxis] + # elif displacement.ndim == 3: + # displacement = displacement[:, np.newaxis] + flow_field_grid = displacement + grid + if isinstance(flow_field_grid, th.Tensor): + flow_field_grid = flow_field_grid.cpu().numpy() + + out_shape = list(flow_field_grid.shape) # tuples are immutable + out_shape[-1] = 3 # change 2 to 3 in last axis + out = al.image_from_numpy(np.empty(out_shape), (), (), device=image.device, dtype=image.dtype) + + if out.ndim == 2: + im1, im2, _, _ = match_image_size(image.image, flow_field_grid) + out.image = al.transformation.utils.F.grid_sample(im1, im2) + else: + for i in range(out_shape[-1]): # iterate over last axis + flow_field_grid, image_slice, _, _ = match_image_size(np.squeeze(flow_field_grid), np.squeeze(image.image[...,i].cpu().numpy())) + # https://numpy.org/doc/stable/reference/generated/numpy.expand_dims.html + # https://pytorch.org/docs/stable/nn.functional.html#grid-sample + out.image[..., i] = al.transformation.utils.F.grid_sample(th.from_numpy(image_slice[np.newaxis][np.newaxis]), th.from_numpy(flow_field_grid[np.newaxis])) + # for i in range(3): + # + # im=al.utils.image.create_tensor_image_from_itk_image(im, dtype=dtype, device=('cuda:{}'.format(gpu_device) if gpu_device>=0 else 'cpu')) + # channels.append(al.transformation.utils.warp_image(im, displacement).numpy()) + img=np.uint8(out.image.detach().cpu().numpy()) + # print(img) + # print(img.max()) + # print(img.shape) + return np.squeeze(img) # img[0][0] # np.stack(channels).transpose((1,2,0)) # Copyright 2018 University of Basel, Center for medical Image Analysis and Navigation # @@ -83,498 +123,736 @@ def displace_image(img, displacement, gpu_device, dtype=th.float32): # See the License for the specific language governing permissions and # limitations under the License. - -def affine_register(im1, im2, iterations=1000, lr=0.01, transform_type='similarity', gpu_device=0, opt_cm=True, sigma=[[11,11],[11,11],[3,3]], order=2, pyramid=[[4,4],[2,2]], loss_fn='mse', use_mask=False, interpolation='bicubic'): - assert use_mask==False, "Masking not implemented" - assert transform_type in ['similarity', 'affine', 'rigid', 'non_parametric','bspline','wendland'] - start = time.perf_counter() - - # set the used data type - dtype = th.float32 - # set the device for the computaion to CPU - device = th.device("cuda:{}".format(gpu_device) if gpu_device >=0 else 'cpu') - - # In order to use a GPU uncomment the following line. The number is the device index of the used GPU - # Here, the GPU with the index 0 is used. - # device = th.device("cuda:0") - - # load the image data and normalize to [0, 1] - # add mask to loss function - fixed_image = al.utils.image.create_tensor_image_from_itk_image(sitk.GetImageFromArray(cv2.cvtColor(im1, cv2.COLOR_BGR2GRAY)), dtype=th.float32, device=device)#al.read_image_as_tensor("./practice_reg/1.png", dtype=dtype, device=device)#th.tensor(img1,device='cuda',dtype=dtype)# - moving_image = al.utils.image.create_tensor_image_from_itk_image(sitk.GetImageFromArray(cv2.cvtColor(im2, cv2.COLOR_BGR2GRAY)), dtype=th.float32, device=device)#al.read_image_as_tensor("./practice_reg/2.png", dtype=dtype, device=device)#th.tensor(img2,device='cuda',dtype=dtype)# - - fixed_image, moving_image = al.utils.normalize_images(fixed_image, moving_image) - - # convert intensities so that the object intensities are 1 and the background 0. This is important in order to - # calculate the center of mass of the object - fixed_image.image = 1 - fixed_image.image - moving_image.image = 1 - moving_image.image - - # create pairwise registration object - registration = al.PairwiseRegistration() - - transforms=dict(similarity=al.transformation.pairwise.SimilarityTransformation, - affine=al.transformation.pairwise.AffineTransformation, - rigid=al.transformation.pairwise.RigidTransformation, - non_parametric=al.transformation.pairwise.NonParametricTransformation, - wendland=al.transformation.pairwise.WendlandKernelTransformation, - bspline=al.transformation.pairwise.BsplineTransformation) - constant_flow=None - - if transform_type in ['similarity', 'affine', 'rigid']: - transform_opts=dict(opt_cm=opt_cm) - transform_args=[moving_image] - sigma,fixed_image_pyramid,moving_image_pyramid=[[]],[[]],[[]] - else: - transform_opts=dict(diffeomorphic=opt_cm, device=('cuda:{}'.format(gpu_device) if gpu_device>=0 else 'cpu')) - transform_args=[moving_image.size] - if transform_type in ['bspline','wendland']: - transform_opts['sigma']=sigma - fixed_image_pyramid = al.create_image_pyramid(fixed_image, pyramid) - moving_image_pyramid = al.create_image_pyramid(moving_image, pyramid) - else: - sigma,fixed_image_pyramid,moving_image_pyramid=[[]],[[fixed_image]],[[moving_image]] - if transform_type=='bspline': - transform_opts['order']=order - if transform_type=='wendland': - transform_opts['cp_scale']=order - - for level, (mov_im_level, fix_im_level) in enumerate(zip(moving_image_pyramid, fixed_image_pyramid)): - - # choose the affine transformation model - if transform_type == 'non_parametric': - transform_args[0]=mov_im_level[level].size - elif transform_type in ['bspline','wendland']: - # for bspline, sigma must be positive tuple of ints - # for bspline, smaller sigma tuple means less loss of - # microarchitectural details - - # transform_opts['sigma'] = sigma[level] - transform_opts['sigma'] = (1, 1) - - transformation = transforms[transform_type](*transform_args,**transform_opts) - - # if level > 0 and transform_type=='bspline': - # constant_flow = al.transformation.utils.upsample_displacement(constant_flow, - # mov_im_level.size, - # interpolation=interpolation) - # transformation.set_constant_flow(constant_flow) - - if transform_type in ['similarity', 'affine', 'rigid']: - # initialize the translation with the center of mass of the fixed image - transformation.init_translation(fixed_image) - - registration.set_transformation(transformation) - - loss_fns=dict(mse=al.loss.pairwise.MSE, - ncc=al.loss.pairwise.NCC, - lcc=al.loss.pairwise.LCC, - mi=al.loss.pairwise.MI, - mgf=al.loss.pairwise.NGF, - ssim=al.loss.pairwise.SSIM) - - # choose the Mean Squared Error as image loss - image_loss = loss_fns[loss_fn](fixed_image, moving_image) - - registration.set_image_loss([image_loss]) - - # choose the Adam optimizer to minimize the objective - optimizer = th.optim.Adam(transformation.parameters(), lr=lr, amsgrad=True) - - registration.set_optimizer(optimizer) - registration.set_number_of_iterations(iterations) - - # start the registration - registration.start() - - # if transform_type == 'bspline': - # constant_flow = transformation.get_flow() - - # set the intensities back to the original for the visualisation - fixed_image.image = 1 - fixed_image.image - moving_image.image = 1 - moving_image.image - - # warp the moving image with the final transformation result - displacement = transformation.get_displacement() - warped_image = al.transformation.utils.warp_image(moving_image, displacement) - - end = time.perf_counter() - - print("=================================================================") - - print("Registration done in:", end - start, "s") - if transform_type in ['similarity', 'affine', 'rigid']: - print("Result parameters:") - transformation.print() - - # plot the results - plt.subplot(131) - plt.imshow(fixed_image.numpy(), cmap='gray') - plt.title('Fixed Image') - - plt.subplot(132) - plt.imshow(moving_image.numpy(), cmap='gray') - plt.title('Moving Image') - - plt.subplot(133) - plt.imshow(warped_image.numpy(), cmap='gray') - plt.title('Warped Moving Image') - - if transform_type in ['similarity', 'affine', 'rigid']: - transformation_param = transformation._phi_z - elif transform_type == 'non_parametric': - transformation_param = transformation.trans_parameters - elif transform_type == 'bspline' or transform_type == 'wendland': - transformation_param = transformation._kernel - else: - pass - return displacement, warped_image, transformation_param, registration.loss.data.item() +def img2h(img): + hed=rgb2hed(img) + return rescale_intensity(hed[:, :, 0], out_range=(0, 1)) + +def affine_register(im1, im2, + iterations=1000, + lr=0.01, + transform_type='similarity', + gpu_device=0, + opt_cm=True, + sigma=[[11,11],[11,11],[3,3]], + order=2, + pyramid=[[4,4],[2,2]], + loss_fn='mse', + use_mask=False, + interpolation='bicubic', + half=False, + regularisation_weight=[1,5,50], + moving_image=None, + register_joint_domain=False, + use_hematoxylin=False, + mask1=None, + mask2=None): +# assert use_mask==False, "Masking not implemented" + if use_mask: assert transform_type in ['similarity', 'affine', 'rigid'] + assert transform_type in ['similarity', 'affine', 'rigid', 'non_parametric','bspline','wendland'] + if half: + raise NotImplementedError + start = time.perf_counter() + + # set the used data type + dtype = th.float32# if not half else th.half + # set the device for the computaion to CPU + device = th.device("cuda:{}".format(gpu_device) if gpu_device >=0 else 'cpu') + + # In order to use a GPU uncomment the following line. The number is the device index of the used GPU + # Here, the GPU with the index 0 is used. + # device = th.device("cuda:0") + + # load the image data and normalize to [0, 1] + # add mask to loss function + + mask1,mask2=None,None + + if use_mask and mask1 is None and mask2 is None: + from pathflowai.utils import generate_tissue_mask + mask1,mask2=th.tensor(generate_tissue_mask(im1).astype(float)),th.tensor(generate_tissue_mask(im2).astype(float)) + + fixed_image = al.utils.image.create_tensor_image_from_itk_image(sitk.GetImageFromArray(cv2.cvtColor(im1, cv2.COLOR_BGR2GRAY) if not use_hematoxylin else img2h(im1)), dtype=th.float32, device='cpu')#device,al.read_image_as_tensor("./practice_reg/1.png", dtype=dtype, device=device)#th.tensor(img1,device='cuda',dtype=dtype)# + cm_displacement = None + if not moving_image: + if isinstance(im2, np.ndarray): + moving_image = al.utils.image.create_tensor_image_from_itk_image(sitk.GetImageFromArray(cv2.cvtColor(im2, cv2.COLOR_BGR2GRAY) if not use_hematoxylin else img2h(im2)), dtype=th.float32, device='cpu') + elif isinstance(im2, al.utils.image.Image): + moving_image = im2 + fixed_image, moving_image = al.utils.normalize_images(fixed_image, moving_image) + + # convert intensities so that the object intensities are 1 and the background 0. This is important in order to + # calculate the center of mass of the object + fixed_image.image = 1. - fixed_image.image + moving_image.image = 1. - moving_image.image + + if register_joint_domain: + # TODO: check this control flow + if transform_type == 'bspline': + joint_domain_interpolation = 3 + elif interpolation == 'linear': + joint_domain_interpolation = 2 + else: + joint_domain_interpolation = 1 + fixed_image, _, moving_image, _, cm_displacement = al.get_joint_domain_images(fixed_image, moving_image, default_value=0, interpolator=joint_domain_interpolation, cm_alignment=True, compute_masks=False) + if cm_displacement is not None: + # the domains are not equal and the images were resampled + # https://github.com/airlab-unibas/airlab/blob/80c9d487c012892c395d63c6d937a67303c321d1/airlab/utils/domain.py#L124-L133 + if isinstance(cm_displacement, np.ndarray): + # cm_displacement = th.as_tensor(cm_displacement, dtype=im2.dtype, device=device) + cm_displacement = th.as_tensor(cm_displacement, dtype=th.float, device=device) + im2 = displace_image(im2, cm_displacement, gpu_device) + + if use_mask: + mask1,mask2=Image(mask1, fixed_image.size, fixed_image.spacing, fixed_image.origin),Image(mask2, moving_image.size, moving_image.spacing, moving_image.origin) + + # create pairwise registration object + registration = (al.PairwiseRegistration if transform_type != 'non_parametric' else al.DemonsRegistraion)()#half=half + + transforms=dict(similarity=al.transformation.pairwise.SimilarityTransformation, + affine=al.transformation.pairwise.AffineTransformation, + rigid=al.transformation.pairwise.RigidTransformation, + non_parametric=al.transformation.pairwise.NonParametricTransformation, + wendland=al.transformation.pairwise.WendlandKernelTransformation, + bspline=al.transformation.pairwise.BsplineTransformation) + constant_flow=None + + if transform_type in ['similarity', 'affine', 'rigid']: + transform_opts=dict(opt_cm=opt_cm) + transform_args=[0] + sigma,fixed_image_pyramid,moving_image_pyramid=[[]],[fixed_image],[moving_image] + else: + transform_opts=dict(diffeomorphic=opt_cm, device=('cuda:{}'.format(gpu_device) if gpu_device>=0 else 'cpu')) + transform_args=[moving_image.size] + if transform_type in ['bspline','wendland']: + transform_opts['sigma']=sigma + fixed_image_pyramid = al.create_image_pyramid(fixed_image, pyramid) + moving_image_pyramid = al.create_image_pyramid(moving_image, pyramid) + else: + fixed_image_pyramid,moving_image_pyramid=[fixed_image],[moving_image] + if transform_type=='bspline': + transform_opts['order']=order + if transform_type=='wendland': + transform_opts['cp_scale']=order + + # transform_opts['half']=half + + for level, (mov_im_level, fix_im_level) in enumerate(zip(moving_image_pyramid, fixed_image_pyramid)): + mov_im_level, fix_im_level, _, _ = match_image_size(mov_im_level, fix_im_level) + + mov_im_level = mov_im_level.to(dtype=th.float32, device=device) + fix_im_level = fix_im_level.to(dtype=th.float32, device=device) + + # choose the affine transformation model + if transform_type == 'non_parametric': + transform_args[0]=mov_im_level.size + elif transform_type in ['bspline','wendland']: + # for bspline, sigma must be positive tuple of ints + # for bspline, smaller sigma tuple means less loss of + # microarchitectural details + transform_args[0]=mov_im_level.size + transform_opts['sigma'] = sigma[level] + else: + transform_args[0]=mov_im_level + + transformation = transforms[transform_type](*transform_args,**transform_opts) + + # if half: + # mov_im_level=mov_im_level.to(dtype=th.float16, device=device) + + transformation=transformation.to(dtype=th.float32, device=device) # dtype=th.float32, if not half else th.float16 + + if level > 0 and transform_type in ['bspline','wendland']: + print(interpolation) + constant_flow = al.transformation.utils.upsample_displacement(constant_flow, + mov_im_level.size, + interpolation=interpolation) + constant_flow, fix_im_level, _, _ = match_image_size(constant_flow, fix_im_level) + transformation.set_constant_flow(constant_flow.to(device=device)) + + # + if transform_type in ['similarity', 'affine', 'rigid']: + # initialize the translation with the center of mass of the fixed image + transformation.init_translation(fix_im_level) + # if half: + # fix_im_level=fix_im_level.to(dtype=th.float16, device=device) + # transformation._dtype=th.float16 + # transformation._device=device + + + if isinstance(lr, float): + optim_lr = lr + elif isinstance(lr, (list, tuple)): + optim_lr = lr[level] + else: + pass + optimizer = th.optim.Adam(transformation.parameters(), lr=optim_lr, amsgrad=True) + # opt_level = "O2" if half else "O1" + # transformation, optimizer = amp.initialize(transformation, optimizer, opt_level=opt_level) + + registration.set_transformation(transformation) + + loss_fns=dict(mse=al.loss.pairwise.MSE, + ncc=al.loss.pairwise.NCC, + lcc=al.loss.pairwise.LCC, + mi=al.loss.pairwise.MI, + mgf=al.loss.pairwise.NGF, + ssim=al.loss.pairwise.SSIM) + + # the loss function requires that both images are on the same device + # https://github.com/airlab-unibas/airlab/blob/80c9d487c012892c395d63c6d937a67303c321d1/airlab/loss/pairwise.py#L47 + mov_im_level = mov_im_level.to(dtype=th.float32, device=device) + fix_im_level = fix_im_level.to(dtype=th.float32, device=device) + + if use_mask: + mask1 = mask1.to(dtype=th.float32, device=device) + mask2 = mask2.to(dtype=th.float32, device=device) + + image_loss = loss_fns[loss_fn](fix_im_level, mov_im_level, fixed_mask=mask1, moving_mask=mask2) + + registration.set_image_loss([image_loss]) + + if transform_type in ['non_parametric','wendland','bspline']: + if transform_type=='non_parametric': + regulariser = al.regulariser.demons.GaussianRegulariser(mov_im_level.spacing,sigma=sigma[level], dtype=dtype, device=device) + registration.set_regulariser([regulariser]) + else: + regulariser = al.regulariser.displacement.DiffusionRegulariser(mov_im_level.spacing) + regulariser.set_weight(regularisation_weight[level]) + registration.set_regulariser_displacement([regulariser]) + + # choose the Adam optimizer to minimize the objective + + registration.set_optimizer(optimizer) + registration.set_number_of_iterations(iterations) + + # start the registration + registration.start() + + if transform_type in ['bspline','wendland']: + constant_flow = transformation.get_flow() + constant_flow, fix_im_level, _, _ = match_image_size(constant_flow, fix_im_level) + + # set the intensities back to the original for the visualisation + fixed_image.image = 1. - fixed_image.image + moving_image.image = 1. - moving_image.image + + # warp the moving image with the final transformation result + displacement = transformation.get_displacement() + warped_image = al.transformation.utils.warp_image(moving_image.to(device=displacement.device), displacement) + + end = time.perf_counter() + + print("=================================================================") + + print("Registration done in:", end - start, "s") + if transform_type in ['similarity', 'affine', 'rigid']: + print("Result parameters:") + transformation.print() + + # # plot the results + # plt.subplot(131) + # plt.imshow(fixed_image.numpy(), cmap='gray') + # plt.title('Fixed Image') + # + # plt.subplot(132) + # plt.imshow(moving_image.numpy(), cmap='gray') + # plt.title('Moving Image') + # + # plt.subplot(133) + # plt.imshow(warped_image.numpy(), cmap='gray') + # plt.title('Warped Moving Image') + + if transform_type in ['similarity', 'affine', 'rigid']: + transformation_param = list(transformation.named_parameters())#._phi_z + elif transform_type == 'non_parametric': + transformation_param = transformation.trans_parameters + elif transform_type == 'bspline' or transform_type == 'wendland': + transformation_param = transformation._kernel + else: + pass + return { + 'displacement': displacement, + 'warped_image': warped_image, + 'moving_image': moving_image, + 'transformation_param': transformation_param, + 'loss': registration.loss, + 'cm_displacement': cm_displacement + } def get_loss(im1,im2,gpu_device): - dh=int(np.abs((im1.shape[0]-im2.shape[0]))) - if im1.shape[0]>im2.shape[0]: - img2 = cv2.copyMakeBorder(im2, dh, 0, 0, 0, cv2.BORDER_CONSTANT, value=[255,255,255]) - elif im1.shape[0]im2.shape[1]: - im2 = cv2.copyMakeBorder(im2, 0, 0, 0, dw, cv2.BORDER_CONSTANT, value=[255,255,255]) - elif im1.shape[1]=0: - th.cuda.empty_cache() - return loss + dh=int(np.abs((im1.shape[0]-im2.shape[0]))) + if im1.shape[0]>im2.shape[0]: + img2 = cv2.copyMakeBorder(im2, dh, 0, 0, 0, cv2.BORDER_CONSTANT, value=[255,255,255]) + elif im1.shape[0]im2.shape[1]: + im2 = cv2.copyMakeBorder(im2, 0, 0, 0, dw, cv2.BORDER_CONSTANT, value=[255,255,255]) + elif im1.shape[1]=0: + th.cuda.empty_cache() + return affine_output['loss'] def rotate_detector(im1,im2,gpu_device): - angles={} - for k,angle in enumerate([0, 90, 180, 270]): - im_test=np.rot90(im2,k)#rotate(im2,angle, resize=True, center=None, order=1, mode='constant', cval=255., clip=False, preserve_range=False) - angles[angle]=get_loss(im1,im_test,gpu_device) - return angles + angles={} + for k,angle in enumerate([0, 90, 180, 270]): + im_test=np.rot90(im2,k)#rotate(im2,angle, resize=True, center=None, order=1, mode='constant', cval=255., clip=False, preserve_range=False) + angles[angle]=get_loss(im1,im_test,gpu_device) + return angles def correct_rotation(im1, im2, scaling_factor=4, gpu_device=0): - print("Resizing image 1 for rotation check.") - # im1_small=resize(im1, (im1.shape[0] // scaling_factor, im1.shape[1] // scaling_factor), - # anti_aliasing=True)*255. - im1_small=cv2.resize(im1,(int(im1.shape[1] // scaling_factor), int(im1.shape[0] // scaling_factor))) - print("Resizing image 2 for rotation check.") - # im2_small=resize(im2, (im2.shape[0] // scaling_factor, im2.shape[1] // scaling_factor), - # anti_aliasing=True)*255. - im2_small=cv2.resize(im2,(int(im2.shape[1] // scaling_factor), int(im2.shape[0] // scaling_factor))) - - print("Detecting ideal rotation.") - with suppress_stdout(): - rotation_loss_dict=rotate_detector(im1_small,im2_small,gpu_device) - - if gpu_device>=0: - th.cuda.empty_cache() - - rotations=np.array(list(rotation_loss_dict.items())) - angle=int(rotations[np.argmin(rotations[:,1]),0]) - print("Ideal rotation angle: {}.".format(angle)) - return np.rot90(im2,angle//90) - -def rotate_image(mat, angle): - """ - https://stackoverflow.com/questions/43892506/opencv-python-rotate-image-without-cropping-sides/47248339 - Rotates an image (angle in degrees) and expands image to avoid cropping - """ - - height, width = mat.shape[:2] # image shape has 3 dimensions - image_center = (width/2, height/2) # getRotationMatrix2D needs coordinates in reverse order (width, height) compared to shape - - rotation_mat = cv2.getRotationMatrix2D(image_center, angle, 1.) - - # rotation calculates the cos and sin, taking absolutes of those. - abs_cos = np.abs(rotation_mat[0,0]) - abs_sin = np.abs(rotation_mat[0,1]) - - # find the new width and height bounds - bound_w = int(height * abs_sin + width * abs_cos) - bound_h = int(height * abs_cos + width * abs_sin) - - # subtract old image center (bringing image back to origo) and adding the new image center coordinates - rotation_mat[0, 2] += bound_w/2 - image_center[0] - rotation_mat[1, 2] += bound_h/2 - image_center[1] - - # rotate image with the new bounds and translated rotation matrix - rotated_mat = cv2.warpAffine(mat, rotation_mat, (bound_w, bound_h), borderValue=(255,255,255)) - return rotated_mat + print("Resizing image 1 for rotation check.") + # im1_small=resize(im1, (im1.shape[0] // scaling_factor, im1.shape[1] // scaling_factor), + # anti_aliasing=True)*255. + im1_small=cv2.resize(im1,(int(im1.shape[1] // scaling_factor), int(im1.shape[0] // scaling_factor))) + print("Resizing image 2 for rotation check.") + # im2_small=resize(im2, (im2.shape[0] // scaling_factor, im2.shape[1] // scaling_factor), + # anti_aliasing=True)*255. + im2_small=cv2.resize(im2,(int(im2.shape[1] // scaling_factor), int(im2.shape[0] // scaling_factor))) + + print("Detecting ideal rotation.") + with suppress_stdout(): + rotation_loss_dict=rotate_detector(im1_small,im2_small,gpu_device) + + if gpu_device>=0: + th.cuda.empty_cache() + + rotations=np.array(list(rotation_loss_dict.items())) + angle=int(rotations[np.argmin(rotations[:,1]),0]) + print("Ideal rotation angle: {}.".format(angle)) + return np.rot90(im2,angle//90) + +def rotate_image(mat, angle, center=None): + """ + https://stackoverflow.com/questions/43892506/opencv-python-rotate-image-without-cropping-sides/47248339 + Rotates an image (angle in degrees) and expands image to avoid cropping + """ + + height, width = mat.shape[:2] # image shape has 3 dimensions + image_center = (width/2, height/2) # getRotationMatrix2D needs coordinates in reverse order (width, height) compared to shape + + rotation_mat = cv2.getRotationMatrix2D(image_center if center is None else center, angle, 1.) + + # rotation calculates the cos and sin, taking absolutes of those. + abs_cos = np.abs(rotation_mat[0,0]) + abs_sin = np.abs(rotation_mat[0,1]) + + # find the new width and height bounds + bound_w = int(height * abs_sin + width * abs_cos) + bound_h = int(height * abs_cos + width * abs_sin) + + # subtract old image center (bringing image back to origo) and adding the new image center coordinates + rotation_mat[0, 2] += bound_w/2 - image_center[0] + rotation_mat[1, 2] += bound_h/2 - image_center[1] + + # rotate image with the new bounds and translated rotation matrix + rotated_mat = cv2.warpAffine(mat, rotation_mat, (bound_w, bound_h), borderValue=(255,255,255)) + return rotated_mat def match_image_size(img1,img2,black_background=False): - white=int(black_background==False) - fill_color=(np.array([255,255,255])*white).astype(int).tolist() - dh=int(np.abs((img1.shape[0]-img2.shape[0]))) - if img1.shape[0]>img2.shape[0]: - img2 = cv2.copyMakeBorder(img2, dh//2+dh%2, dh//2, 0, 0, cv2.BORDER_CONSTANT, value=fill_color) - elif img1.shape[0]img2.shape[1]: - img2 = cv2.copyMakeBorder(img2, 0, 0, dw//2+dw%2, dw//2, cv2.BORDER_CONSTANT, value=fill_color) - elif img1.shape[1]img2.shape[0]: + img2 = cv2.copyMakeBorder(img2, dh//2+dh%2, dh//2, 0, 0, cv2.BORDER_CONSTANT, value=fill_color) + elif img1.shape[0]img2.shape[1]: + img2 = cv2.copyMakeBorder(img2, 0, 0, dw//2+dw%2, dw//2, cv2.BORDER_CONSTANT, value=fill_color) + elif img1.shape[1]=0: - th.cuda.empty_cache() - - print("Locating Sections.") - - mask,labels=label_objects(im1, connectivity=connectivity, min_object_size=min_object_size, apply_watershed=apply_watershed) - mask2,labels2=label_objects(im2, connectivity=connectivity, min_object_size=min_object_size, apply_watershed=apply_watershed) - - print("Estimating section properties.") - - props = pd.DataFrame(measure.regionprops_table(labels, cv2.cvtColor(im1, cv2.COLOR_BGR2GRAY),properties=properties)) - props2 = pd.DataFrame(measure.regionprops_table(labels2, cv2.cvtColor(im2, cv2.COLOR_BGR2GRAY),properties=properties)) - - bboxes=pd.DataFrame(measure.regionprops_table(labels,cv2.cvtColor(im1, cv2.COLOR_BGR2GRAY),properties=['bbox'])) - bboxes2=pd.DataFrame(measure.regionprops_table(labels2,cv2.cvtColor(im2, cv2.COLOR_BGR2GRAY),properties=['bbox'])) - - print("Matching tissue sections.") - - sections=list(get_matched_tissue(props,props2)) - - N=len(sections) - - if file_ext=='.npy': - - del im1, im2 - - im1=np.load(im1_fname,mmap_mode='r+') - im2=np.load(im2_fname,mmap_mode='r+') - - print("Performing alignments on pairs of sections.") - - for idx,(i,j) in enumerate(sections): - if os.path.exists(img_out1)==False or os.path.exists(img_out1)==False or overwrite: - - try: - print("[{}/{}] - Extracting and rotating sections {} and {}".format(idx+1,N,i,j)) - angle=-props.loc[i,'orientation']*180./(np.pi) - img=im1[bboxes.iloc[i,0]:bboxes.iloc[i,2],bboxes.iloc[i,1]:bboxes.iloc[i,3]].copy()#.copy() - img[labels[bboxes.iloc[i,0]:bboxes.iloc[i,2],bboxes.iloc[i,1]:bboxes.iloc[i,3]]!=(i+1)]=255. - img1=rotate_image(img, angle)#np.uint8(rotate(img, -props.loc[i,'orientation']*180./(np.pi), resize=True, center=None, order=1, mode='constant', cval=1., clip=clip, preserve_range=False)*255.) - - angle=-props2.loc[j,'orientation']*180./(np.pi) - img=im2[bboxes2.iloc[j,0]:bboxes2.iloc[j,2],bboxes2.iloc[j,1]:bboxes2.iloc[j,3]].copy()#.copy() - img[labels2[bboxes2.iloc[j,0]:bboxes2.iloc[j,2],bboxes2.iloc[j,1]:bboxes2.iloc[j,3]]!=(j+1)]=255. - img2=rotate_image(img, angle)#np.uint8(rotate(img, -props2.loc[j,'orientation']*180./(np.pi), resize=True, center=None, order=1, mode='constant', cval=1., clip=clip, preserve_range=False)*255.) - - print("[{}/{}] - Trimming sections to proper W & H".format(idx+1,N)) - - c=int(img1.shape[1]/2.) - w=max(int(props.loc[i,'minor_axis_length']/2),int(props2.loc[j,'minor_axis_length']/2*mult_factor)) - h=max(int(img1.shape[0]/2.),int(img2.shape[1]/2.)) - img1=img1[:,c-w:c+w] - - c=int(img2.shape[1]/2.) - img2=img2[:,c-w:c+w] - - img1,img2=match_image_size(img1,img2,black_background=black_background) - - print("[{}/{}] - Begin alignment of sections.".format(idx+1,N)) - - with (suppress_stdout() if not verbose else contextlib.suppress()): - new_img=displace_image(img2,affine_register(img1, img2, gpu_device=gpu_device, lr=lr, loss_fn=loss_fn, transform_type=transform_type, iterations=iterations, opt_cm=opt_cm, sigma=sigma, order=order, pyramid=pyramid,interpolation=interpolation)[0],gpu_device=gpu_device) # new tri, output he as well - - if gpu_device>=0: - th.cuda.empty_cache() - - print("[{}/{}] - Writing registered sections to file.".format(idx+1,N)) - - cv2.imwrite(img_out1,cv2.cvtColor(img1,cv2.COLOR_BGR2RGB)) - cv2.imwrite(img_out2,cv2.cvtColor(new_img,cv2.COLOR_BGR2RGB)) - - except Exception as e: - print(str(e)) - - else: - - im1,im2=match_image_size(im1,im2,black_background=black_background) + im2_fname='B.npy', + connectivity=8, + apply_watershed=False, + overwrite=False, + min_object_size=50000, + clip=True, + fix_rotation=True, + mult_factor=2.0, + scaling_factor=4., + output_dir='practice_reg_results/images', + properties=['area', + 'convex_area', + 'eccentricity', + 'major_axis_length', + 'minor_axis_length', + 'inertia_tensor', + 'inertia_tensor_eigvals', + 'perimeter', + 'orientation'], + gpu_device=0, + max_rotation_vertical_px=0, + loss_fn='mse', + lr=[0.01]*3, + transform_type='rigid', + iterations=1000, + no_segment_analysis=False, + black_background=False, + verbose=False, + opt_cm=True, + sigma=[[11,11],[11,11],[3,3]], + order=2, + pyramid=[[4,4],[2,2]], + interpolation='bicubic', + half=False, + regularisation_weight=[1,5,50], + points1='', + points2='', + pre_transform=''): + + print("Loading images.") + + _, file_ext = os.path.splitext(im1_fname) + + # TODO: don't rerun splitext, see variable file_ext + os.makedirs(output_dir, exist_ok=True) + # img_splitext1 and img_splitext2 are pairs, (root, ext) + img_splitext1 = os.path.splitext(os.path.basename(im1_fname)) + img_out1 = os.path.join( + output_dir, + f"{img_splitext1[0]}_registered{img_splitext1[1]}" + ) + img_splitext2 = os.path.splitext(os.path.basename(im2_fname)) + img_out2 = os.path.join( + output_dir, + f"{img_splitext2[0]}_registered{img_splitext2[1]}" + ) + + # TODO: load images with airlab.utils.image.read_image_as_tensor + # https://airlab.readthedocs.io/en/latest/_modules/airlab/utils/image.html#read_image_as_tensor + + if file_ext=='.npy': + im1=np.load(im1_fname) + im2=np.load(im2_fname) + else: + im1=cv2.imread(im1_fname) + im2=cv2.imread(im2_fname) + + tre=-1 + if not no_segment_analysis: + + if max_rotation_vertical_px: + scaling_factor=((sum(im1.shape[:2])+sum(im2.shape[:2]))/4.)/float(max_rotation_vertical_px) + print("New scaling factor: {}".format(scaling_factor)) + + if fix_rotation: + print("Adjusting 90 degree rotations.") + + im2=correct_rotation(im1, im2, scaling_factor=scaling_factor, gpu_device=gpu_device) + + if gpu_device>=0: + th.cuda.empty_cache() + + print("Locating Sections.") + + mask,labels=label_objects(im1, connectivity=connectivity, min_object_size=min_object_size, apply_watershed=apply_watershed) + mask2,labels2=label_objects(im2, connectivity=connectivity, min_object_size=min_object_size, apply_watershed=apply_watershed) + + print("Estimating section properties.") + + props = pd.DataFrame(measure.regionprops_table(labels, cv2.cvtColor(im1, cv2.COLOR_BGR2GRAY),properties=properties)) + props2 = pd.DataFrame(measure.regionprops_table(labels2, cv2.cvtColor(im2, cv2.COLOR_BGR2GRAY),properties=properties)) + + bboxes=pd.DataFrame(measure.regionprops_table(labels,cv2.cvtColor(im1, cv2.COLOR_BGR2GRAY),properties=['bbox'])) + bboxes2=pd.DataFrame(measure.regionprops_table(labels2,cv2.cvtColor(im2, cv2.COLOR_BGR2GRAY),properties=['bbox'])) + + print("Matching tissue sections.") + + sections=list(get_matched_tissue(props,props2)) + + N=len(sections) + + if file_ext=='.npy': + + del im1, im2 + + im1=np.load(im1_fname,mmap_mode='r+') + im2=np.load(im2_fname,mmap_mode='r+') + + print("Performing alignments on pairs of sections.") + + for idx,(i,j) in enumerate(sections): + if os.path.exists(img_out1)==False or os.path.exists(img_out1)==False or overwrite: + + try: + print("[{}/{}] - Extracting and rotating sections {} and {}".format(idx+1,N,i,j)) + angle=-props.loc[i,'orientation']*180./(np.pi) + img=im1[bboxes.iloc[i,0]:bboxes.iloc[i,2],bboxes.iloc[i,1]:bboxes.iloc[i,3]].copy()#.copy() + img[labels[bboxes.iloc[i,0]:bboxes.iloc[i,2],bboxes.iloc[i,1]:bboxes.iloc[i,3]]!=(i+1)]=255. + img1=rotate_image(img, angle)#np.uint8(rotate(img, -props.loc[i,'orientation']*180./(np.pi), resize=True, center=None, order=1, mode='constant', cval=1., clip=clip, preserve_range=False)*255.) + + angle=-props2.loc[j,'orientation']*180./(np.pi) + img=im2[bboxes2.iloc[j,0]:bboxes2.iloc[j,2],bboxes2.iloc[j,1]:bboxes2.iloc[j,3]].copy()#.copy() + img[labels2[bboxes2.iloc[j,0]:bboxes2.iloc[j,2],bboxes2.iloc[j,1]:bboxes2.iloc[j,3]]!=(j+1)]=255. + img2=rotate_image(img, angle)#np.uint8(rotate(img, -props2.loc[j,'orientation']*180./(np.pi), resize=True, center=None, order=1, mode='constant', cval=1., clip=clip, preserve_range=False)*255.) + + print("[{}/{}] - Trimming sections to proper W & H".format(idx+1,N)) + + c=int(img1.shape[1]/2.) + w=max(int(props.loc[i,'minor_axis_length']/2),int(props2.loc[j,'minor_axis_length']/2*mult_factor)) + h=max(int(img1.shape[0]/2.),int(img2.shape[1]/2.)) + img1=img1[:,c-w:c+w] + + c=int(img2.shape[1]/2.) + img2=img2[:,c-w:c+w] + + img1,img2=match_image_size(img1,img2,black_background=black_background)[:2] + + print("[{}/{}] - Begin alignment of sections.".format(idx+1,N)) + + with (suppress_stdout() if not verbose else contextlib.suppress()): + new_img=displace_image(img2,affine_register(img1, img2, gpu_device=gpu_device, lr=lr, loss_fn=loss_fn, transform_type=transform_type, iterations=iterations, opt_cm=opt_cm, sigma=sigma, order=order, pyramid=pyramid,interpolation=interpolation, half=half, regularisation_weight=regularisation_weight)['displacement'],gpu_device=gpu_device, dtype=th.float32 if not half else th.half) # new tri, output he as well + + if gpu_device>=0: + th.cuda.empty_cache() + + print("[{}/{}] - Writing registered sections to file.".format(idx+1,N)) + + cv2.imwrite(img_out1,cv2.cvtColor(img1,cv2.COLOR_BGR2RGB)) + cv2.imwrite(img_out2,cv2.cvtColor(new_img,cv2.COLOR_BGR2RGB)) + + except Exception as e: + print(str(e)) + + else: - print("Performing registration.") + im1,im2,dw,dh=match_image_size(im1,im2,black_background=black_background) - with (suppress_stdout() if not verbose else contextlib.suppress()): - new_img=displace_image(im2,affine_register(im1, im2, gpu_device=gpu_device, lr=lr, loss_fn=loss_fn, transform_type=transform_type, iterations=iterations, opt_cm=opt_cm, sigma=sigma, order=order, pyramid=pyramid,interpolation=interpolation)[0],gpu_device=gpu_device) # new tri, output he as well + print("Performing registration.") - if gpu_device>=0: - th.cuda.empty_cache() + displacements=[] + with (suppress_stdout() if not verbose else contextlib.suppress()): + if pre_transform: + # tensor_img1 = al.utils.image.create_tensor_image_from_itk_image(sitk.GetImageFromArray(cv2.cvtColor(im1, cv2.COLOR_BGR2GRAY)), dtype=th.float32) + # tensor_img2 = al.utils.image.create_tensor_image_from_itk_image(sitk.GetImageFromArray(cv2.cvtColor(im2, cv2.COLOR_BGR2GRAY)), dtype=th.float32) + # displacement_1,warp_mv_im_1=affine_register(im1, im2, gpu_device=gpu_device, lr=lr, loss_fn=loss_fn, transform_type=pre_transform, iterations=iterations, opt_cm=opt_cm, sigma=sigma, order=order, pyramid=pyramid,interpolation=interpolation, half=half, regularisation_weight=regularisation_weight)[:2] + # fixed_image, moving_image = al.utils.normalize_images(al.utils.image.image_from_numpy(im1), al.utils.image.image_from_numpy(im2)) + # fixed_image, f_mask, moving_image, m_mask, cm_displacement = al.get_joint_domain_images(im1, im2, default_value=0, interpolator=2, cm_alignment=True, compute_masks=False) + # output_reg=affine_register(im1, im2, gpu_device=gpu_device, lr=lr, loss_fn=loss_fn, transform_type=pre_transform, iterations=iterations, opt_cm=opt_cm, sigma=sigma, order=order, pyramid=pyramid,interpolation=interpolation, half=half, regularisation_weight=regularisation_weight,register_joint_domain=True) + # cm_displacement,displacement_1,warp_mv_im_1=output_reg['cm_displacement'],output_reg['displacement'], output_reg['warped_image']) + # new_img = warp_mv_im_1 - print("Writing registered section to file.") + print('Using pre_transform {}'.format(pre_transform)) + reg1 = affine_register(im1, im2, gpu_device=gpu_device, lr=lr, loss_fn=loss_fn, transform_type=pre_transform, iterations=iterations, opt_cm=opt_cm, sigma=sigma, order=order, pyramid=pyramid,interpolation=interpolation, half=half, regularisation_weight=regularisation_weight, register_joint_domain=True) + displacement_1 = reg1['displacement'] + new_img = displace_image(im2, displacement_1, gpu_device) - cv2.imwrite(img_out1, cv2.cvtColor(im1,cv2.COLOR_BGR2RGB)) - cv2.imwrite(img_out2, cv2.cvtColor(new_img,cv2.COLOR_BGR2RGB)) + reg2 = affine_register(im1, im2, moving_image=reg1['moving_image'], gpu_device=gpu_device, lr=lr, loss_fn=loss_fn, transform_type=transform_type, iterations=iterations, opt_cm=opt_cm, sigma=sigma, order=order, pyramid=pyramid, interpolation=interpolation, half=half, regularisation_weight=regularisation_weight, register_joint_domain=True) + displacement_2 = reg2['displacement'] + new_img = displace_image(new_img, displacement_2, gpu_device) + displacements = [ + [displacement_2, (reg2['warped_image'], reg2['moving_image'])], + [displacement_1, (reg1['warped_image'], reg1['moving_image'])], + ] + + # print(cm_displacement) + # print(type(cm_displacement)) + + # displacement_1,warp_mv_im_1=affine_register(im1, im2, gpu_device=gpu_device, lr=lr, loss_fn=loss_fn, transform_type=pre_transform, iterations=iterations, opt_cm=opt_cm, sigma=sigma, order=order, pyramid=pyramid,interpolation=interpolation, half=half, regularisation_weight=regularisation_weight)[:2] + # displacements.append([displacement_1,warp_mv_im_1]) + # displacement_2,warp_mv_im_2=affine_register(im1, warp_mv_im_1[0], gpu_device=gpu_device, lr=lr, loss_fn=loss_fn, transform_type=transform_type, iterations=iterations, opt_cm=opt_cm, sigma=sigma, order=order, pyramid=pyramid,interpolation=interpolation, half=half, regularisation_weight=regularisation_weight, moving_image=warp_mv_im_1[1])[:2] + # displacements.append([displacement_2,warp_mv_im_2]) + # new_img=displace_image(fixed_image.numpy(),cm_displacement,gpu_device=gpu_device, dtype=th.float32 if not half else th.half) # new tri, output he as well + # new_img=displace_image(new_img,displacement_2,gpu_device=gpu_device, dtype=th.float32 if not half else th.half) # new tri, output he as well + + # mask,labels=label_objects(im1, connectivity=connectivity, min_object_size=min_object_size, apply_watershed=apply_watershed) + # mask2,labels2=label_objects(im2, connectivity=connectivity, min_object_size=min_object_size, apply_watershed=apply_watershed) + # + # print("Estimating section properties.") + # + # props = pd.DataFrame(measure.regionprops_table(labels, cv2.cvtColor(im1, cv2.COLOR_BGR2GRAY),properties=['bbox', 'area', 'orientation'])) + # props2 = pd.DataFrame(measure.regionprops_table(labels2, cv2.cvtColor(im2, cv2.COLOR_BGR2GRAY),properties=['bbox', 'area', 'orientation'])) + # print(props) + # print(props2) + # + # + # bboxes=pd.DataFrame(measure.regionprops_table(labels,cv2.cvtColor(im1, cv2.COLOR_BGR2GRAY),properties=['bbox'])) + # bboxes2=pd.DataFrame(measure.regionprops_table(labels2,cv2.cvtColor(im2, cv2.COLOR_BGR2GRAY),properties=['bbox'])) + else: + reg = affine_register(im1, im2, transform_type=transform_type, gpu_device=gpu_device, lr=lr, loss_fn=loss_fn, iterations=iterations, opt_cm=opt_cm, sigma=sigma, order=order, pyramid=pyramid,interpolation=interpolation, half=half, regularisation_weight=regularisation_weight, register_joint_domain=True) + displacement = reg['displacement'] + displacements.append([displacement, (reg['warped_image'], reg['moving_image'])]) + new_img = displace_image(im2, displacement, gpu_device) + + if points1 and points2 and os.path.exists(points1) and os.path.exists(points2): + displacements_final=[] + for displacement,(warp,m_im) in displacements: + displacement = al.transformation.utils.unit_displacement_to_displacement(displacement) # unit measures to image domain measures + displacement = al.create_displacement_image_from_image(displacement, m_im) + displacements_final.append(displacement) + points1=pd.read_csv(points1,index_col=0).values + points2=pd.read_csv(points2,index_col=0).values + points2[:,0]+=dw + points2[:,1]+=dh + for displacement in displacements_final: + points2=al.utils.points.Points.transform(points2,displacement) + tre=al.utils.points.Points.TRE(points1,points2) + else: + tre=-1 + if gpu_device>=0: + th.cuda.empty_cache() + + print("Writing registered section to file.") + + if isinstance(new_img, al.utils.image.Image): + new_img = new_img.numpy() + + cv2.imwrite(img_out1, cv2.cvtColor(im1,cv2.COLOR_BGR2RGB)) + cv2.imwrite(img_out2, cv2.cvtColor(new_img,cv2.COLOR_BGR2RGB)) + return tre class Commands(object): - def __init__(self): - pass - - def apply_drop2_transform(self, - source_image='A.png', - ref_image='B.png', - dx='warp_field_x.nii.gz', - dy='warp_field_y.nii.gz', - gpu_device=-1, - output_file=''):# list of -1s and 1s - import nibabel - assert source_image.split('.')[-1]=='png' and ref_image.split('.')[-1]=='png' - assert os.path.exists(source_image) - assert os.path.exists(ref_image) - flip_pattern=[1,1,1,1] - flip_xy=False - source_img=cv2.imread(source_image)#cv2.cvtColor(,cv2.COLOR_BGR2RGB) - ref_img=cv2.imread(ref_image) - source_img=cv2.resize(source_img,ref_img.shape[-2::-1]) - dx,dy=nibabel.load(dx).get_fdata(),nibabel.load(dy).get_fdata() - displacement=th.tensor(np.concatenate([dx[::flip_pattern[0],::flip_pattern[1]],dy[::flip_pattern[2],::flip_pattern[3]]][::(-1 if flip_xy else 1)],-1)[::1,::1,...].copy()).unsqueeze(0).float().permute(0,2,1,3) - for dim in range(displacement.shape[-1]): - displacement[...,dim]=2.0*displacement[...,dim]/float(displacement.shape[-dim - 2] - 1) - if gpu_device >= 0: - displacement=displacement.cuda() - new_img=displace_image(source_img, displacement, gpu_device) - cv2.imwrite(source_image.replace('.png','_warped.png') if not output_file else output_file,new_img) - - def register_images(self, - im1='A.npy', - im2='B.npy', - connectivity=8, - overwrite=False, - clip=True, - min_object_size=50000, - fix_rotation=True, - mult_factor=1.8, - scaling_factor=4., - max_rotation_vertical_px=0, - output_dir='reg_results/images', - properties=['area', - 'convex_area', - 'eccentricity', - 'major_axis_length', - 'minor_axis_length', - 'inertia_tensor', - 'inertia_tensor_eigvals', - 'perimeter', - 'orientation'], - gpu_device=0, - loss_fn='mse', - lr=0.01, - transform_type='rigid', - iterations=1000, - no_segment_analysis=False, - black_background=False, - verbose=False, - opt_cm=True, - sigma=[[11,11],[11,11],[3,3]], - order=2, - pyramid=[[4,4],[2,2]], - interpolation='bicubic'): - register_images_(im1_fname=im1, - im2_fname=im2, - connectivity=connectivity, - apply_watershed=False, - overwrite=overwrite, - min_object_size=min_object_size, - clip=clip, - fix_rotation=fix_rotation, - mult_factor=mult_factor, - scaling_factor=scaling_factor, - output_dir=output_dir, - properties=properties, - gpu_device=gpu_device, - max_rotation_vertical_px=max_rotation_vertical_px, - loss_fn=loss_fn, - lr=lr, - transform_type=transform_type, - iterations=iterations, - no_segment_analysis=no_segment_analysis, - black_background=black_background, - verbose=verbose, - opt_cm=opt_cm, - sigma=sigma, - order=order, - pyramid=pyramid, - interpolation=interpolation) + def __init__(self): + pass + + def apply_drop2_transform(self, + source_image='A.png', + ref_image='B.png', + dx='warp_field_x.nii.gz', + dy='warp_field_y.nii.gz', + gpu_device=-1, + output_file=''):# list of -1s and 1s + import nibabel + assert source_image.split('.')[-1]=='png' and ref_image.split('.')[-1]=='png' + assert os.path.exists(source_image) + assert os.path.exists(ref_image) + flip_pattern=[1,1,1,1] + flip_xy=False + source_img=cv2.imread(source_image)#cv2.cvtColor(,cv2.COLOR_BGR2RGB) + ref_img=cv2.imread(ref_image) + source_img=cv2.resize(source_img,ref_img.shape[-2::-1]) + dx,dy=nibabel.load(dx).get_fdata(),nibabel.load(dy).get_fdata() + displacement=th.tensor(np.concatenate([dx[::flip_pattern[0],::flip_pattern[1]],dy[::flip_pattern[2],::flip_pattern[3]]][::(-1 if flip_xy else 1)],-1)[::1,::1,...].copy()).unsqueeze(0).float().permute(0,2,1,3) + for dim in range(displacement.shape[-1]): + displacement[...,dim]=2.0*displacement[...,dim]/float(displacement.shape[-dim - 2] - 1) + if gpu_device >= 0: + displacement=displacement.cuda() + new_img=displace_image(source_img, displacement, gpu_device) + cv2.imwrite(source_image.replace('.png','_warped.png') if not output_file else output_file,new_img) + + def compress_image(self, + im='A.png', + compression_factor=2., + im_out='A.compressed.png'): + fx=fy=1/compression_factor + cv2.imwrite(im_out,cv2.resize(cv2.imread(im),None,fx=fx,fy=fy,interpolation=cv2.INTER_CUBIC)) + + def register_images(self, + im1='A.npy', + im2='B.npy', + connectivity=8, + overwrite=False, + clip=True, + min_object_size=50000, + fix_rotation=True, + mult_factor=1.8, + scaling_factor=4., + max_rotation_vertical_px=0, + output_dir='reg_results/images', + properties=['area', + 'convex_area', + 'eccentricity', + 'major_axis_length', + 'minor_axis_length', + 'inertia_tensor', + 'inertia_tensor_eigvals', + 'perimeter', + 'orientation'], + gpu_device=0, + loss_fn='mse', + lr=0.01, + transform_type='rigid', + iterations=1000, + no_segment_analysis=False, + black_background=False, + verbose=False, + opt_cm=True, + sigma=[[11,11],[11,11],[3,3]], + order=2, + pyramid=[[4,4],[2,2]], + interpolation='bicubic', + half=False, + regularisation_weight=[1,5,50], + points1='', + points2='', + tre_dictionary='results.p', + pre_transform='affine'): + tre=register_images_(im1_fname=im1, + im2_fname=im2, + connectivity=connectivity, + apply_watershed=False, + overwrite=overwrite, + min_object_size=min_object_size, + clip=clip, + fix_rotation=fix_rotation, + mult_factor=mult_factor, + scaling_factor=scaling_factor, + output_dir=output_dir, + properties=properties, + gpu_device=gpu_device, + max_rotation_vertical_px=max_rotation_vertical_px, + loss_fn=loss_fn, + lr=lr, + transform_type=transform_type, + iterations=iterations, + no_segment_analysis=no_segment_analysis, + black_background=black_background, + verbose=verbose, + opt_cm=opt_cm, + sigma=sigma, + order=order, + pyramid=pyramid, + interpolation=interpolation, + half=half, + regularisation_weight=regularisation_weight, + points1=points1, + points2=points2, + pre_transform=pre_transform) + if os.path.exists(tre_dictionary): + tre_dict=pickle.load(open(tre_dictionary,'rb')) + else: + tre_dict=dict() + tre_dict[transform_type,loss_fn]=tre + pickle.dump(tre_dict,open(tre_dictionary,'wb')) def main(): - fire.Fire(Commands) + fire.Fire(Commands) if __name__=='__main__': - main() + main() diff --git a/setup.py b/setup.py index c55b0d1..9475338 100644 --- a/setup.py +++ b/setup.py @@ -13,7 +13,8 @@ 'torch', 'airlab', 'fire', - 'nibabel'] + 'nibabel', + 'pathflowai'] with open('README.md','r', encoding='utf-8') as f: long_description = f.read() @@ -25,7 +26,7 @@ author='Joshua Levy', author_email='joshualevy44@berkeley.edu', license='MIT', - scripts=[], + scripts=['bin/install_mixmatch_dependencies'], entry_points={ 'console_scripts':['pathflow-mixmatch=pathflow_mixmatch.cli:main'] },