From d10f5aada38215026e1187374dd5e2ca19b14cd7 Mon Sep 17 00:00:00 2001 From: Richard Gerum <14153051+rgerum@users.noreply.github.com> Date: Fri, 19 Apr 2024 18:33:21 -0400 Subject: [PATCH] used black to format pyTFM code and removed some linter warnings --- poetry.lock | 88 ++- pyproject.toml | 1 + saenopy/pyTFM/TFM_functions.py | 38 +- saenopy/pyTFM/TFM_tractions.py | 209 ++++++-- saenopy/pyTFM/calculate_deformation.py | 21 +- saenopy/pyTFM/calculate_forces.py | 8 +- saenopy/pyTFM/calculate_strain_energy.py | 23 +- saenopy/pyTFM/calculate_stress.py | 130 +++-- saenopy/pyTFM/correct_stage_drift.py | 16 +- .../pyTFM/graph_theory_for_cell_boundaries.py | 203 ++++--- saenopy/pyTFM/grid_setup_solids_py.py | 503 +++++++++++++----- saenopy/pyTFM/plotting.py | 313 +++++++---- saenopy/pyTFM/stress_functions.py | 200 ++++--- saenopy/pyTFM/suppress_warnings.py | 5 +- saenopy/pyTFM/utilities_TFM.py | 19 +- 15 files changed, 1252 insertions(+), 525 deletions(-) diff --git a/poetry.lock b/poetry.lock index 60d94cc..b757596 100644 --- a/poetry.lock +++ b/poetry.lock @@ -114,6 +114,29 @@ charset-normalizer = ["charset-normalizer"] html5lib = ["html5lib"] lxml = ["lxml"] +[[package]] +name = "black" +version = "24.4.0" +description = "The uncompromising code formatter." +category = "main" +optional = false +python-versions = ">=3.8" + +[package.dependencies] +click = ">=8.0.0" +mypy-extensions = ">=0.4.3" +packaging = ">=22.0" +pathspec = ">=0.9.0" +platformdirs = ">=2" +tomli = {version = ">=1.1.0", markers = "python_version < \"3.11\""} +typing-extensions = {version = ">=4.0.1", markers = "python_version < \"3.11\""} + +[package.extras] +colorama = ["colorama (>=0.4.3)"] +d = ["aiohttp (>=3.7.4)", "aiohttp (>=3.7.4,!=3.9.0)"] +jupyter = ["ipython (>=7.8.0)", "tokenize-rt (>=3.2.0)"] +uvloop = ["uvloop (>=0.15.2)"] + [[package]] name = "bleach" version = "6.1.0" @@ -176,6 +199,17 @@ category = "main" optional = false python-versions = ">=3.7.0" +[[package]] +name = "click" +version = "8.1.7" +description = "Composable command line interface toolkit" +category = "main" +optional = false +python-versions = ">=3.7" + +[package.dependencies] +colorama = {version = "*", markers = "platform_system == \"Windows\""} + [[package]] name = "colorama" version = "0.4.6" @@ -934,6 +968,14 @@ docs = ["sphinx"] gmpy = ["gmpy2 (>=2.1.0a4)"] tests = ["pytest (>=4.6)"] +[[package]] +name = "mypy-extensions" +version = "1.0.0" +description = "Type system extensions for programs checked with the mypy type checker." +category = "main" +optional = false +python-versions = ">=3.5" + [[package]] name = "natsort" version = "8.4.0" @@ -1190,6 +1232,14 @@ python-versions = ">=3.6" qa = ["flake8 (==5.0.4)", "mypy (==0.971)", "types-setuptools (==67.2.0.1)"] testing = ["docopt", "pytest"] +[[package]] +name = "pathspec" +version = "0.12.1" +description = "Utility library for gitignore style pattern matching of file paths." +category = "main" +optional = false +python-versions = ">=3.8" + [[package]] name = "pefile" version = "2023.2.7" @@ -2169,7 +2219,7 @@ docs = ["sphinx", "sphinx-rtd-theme", "nbsphinx", "sphinx-gallery"] [metadata] lock-version = "1.1" python-versions = "^3.8,<3.12" -content-hash = "d589fccf60ad95f0fa107054befc389a3751742c82aee84da2bdf3250f98b98c" +content-hash = "d5fdaaf000076e40eca09fafc8746ef814ff67babe87a2efc42ef12456e3dbed" [metadata.files] alabaster = [ @@ -2212,6 +2262,30 @@ beautifulsoup4 = [ {file = "beautifulsoup4-4.12.3-py3-none-any.whl", hash = "sha256:b80878c9f40111313e55da8ba20bdba06d8fa3969fc68304167741bbf9e082ed"}, {file = "beautifulsoup4-4.12.3.tar.gz", hash = "sha256:74e3d1928edc070d21748185c46e3fb33490f22f52a3addee9aee0f4f7781051"}, ] +black = [ + {file = "black-24.4.0-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:6ad001a9ddd9b8dfd1b434d566be39b1cd502802c8d38bbb1ba612afda2ef436"}, + {file = "black-24.4.0-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:e3a3a092b8b756c643fe45f4624dbd5a389f770a4ac294cf4d0fce6af86addaf"}, + {file = "black-24.4.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:dae79397f367ac8d7adb6c779813328f6d690943f64b32983e896bcccd18cbad"}, + {file = "black-24.4.0-cp310-cp310-win_amd64.whl", hash = "sha256:71d998b73c957444fb7c52096c3843875f4b6b47a54972598741fe9a7f737fcb"}, + {file = "black-24.4.0-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:8e5537f456a22cf5cfcb2707803431d2feeb82ab3748ade280d6ccd0b40ed2e8"}, + {file = "black-24.4.0-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:64e60a7edd71fd542a10a9643bf369bfd2644de95ec71e86790b063aa02ff745"}, + {file = "black-24.4.0-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:5cd5b4f76056cecce3e69b0d4c228326d2595f506797f40b9233424e2524c070"}, + {file = "black-24.4.0-cp311-cp311-win_amd64.whl", hash = "sha256:64578cf99b6b46a6301bc28bdb89f9d6f9b592b1c5837818a177c98525dbe397"}, + {file = "black-24.4.0-cp312-cp312-macosx_10_9_x86_64.whl", hash = "sha256:f95cece33329dc4aa3b0e1a771c41075812e46cf3d6e3f1dfe3d91ff09826ed2"}, + {file = "black-24.4.0-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:4396ca365a4310beef84d446ca5016f671b10f07abdba3e4e4304218d2c71d33"}, + {file = "black-24.4.0-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:44d99dfdf37a2a00a6f7a8dcbd19edf361d056ee51093b2445de7ca09adac965"}, + {file = "black-24.4.0-cp312-cp312-win_amd64.whl", hash = "sha256:21f9407063ec71c5580b8ad975653c66508d6a9f57bd008bb8691d273705adcd"}, + {file = "black-24.4.0-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:652e55bb722ca026299eb74e53880ee2315b181dfdd44dca98e43448620ddec1"}, + {file = "black-24.4.0-cp38-cp38-macosx_11_0_arm64.whl", hash = "sha256:7f2966b9b2b3b7104fca9d75b2ee856fe3fdd7ed9e47c753a4bb1a675f2caab8"}, + {file = "black-24.4.0-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:1bb9ca06e556a09f7f7177bc7cb604e5ed2d2df1e9119e4f7d2f1f7071c32e5d"}, + {file = "black-24.4.0-cp38-cp38-win_amd64.whl", hash = "sha256:d4e71cdebdc8efeb6deaf5f2deb28325f8614d48426bed118ecc2dcaefb9ebf3"}, + {file = "black-24.4.0-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:6644f97a7ef6f401a150cca551a1ff97e03c25d8519ee0bbc9b0058772882665"}, + {file = "black-24.4.0-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:75a2d0b4f5eb81f7eebc31f788f9830a6ce10a68c91fbe0fade34fff7a2836e6"}, + {file = "black-24.4.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:eb949f56a63c5e134dfdca12091e98ffb5fd446293ebae123d10fc1abad00b9e"}, + {file = "black-24.4.0-cp39-cp39-win_amd64.whl", hash = "sha256:7852b05d02b5b9a8c893ab95863ef8986e4dda29af80bbbda94d7aee1abf8702"}, + {file = "black-24.4.0-py3-none-any.whl", hash = "sha256:74eb9b5420e26b42c00a3ff470dc0cd144b80a766128b1771d07643165e08d0e"}, + {file = "black-24.4.0.tar.gz", hash = "sha256:f07b69fda20578367eaebbd670ff8fc653ab181e1ff95d84497f9fa20e7d0641"}, +] bleach = [ {file = "bleach-6.1.0-py3-none-any.whl", hash = "sha256:3225f354cfc436b9789c66c4ee030194bee0568fbf9cbdad3bc8b5c26c5f12b6"}, {file = "bleach-6.1.0.tar.gz", hash = "sha256:0a31f1837963c41d46bbf1331b8778e1308ea0791db03cc4e7357b97cf42a8fe"}, @@ -2373,6 +2447,10 @@ charset-normalizer = [ {file = "charset_normalizer-3.3.2-cp39-cp39-win_amd64.whl", hash = "sha256:b01b88d45a6fcb69667cd6d2f7a9aeb4bf53760d7fc536bf679ec94fe9f3ff3d"}, {file = "charset_normalizer-3.3.2-py3-none-any.whl", hash = "sha256:3e4d1f6587322d2788836a99c69062fbb091331ec940e02d12d179c1d53e25fc"}, ] +click = [ + {file = "click-8.1.7-py3-none-any.whl", hash = "sha256:ae74fb96c20a0277a1d615f1e4d73c8414f5a98db8b799a7931d1582f3390c28"}, + {file = "click-8.1.7.tar.gz", hash = "sha256:ca9853ad459e787e2192211578cc907e7594e294c7ccc834310722b41b9ca6de"}, +] colorama = [ {file = "colorama-0.4.6-py2.py3-none-any.whl", hash = "sha256:4f1d9991f5acc0ca119f9d443620b77f9d6b33703e51011c16baf57afb285fc6"}, {file = "colorama-0.4.6.tar.gz", hash = "sha256:08695f5cb7ed6e0531a20572697297273c47b8cae5a63ffc6d6ed5c201be6e44"}, @@ -3101,6 +3179,10 @@ mpmath = [ {file = "mpmath-1.3.0-py3-none-any.whl", hash = "sha256:a0b2b9fe80bbcd81a6647ff13108738cfb482d481d826cc0e02f5b35e5c88d2c"}, {file = "mpmath-1.3.0.tar.gz", hash = "sha256:7a28eb2a9774d00c7bc92411c19a89209d5da7c4c9a9e227be8330a23a25b91f"}, ] +mypy-extensions = [ + {file = "mypy_extensions-1.0.0-py3-none-any.whl", hash = "sha256:4392f6c0eb8a5668a69e23d168ffa70f0be9ccfd32b5cc2d26a34ae5b844552d"}, + {file = "mypy_extensions-1.0.0.tar.gz", hash = "sha256:75dbf8955dc00442a438fc4d0666508a9a97b6bd41aa2f0ffe9d2f2725af0782"}, +] natsort = [ {file = "natsort-8.4.0-py3-none-any.whl", hash = "sha256:4732914fb471f56b5cce04d7bae6f164a592c7712e1c85f9ef585e197299521c"}, {file = "natsort-8.4.0.tar.gz", hash = "sha256:45312c4a0e5507593da193dedd04abb1469253b601ecaf63445ad80f0a1ea581"}, @@ -3236,6 +3318,10 @@ parso = [ {file = "parso-0.8.4-py2.py3-none-any.whl", hash = "sha256:a418670a20291dacd2dddc80c377c5c3791378ee1e8d12bffc35420643d43f18"}, {file = "parso-0.8.4.tar.gz", hash = "sha256:eb3a7b58240fb99099a345571deecc0f9540ea5f4dd2fe14c2a99d6b281ab92d"}, ] +pathspec = [ + {file = "pathspec-0.12.1-py3-none-any.whl", hash = "sha256:a0d503e138a4c123b27490a4f7beda6a01c6f288df0e4a8b79c7eb0dc7b4cc08"}, + {file = "pathspec-0.12.1.tar.gz", hash = "sha256:a482d51503a1ab33b1c67a6c3813a26953dbdc71c31dacaef9a838c4e29f5712"}, +] pefile = [ {file = "pefile-2023.2.7-py3-none-any.whl", hash = "sha256:da185cd2af68c08a6cd4481f7325ed600a88f6a813bad9dea07ab3ef73d8d8d6"}, {file = "pefile-2023.2.7.tar.gz", hash = "sha256:82e6114004b3d6911c77c3953e3838654b04511b8b66e8583db70c65998017dc"}, diff --git a/pyproject.toml b/pyproject.toml index 224cdac..34e7d41 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -34,6 +34,7 @@ sphinx-rtd-theme = { version = "^1.2.2", optional = true } nbsphinx = { version = "^0.9.2", optional = true } sphinx-gallery = {version = "^0.13.0", optional = true } solidspy = "*" +black = { version="^24.4.0", optional = true } [tool.poetry.extras] docs = ["sphinx", "sphinx-rtd-theme", "nbsphinx", "sphinx-gallery"] diff --git a/saenopy/pyTFM/TFM_functions.py b/saenopy/pyTFM/TFM_functions.py index d2a1828..d15e02c 100644 --- a/saenopy/pyTFM/TFM_functions.py +++ b/saenopy/pyTFM/TFM_functions.py @@ -1,18 +1,12 @@ import numpy as np -import openpiv.filters -from openpiv.pyprocess import extended_search_area_piv -import openpiv.scaling -import openpiv.tools -import openpiv.validation - -def strain_energy_points(u, v, tx, ty, pixelsize1, pixelsize2): - pixelsize1 *= 10 ** -6 - pixelsize2 *= 10 ** -6 # conversion to m - # u is given in pixels/minutes where a pixel is from the original image (pixelsize1) - # tx is given in forces/pixels**2 where a pixel is from the deformation/traction field (pixelsize2) - energy_points = 0.5 * (pixelsize2 ** 2) * (tx * u * pixelsize1 + ty * v * pixelsize1) +def strain_energy_points(u, v, tx, ty, pixel_size1, pixel_size2): + pixel_size1 *= 10 ** -6 + pixel_size2 *= 10 ** -6 # conversion to m + # u is given in pixels/minutes where a pixel is from the original image (pixel_size1) + # tx is given in forces/pixels**2 where a pixel is from the deformation/traction field (pixel_size2) + energy_points = 0.5 * (pixel_size2 ** 2) * (tx * u * pixel_size1 + ty * v * pixel_size1) # value of a background point bg = np.percentile(energy_points, 20) energy_points -= bg @@ -34,14 +28,13 @@ def get_xy_for_quiver(u): return xs, ys -def contractillity(tx, ty, pixelsize, mask): - +def contractility(tx, ty, pixel_size, mask): """ Calculation of contractile force and force epicenter.Contractile force is the sum of all projection of traction forces (in N) towards the force epicenter. The force epicenter is the point that maximizes the contractile force. :param tx: traction forces in x direction in Pa :param ty: traction forces in y direction in Pa - :param pixelsize: pixelsize of the traction field + :param pixel_size: pixel size of the traction field :param mask: mask of which values to use for calculation :return: contractile_force,contractile force in N proj_x, projection of traction forces towards the force epicenter, x component @@ -56,19 +49,19 @@ def contractillity(tx, ty, pixelsize, mask): ty_filter = np.zeros(np.shape(ty)) ty_filter[mask] = ty[mask] - tract_abs = np.sqrt(tx_filter ** 2 + ty_filter ** 2) + tract_abs = np.sqrt(tx_filter**2 + ty_filter**2) - area = (pixelsize * (10 ** -6)) ** 2 # in meter + area = (pixel_size * (10 ** -6)) ** 2 # in meter fx = tx_filter * area # calculating forces (in Newton) by multiplying with area fy = ty_filter * area x, y = get_xy_for_quiver(tx) - bx = np.sum(x * (tract_abs ** 2) + fx * (tx_filter * fx + ty_filter * fy)) - by = np.sum(y * (tract_abs ** 2) + fy * (tx_filter * fx + ty_filter * fy)) + bx = np.sum(x * (tract_abs**2) + fx * (tx_filter * fx + ty_filter * fy)) + by = np.sum(y * (tract_abs**2) + fy * (tx_filter * fx + ty_filter * fy)) - axx = np.sum(tract_abs ** 2 + fx ** 2) + axx = np.sum(tract_abs**2 + fx**2) axy = np.sum(fx * fy) - ayy = np.sum(tract_abs ** 2 + fy ** 2) + ayy = np.sum(tract_abs**2 + fy**2) # ayx=np.sum(tx*ty) A = np.array([[axx, axy], [axy, ayy]]) @@ -85,7 +78,7 @@ def contractillity(tx, ty, pixelsize, mask): dist_x = center[0] - x dist_y = center[1] - y - dist_abs = np.sqrt(dist_y ** 2 + dist_x ** 2) + dist_abs = np.sqrt(dist_y**2 + dist_x**2) proj_abs = (fx * dist_x + fy * dist_y) / dist_abs contractile_force = np.nansum(proj_abs) @@ -94,4 +87,3 @@ def contractillity(tx, ty, pixelsize, mask): proj_y = proj_abs * dist_y / dist_abs return contractile_force, proj_x, proj_y, center # unit of contractile force is N - diff --git a/saenopy/pyTFM/TFM_tractions.py b/saenopy/pyTFM/TFM_tractions.py index ec7b41e..0453d9b 100644 --- a/saenopy/pyTFM/TFM_tractions.py +++ b/saenopy/pyTFM/TFM_tractions.py @@ -6,10 +6,12 @@ from .suppress_warnings import suppress_warnings -def ffttc_traction(u, v, pixelsize1, pixelsize2, young, sigma=0.49, spatial_filter="gaussian", fs=None): +def ffttc_traction( + u, v, pixel_size1, pixel_size2, young, sigma=0.49, spatial_filter="gaussian", fs=None +): """ fourier transform based calculation of the traction force. U and v must be given as deformations in pixel. Size of - these pixels must be the pixelsize (size of a pixel in the deformation field u or v). Note that thePiv deformation + these pixels must be the pixel size (size of a pixel in the deformation field u or v). Note that thePiv deformation returns deformation in pixel of the size of pixels in the images of beads before and after. If bf_image is provided this script will return a traction field that is zoomed to the size of the bright field image, by interpolation. It is not recommended to use this for any calculations. @@ -18,17 +20,23 @@ def ffttc_traction(u, v, pixelsize1, pixelsize2, young, sigma=0.49, spatial_filt :param u:deformation field in x direction in pixel of the deformation image :param v:deformation field in y direction in pixel of the deformation image :param young: Young's modulus in Pa - :param pixelsize1: pixelsize in m/pixel of the original image, needed because u and v is given as + :param pixel_size1: pixel size in m/pixel of the original image, needed because u and v is given as displacement of these pixels - :param pixelsize2: pixelsize of m/pixel the deformation image + :param pixel_size2: pixel size of m/pixel the deformation image :param sigma: poisson ratio of the gel :param spatial_filter: str, values: "mean","gaussian","median". Different smoothing methods for the traction field :return: tx_filter,ty_filter: traction forces in x and y direction in Pa + + Parameters + ---------- + fs """ # 0) subtracting mean(better name for this step) - u_shift = (u - np.mean(u)) # shifting to zero mean (translating to pixelsize of u-image is done later) - v_shift = (v - np.mean(v)) + u_shift = u - np.mean( + u + ) # shifting to zero mean (translating to pixel size of u-image is done later) + v_shift = v - np.mean(v) # Ben's algorithm: # 1)Zero padding to get square array with even index number @@ -45,35 +53,51 @@ def ffttc_traction(u, v, pixelsize1, pixelsize2, young, sigma=0.49, spatial_filt # 2) producing wave vectors # form 1:max_ind/2 then -(max_ind/2:1) - kx1 = np.array([list(range(0, int(max_ind / 2), 1)), ] * int(max_ind)) - kx2 = np.array([list(range(-int(max_ind / 2), 0, 1)), ] * int(max_ind)) - kx = np.append(kx1, kx2, axis=1) * 2 * np.pi # fourier transform in this case is defined as + kx1 = np.array( + [ + list(range(0, int(max_ind / 2), 1)), + ] + * int(max_ind) + ) + kx2 = np.array( + [ + list(range(-int(max_ind / 2), 0, 1)), + ] + * int(max_ind) + ) + kx = ( + np.append(kx1, kx2, axis=1) * 2 * np.pi + ) # fourier transform in this case is defined as # F(kx)=1/2pi integral(exp(i*kx*x)dk therefore kx must be expressed as a spatial frequency in distance*2*pi ky = np.transpose(kx) - k = np.sqrt(kx ** 2 + ky ** 2) / (pixelsize2 * max_ind) + k = np.sqrt(kx**2 + ky**2) / (pixel_size2 * max_ind) # 2.1) calculating angle between k and kx with atan2 function alpha = np.arctan2(ky, kx) alpha[0, 0] = np.pi / 2 - # 3) calculation of K --> Tensor to calculate displacements from tractions. We calculate inverse of K + # 3) calculation of K --> Tensor to calculate displacements from traction forces. We calculate inverse of K # K⁻¹=[[kix kid], # [kid,kiy]] ,,, so is "diagonal, kid appears two times - kix = ((k * young) / (2 * (1 - sigma ** 2))) * (1 - sigma + sigma * np.cos(alpha) ** 2) - kiy = ((k * young) / (2 * (1 - sigma ** 2))) * (1 - sigma + sigma * np.sin(alpha) ** 2) - kid = ((k * young) / (2 * (1 - sigma ** 2))) * (sigma * np.sin(alpha) * np.cos(alpha)) + kix = ((k * young) / (2 * (1 - sigma**2))) * ( + 1 - sigma + sigma * np.cos(alpha) ** 2 + ) + kiy = ((k * young) / (2 * (1 - sigma**2))) * ( + 1 - sigma + sigma * np.sin(alpha) ** 2 + ) + kid = ((k * young) / (2 * (1 - sigma**2))) * (sigma * np.sin(alpha) * np.cos(alpha)) # adding zeros in kid diagonals kid[:, int(max_ind / 2)] = np.zeros(max_ind) kid[int(max_ind / 2), :] = np.zeros(max_ind) # 4) calculate Fourier transform of displacement - # u_ft=np.fft.fft2(u_expand*pixelsize1*2*np.pi) - # v_ft=np.fft.fft2(v_expand*pixelsize1*2*np.pi) - u_ft = scipy.fft.fft2(u_expand * pixelsize1) - v_ft = scipy.fft.fft2(v_expand * pixelsize1) + # u_ft=np.fft.fft2(u_expand*pixel_size1*2*np.pi) + # v_ft=np.fft.fft2(v_expand*pixel_size1*2*np.pi) + u_ft = scipy.fft.fft2(u_expand * pixel_size1) + v_ft = scipy.fft.fft2(v_expand * pixel_size1) - # 4.1) calculate tractions in Fourier space T=K⁻¹*U, U=[u,v] here with individual matrix elements.. + # 4.1) calculate traction forces in Fourier space T=K⁻¹*U, U=[u,v] here with individual matrix elements. tx_ft = kix * u_ft + kid * v_ft ty_ft = kid * u_ft + kiy * v_ft @@ -90,22 +114,44 @@ def ffttc_traction(u, v, pixelsize1, pixelsize2, young, sigma=0.49, spatial_filt ty_filter = ty_cut if spatial_filter == "mean": - fs = fs if isinstance(fs, (float, int)) else int(int(np.max((ax1_length, ax2_length))) / 16) + fs = ( + fs + if isinstance(fs, (float, int)) + else int(int(np.max((ax1_length, ax2_length))) / 16) + ) tx_filter = uniform_filter(tx_cut, size=fs) ty_filter = uniform_filter(ty_cut, size=fs) if spatial_filter == "gaussian": - fs = fs if isinstance(fs, (float, int)) else int(np.max((ax1_length, ax2_length))) / 50 + fs = ( + fs + if isinstance(fs, (float, int)) + else int(np.max((ax1_length, ax2_length))) / 50 + ) tx_filter = gaussian_filter(tx_cut, sigma=fs) ty_filter = gaussian_filter(ty_cut, sigma=fs) if spatial_filter == "median": - fs = fs if isinstance(fs, (float, int)) else int(int(np.max((ax1_length, ax2_length))) / 16) + fs = ( + fs + if isinstance(fs, (float, int)) + else int(int(np.max((ax1_length, ax2_length))) / 16) + ) tx_filter = median_filter(tx_cut, size=fs) ty_filter = median_filter(ty_cut, size=fs) return tx_filter, ty_filter -def ffttc_traction_finite_thickness(u, v, pixelsize1, pixelsize2, h, young, sigma=0.49, spatial_filter="gaussian", fs=None): +def ffttc_traction_finite_thickness( + u, + v, + pixel_size1, + pixel_size2, + h, + young, + sigma=0.49, + spatial_filter="gaussian", + fs=None, +): """ FTTC with correction for finite substrate thickness according to Xavier Trepat, Physical forces during collective cell migration, 2009 @@ -113,8 +159,8 @@ def ffttc_traction_finite_thickness(u, v, pixelsize1, pixelsize2, h, young, sigm :param u:deformation field in x direction in pixel of the deformation image :param v:deformation field in y direction in pixel of the deformation image :param young: Young's modulus in Pa - :param pixelsize1: pixelsize of the original image, needed because u and v is given as displacement of these pixels - :param pixelsize2: pixelsize of the deformation image + :param pixel_size1: pixel size of the original image, needed because u and v is given as displacement of these pixels + :param pixel_size2: pixel size of the deformation image :param h: height of the membrane the cells lie on, in µm :param sigma: Poisson's ratio of the gel :param spatial_filter: str, values: "mean","gaussian","median". Different smoothing methods for the traction field. @@ -123,10 +169,10 @@ def ffttc_traction_finite_thickness(u, v, pixelsize1, pixelsize2, h, young, sigm """ # 0) subtracting mean(better name for this step) - u_shift = (u - np.mean(u)) * pixelsize1 - v_shift = (v - np.mean(v)) * pixelsize1 + u_shift = (u - np.mean(u)) * pixel_size1 + v_shift = (v - np.mean(v)) * pixel_size1 - # Ben's algortithm: + # Ben's algorithm: # 1)Zero padding to get square array with even index number ax1_length = np.shape(u_shift)[0] # u and v must have same dimensions ax2_length = np.shape(u_shift)[1] @@ -135,18 +181,30 @@ def ffttc_traction_finite_thickness(u, v, pixelsize1, pixelsize2, h, young, sigm max_ind += 1 u_expand = np.zeros((max_ind, max_ind)) v_expand = np.zeros((max_ind, max_ind)) - u_expand[max_ind - ax1_length:max_ind, max_ind - ax2_length:max_ind] = u_shift - v_expand[max_ind - ax1_length:max_ind, max_ind - ax2_length:max_ind] = v_shift + u_expand[max_ind - ax1_length : max_ind, max_ind - ax2_length : max_ind] = u_shift + v_expand[max_ind - ax1_length : max_ind, max_ind - ax2_length : max_ind] = v_shift # 2) producing wave vectors (FT-space) # form 1:max_ind/2 then -(max_ind/2:1) - kx1 = np.array([list(range(0, int(max_ind / 2), 1)), ] * int(max_ind), dtype=np.float64) - kx2 = np.array([list(range(-int(max_ind / 2), 0, 1)), ] * int(max_ind), dtype=np.float64) + kx1 = np.array( + [ + list(range(0, int(max_ind / 2), 1)), + ] + * int(max_ind), + dtype=np.float64, + ) + kx2 = np.array( + [ + list(range(-int(max_ind / 2), 0, 1)), + ] + * int(max_ind), + dtype=np.float64, + ) # spatial frequencies: 1/wavelength,in 1/µm in fractions of total length - kx = np.append(kx1, kx2, axis=1) * 2 * np.pi / (pixelsize2 * max_ind) + kx = np.append(kx1, kx2, axis=1) * 2 * np.pi / (pixel_size2 * max_ind) ky = np.transpose(kx) - k = np.sqrt(kx ** 2 + ky ** 2) # matrix with "relative" distances??# + k = np.sqrt(kx**2 + ky**2) # matrix with "relative" distances??# r = k * h c = np.cosh(r) @@ -155,8 +213,9 @@ def ffttc_traction_finite_thickness(u, v, pixelsize1, pixelsize2, h, young, sigm # gamma = ((3 - 4 * sigma) * (c ** 2) + (1 - 2 * sigma) ** 2 + (k * h) ** 2) / ( # (3 - 4 * sigma) * s * c + k * h) ## inf values here because k goes to zero - gamma = ((3 - 4 * sigma) + (((1 - 2 * sigma) ** 2) / (c ** 2)) + ((r ** 2) / (c ** 2))) / ( - (3 - 4 * sigma) * s_c + r / (c ** 2)) + gamma = ( + (3 - 4 * sigma) + (((1 - 2 * sigma) ** 2) / (c**2)) + ((r**2) / (c**2)) + ) / ((3 - 4 * sigma) * s_c + r / (c**2)) # 4) calculate fourier transform of displacements u_ft = scipy.fft.fft2(u_expand) @@ -170,14 +229,16 @@ def ffttc_traction_finite_thickness(u, v, pixelsize1, pixelsize2, h, young, sigm ty_ft=factor3*(v_ft*((k**2)*(1-sigma)+sigma*(ky**2)) + u_ft*kx*ky*sigma) """ - # 4.1) calculate tractions in Fourier space - factor1 = (v_ft * kx - u_ft * ky) - factor2 = (u_ft * kx + v_ft * ky) + # 4.1) calculate traction forces in Fourier space + factor1 = v_ft * kx - u_ft * ky + factor2 = u_ft * kx + v_ft * ky tx_ft = ((-young * ky * c) / (2 * k * s * (1 + sigma))) * factor1 + ( - (young * kx) / (2 * k * (1 - sigma ** 2))) * gamma * factor2 + (young * kx) / (2 * k * (1 - sigma**2)) + ) * gamma * factor2 tx_ft[0, 0] = 0 # zero frequency would represent force everywhere (constant) ty_ft = ((young * kx * c) / (2 * k * s * (1 + sigma))) * factor1 + ( - (young * ky) / (2 * k * (1 - sigma ** 2))) * gamma * factor2 + (young * ky) / (2 * k * (1 - sigma**2)) + ) * gamma * factor2 ty_ft[0, 0] = 0 # 4.2) go back to real space @@ -185,29 +246,43 @@ def ffttc_traction_finite_thickness(u, v, pixelsize1, pixelsize2, h, young, sigm ty = scipy.fft.ifft2(ty_ft.astype(np.complex128)).real # 5.2) cut like in script from ben - tx_cut = tx[max_ind - ax1_length:max_ind, max_ind - ax2_length:max_ind] - ty_cut = ty[max_ind - ax1_length:max_ind, max_ind - ax2_length:max_ind] + tx_cut = tx[max_ind - ax1_length : max_ind, max_ind - ax2_length : max_ind] + ty_cut = ty[max_ind - ax1_length : max_ind, max_ind - ax2_length : max_ind] # 5.3) using filter tx_filter = tx_cut ty_filter = ty_cut if spatial_filter == "mean": - fs = fs if isinstance(fs, (float, int)) else int(int(np.max((ax1_length, ax2_length))) / 16) + fs = ( + fs + if isinstance(fs, (float, int)) + else int(int(np.max((ax1_length, ax2_length))) / 16) + ) tx_filter = uniform_filter(tx_cut, size=fs) ty_filter = uniform_filter(ty_cut, size=fs) if spatial_filter == "gaussian": - fs = fs if isinstance(fs, (float, int)) else int(np.max((ax1_length, ax2_length))) / 50 + fs = ( + fs + if isinstance(fs, (float, int)) + else int(np.max((ax1_length, ax2_length))) / 50 + ) tx_filter = gaussian_filter(tx_cut, sigma=fs) ty_filter = gaussian_filter(ty_cut, sigma=fs) if spatial_filter == "median": - fs = fs if isinstance(fs, (float, int)) else int(int(np.max((ax1_length, ax2_length))) / 16) + fs = ( + fs + if isinstance(fs, (float, int)) + else int(int(np.max((ax1_length, ax2_length))) / 16) + ) tx_filter = median_filter(tx_cut, size=fs) ty_filter = median_filter(ty_cut, size=fs) return tx_filter, ty_filter -def TFM_tractions(u, v, pixelsize1, pixelsize2, h, young, sigma=0.49, spatial_filter="gaussian", fs=6): +def TFM_tractions( + u, v, pixel_size1, pixel_size2, h, young, sigma=0.49, spatial_filter="gaussian", fs=6 +): """ height correction breaks down due to numerical reasons at large gel height and small wavelengths of deformations. In this case the height corrected ffttc-function returns Nans. THis function falls back @@ -215,19 +290,43 @@ def TFM_tractions(u, v, pixelsize1, pixelsize2, h, young, sigma=0.49, spatial_fi :return: """ # translate the filter size to pixels of traction field - fs = fs / pixelsize2 if isinstance(fs, (int, float)) else None + fs = fs / pixel_size2 if isinstance(fs, (int, float)) else None if isinstance(h, (int, float)): with suppress_warnings(RuntimeWarning): - tx, ty = ffttc_traction_finite_thickness(u, v, pixelsize1=pixelsize1, pixelsize2=pixelsize2, - h=h, young=young, sigma=sigma, - spatial_filter=spatial_filter, fs=fs) # unit is N/m**2 + tx, ty = ffttc_traction_finite_thickness( + u, + v, + pixel_size1=pixel_size1, + pixel_size2=pixel_size2, + h=h, + young=young, + sigma=sigma, + spatial_filter=spatial_filter, + fs=fs, + ) # unit is N/m**2 if np.any(np.isnan(tx)) or np.any(np.isnan(ty)): - tx, ty = ffttc_traction(u, v, pixelsize1=pixelsize1, pixelsize2=pixelsize2, young=young, sigma=sigma, - spatial_filter=spatial_filter, fs=fs) + tx, ty = ffttc_traction( + u, + v, + pixel_size1=pixel_size1, + pixel_size2=pixel_size2, + young=young, + sigma=sigma, + spatial_filter=spatial_filter, + fs=fs, + ) elif h == "infinite": - tx, ty = ffttc_traction(u, v, pixelsize1=pixelsize1, pixelsize2=pixelsize2, young=young, sigma=sigma, - spatial_filter=spatial_filter, fs=fs) + tx, ty = ffttc_traction( + u, + v, + pixel_size1=pixel_size1, + pixel_size2=pixel_size2, + young=young, + sigma=sigma, + spatial_filter=spatial_filter, + fs=fs, + ) else: raise ValueError("illegal value for h") return tx, ty diff --git a/saenopy/pyTFM/calculate_deformation.py b/saenopy/pyTFM/calculate_deformation.py index 44205db..9eab92e 100644 --- a/saenopy/pyTFM/calculate_deformation.py +++ b/saenopy/pyTFM/calculate_deformation.py @@ -40,15 +40,20 @@ def calculate_deformation(im1, im2, window_size=64, overlap=32, std_factor=20): else: raise ValueError - u, v, sig2noise = extended_search_area_piv(frame_a, frame_b, window_size=window_size, - overlap=overlap, - dt=1, subpixel_method="gaussian", - search_area_size=window_size, - sig2noise_method='peak2peak') + u, v, sig2noise = extended_search_area_piv( + frame_a, + frame_b, + window_size=window_size, + overlap=overlap, + dt=1, + subpixel_method="gaussian", + search_area_size=window_size, + sig2noise_method="peak2peak", + ) u, v, mask = openpiv.validation.sig2noise_val(u, v, sig2noise, threshold=1.05) - def_abs = np.sqrt(u ** 2 + v ** 2) + def_abs = np.sqrt(u**2 + v**2) m = np.nanmean(def_abs) std = np.nanstd(def_abs) @@ -57,5 +62,7 @@ def calculate_deformation(im1, im2, window_size=64, overlap=32, std_factor=20): u[mask_std] = np.nan v[mask_std] = np.nan - u, v = openpiv.filters.replace_outliers(u, v, method='localmean', max_iter=10, kernel_size=2) + u, v = openpiv.filters.replace_outliers( + u, v, method="localmean", max_iter=10, kernel_size=2 + ) return u, -v, mask, mask_std # return -v because of image inverted axis diff --git a/saenopy/pyTFM/calculate_forces.py b/saenopy/pyTFM/calculate_forces.py index 03dcf95..ac2de9c 100644 --- a/saenopy/pyTFM/calculate_forces.py +++ b/saenopy/pyTFM/calculate_forces.py @@ -1,11 +1,15 @@ import numpy as np from saenopy.pyTFM.TFM_tractions import TFM_tractions + def calculate_forces(u, v, pixel_size, shape, h, young, sigma): ps1 = pixel_size # pixel size of the image of the beads # dimensions of the image of the beads im1_shape = shape ps2 = ps1 * np.mean( - np.array(im1_shape) / np.array(u.shape)) # pixel size of the deformation field - tx, ty = TFM_tractions(u, v, pixelsize1=ps1, pixelsize2=ps2, h=h, young=young, sigma=sigma) + np.array(im1_shape) / np.array(u.shape) + ) # pixel size of the deformation field + tx, ty = TFM_tractions( + u, v, pixel_size1=ps1, pixel_size2=ps2, h=h, young=young, sigma=sigma + ) return tx, ty diff --git a/saenopy/pyTFM/calculate_strain_energy.py b/saenopy/pyTFM/calculate_strain_energy.py index d606a12..d1e2636 100644 --- a/saenopy/pyTFM/calculate_strain_energy.py +++ b/saenopy/pyTFM/calculate_strain_energy.py @@ -1,19 +1,26 @@ import numpy as np + try: from scipy.ndimage import binary_fill_holes except ImportError: from scipy.ndimage.morphology import binary_fill_holes -from saenopy.pyTFM.TFM_functions import strain_energy_points, contractillity -from saenopy.pyTFM.grid_setup_solids_py import interpolation # a simple function to resize the mask +from saenopy.pyTFM.TFM_functions import strain_energy_points, contractility +from saenopy.pyTFM.grid_setup_solids_py import ( + interpolation, +) # a simple function to resize the mask def calculate_strain_energy(mask, pixel_size, shape, u, v, tx, ty): - mask = binary_fill_holes(mask == 1) # the mask should be a single patch without holes + mask = binary_fill_holes( + mask == 1 + ) # the mask should be a single patch without holes # changing the masks dimensions to fit to the deformation and traction fields mask = interpolation(mask, dims=u.shape) ps1 = pixel_size # pixel size of the image of the beads # dimensions of the image of the beads - ps2 = ps1 * np.mean(np.array(shape) / np.array(u.shape)) # pixel size of the deformation field + ps2 = ps1 * np.mean( + np.array(shape) / np.array(u.shape) + ) # pixel size of the deformation field # strain energy: # first we calculate a map of strain energy energy_points = strain_energy_points(u, v, tx, ty, ps1, ps2) # J/pixel @@ -21,11 +28,13 @@ def calculate_strain_energy(mask, pixel_size, shape, u, v, tx, ty): # then we sum all energy points in the area defined by mask strain_energy = np.sum(energy_points[mask]) # 2.14*10**-13 J # contractility - contractile_force, proj_x, proj_y, center = contractillity(tx, ty, ps2, mask) # 2.03*10**-6 N + contractile_force, proj_x, proj_y, center = contractility( + tx, ty, ps2, mask + ) # 2.03*10**-6 N return { "contractility": contractile_force, - "area Traction Area": np.sum(mask) * ((pixel_size * 10 ** -6) ** 2), + "area Traction Area": np.sum(mask) * ((pixel_size * 10**-6) ** 2), "strain energy": strain_energy, "center of object": center, - } \ No newline at end of file + } diff --git a/saenopy/pyTFM/calculate_stress.py b/saenopy/pyTFM/calculate_stress.py index 9c2242e..b10b28e 100644 --- a/saenopy/pyTFM/calculate_stress.py +++ b/saenopy/pyTFM/calculate_stress.py @@ -1,25 +1,17 @@ import numpy as np -import matplotlib.pyplot as plt -from qtpy import QtWidgets -from saenopy.gui.common import QtShortCuts -from tifffile import imread + try: from scipy.ndimage import binary_fill_holes except ImportError: from scipy.ndimage.morphology import binary_fill_holes -from typing import List, Tuple - -from saenopy.gui.common.gui_classes import CheckAbleGroup -from saenopy.gui.common.code_export import get_code -from saenopy.pyTFM.grid_setup_solids_py import interpolation # a simple function to resize the mask +from saenopy.pyTFM.grid_setup_solids_py import ( + interpolation, +) # a simple function to resize the mask from saenopy.pyTFM.grid_setup_solids_py import prepare_forces from saenopy.pyTFM.grid_setup_solids_py import grid_setup, FEM_simulation from saenopy.pyTFM.grid_setup_solids_py import find_borders - -from saenopy.pyTFM.plotting import plot_continuous_boundary_stresses - -from saenopy.pyTFM.stress_functions import lineTension +from saenopy.pyTFM.stress_functions import line_tension from saenopy.pyTFM.stress_functions import all_stress_measures, coefficient_of_variation from saenopy.pyTFM.stress_functions import mean_stress_vector_norm @@ -29,51 +21,68 @@ def calculate_stress(mask, pixel_size, shape, u, tx, ty): ps1 = pixel_size # pixel size of the image of the beads # dimensions of the image of the beads - ps2 = ps1 * np.mean(np.array(shape) / np.array(u.shape)) # pixel size of the deformation field + ps2 = ps1 * np.mean( + np.array(shape) / np.array(u.shape) + ) # pixel size of the deformation field # first mask: The area used for Finite Elements Methods # it should encircle all forces generated by the cell colony - mask_FEM = binary_fill_holes(mask == 1) # the mask should be a single patch without holes + mask_fem = binary_fill_holes( + mask == 1 + ) # the mask should be a single patch without holes # changing the masks dimensions to fit to the deformation and traction field: - mask_FEM = interpolation(mask_FEM, dims=tx.shape) + mask_fem = interpolation(mask_fem, dims=tx.shape) # second mask: The area of the cells. Average stresses and other values are calculated only # on the actual area of the cell, represented by this mask. mask_cells_full = binary_fill_holes(mask == 2) mask_cells = interpolation(mask_cells_full, dims=tx.shape) - # converting tractions (forces per surface area) to actual forces + # converting traction forces (forces per surface area) to actual forces # and correcting imbalanced forces and torques # tx->traction forces in x direction, ty->traction forces in y direction # ps2->pixel size of the traction field, mask_FEM-> mask for FEM - fx, fy = prepare_forces(tx, ty, ps2, mask_FEM) + fx, fy = prepare_forces(tx, ty, ps2, mask_fem) # constructing the FEM grid - nodes, elements, loads, mats = grid_setup(mask_FEM, -fx, -fy, sigma=0.5) + nodes, elements, loads, mats = grid_setup(mask_fem, -fx, -fy, sigma=0.5) # performing the FEM analysis # verbose prints the progress of numerically solving the FEM system of equations. - UG_sol, stress_tensor = FEM_simulation(nodes, elements, loads, mats, mask_FEM, verbose=True) + ug_sol, stress_tensor = FEM_simulation( + nodes, elements, loads, mats, mask_fem, verbose=True + ) stress_tensor = np.nan_to_num(stress_tensor) # UG_sol is a list of deformations for each node. We don't need it here. # analyzing the FEM results with average stresses - sigma_max, sigma_min, sigma_max_abs, tau_max, phi_n, phi_shear, sigma_mean = all_stress_measures(stress_tensor, - px_size=ps2 * 10 ** -6) + sigma_max, sigma_min, sigma_max_abs, tau_max, phi_n, phi_shear, sigma_mean = ( + all_stress_measures(stress_tensor, px_size=ps2 * 10**-6) + ) mask_int = interpolation(mask_cells_full, dims=sigma_mean.shape, min_cell_size=100) res_dict["mean normal stress Cell Area"] = np.mean(np.abs(sigma_mean[mask_int])) res_dict["max normal stress Cell Area"] = np.mean(np.abs(sigma_max_abs[mask_int])) res_dict["max shear stress Cell Area"] = np.mean(np.abs(tau_max[mask_int])) - res_dict["cv mean normal stress Cell Area"] = coefficient_of_variation(mask_int, sigma_mean, 0) - res_dict["cv max normal stress Cell Area"] = coefficient_of_variation(mask_int, sigma_max_abs, 0) - res_dict["cv max shear stress Cell Area"] = coefficient_of_variation(mask_int, tau_max, 0) + res_dict["cv mean normal stress Cell Area"] = coefficient_of_variation( + mask_int, sigma_mean, 0 + ) + res_dict["cv max normal stress Cell Area"] = coefficient_of_variation( + mask_int, sigma_max_abs, 0 + ) + res_dict["cv max shear stress Cell Area"] = coefficient_of_variation( + mask_int, tau_max, 0 + ) # mean normal stress - ms_map = ((stress_tensor[:, :, 0, 0] + stress_tensor[:, :, 1, 1]) / 2) / (ps2 * 10 ** -6) + ms_map = ((stress_tensor[:, :, 0, 0] + stress_tensor[:, :, 1, 1]) / 2) / ( + ps2 * 10**-6 + ) # average on the area of the cell colony. ms = np.mean(ms_map[mask_cells]) # 0.0043 N/m # coefficient of variation - cv = np.nanstd(ms_map[mask_cells]) / np.abs(np.nanmean(ms_map[mask_cells])) # 0.41 no unit + cv = np.nanstd(ms_map[mask_cells]) / np.abs( + np.nanmean(ms_map[mask_cells]) + ) # 0.41 no unit """ Calculating the Line Tension """ # identifying borders, counting cells, performing spline interpolation to smooth the borders @@ -81,71 +90,94 @@ def calculate_stress(mask, pixel_size, shape, u, tx, ty): # we can for example get the number of cells from the "borders" object n_cells = borders.n_cells # 8 res_dict["cell number"] = n_cells - res_dict["area Cell Area"] = np.sum(mask == 2) * ((pixel_size * 10 ** -6) ** 2) + res_dict["area Cell Area"] = np.sum(mask == 2) * ((pixel_size * 10**-6) ** 2) # calculating the line tension along the cell borders - lt, min_v, max_v = lineTension(borders.lines_splines, borders.line_lengths, stress_tensor, pixel_length=ps2) + lt, min_v, max_v = line_tension( + borders.lines_splines, borders.line_lengths, stress_tensor, pixel_length=ps2 + ) # lt is a nested dictionary. The first key is the id of a cell border. # For each cell border the line tension vectors ("t_vecs"), the normal # and shear component of the line tension ("t_shear") and the normal # vectors of the cell border ("n_vecs") are calculated at a large number of points. # average norm of the line tension. Only borders not at colony edge are used - lt_vecs = np.concatenate([lt[l_id]["t_vecs"] for l_id in lt.keys() if l_id not in borders.edge_lines]) - avg_line_tension = np.mean(np.linalg.norm(lt_vecs, axis=1)) # 0.00569 N/m - res_dict["average magnitude line tension"] = np.mean(np.linalg.norm(lt_vecs, axis=1)) + lt_vecs = np.concatenate( + [lt[l_id]["t_vecs"] for l_id in lt.keys() if l_id not in borders.edge_lines] + ) + + res_dict["average magnitude line tension"] = np.mean( + np.linalg.norm(lt_vecs, axis=1) + ) res_dict["std magnitude line tension"] = np.std(np.linalg.norm(lt_vecs, axis=1)) # average normal component of the line tension - lt_normal = np.concatenate([lt[l_id]["t_normal"] for l_id in lt.keys() if l_id not in borders.edge_lines]) - avg_normal_line_tension = np.mean(np.abs(lt_normal)) # 0.00566 N/m, + lt_normal = np.concatenate( + [lt[l_id]["t_normal"] for l_id in lt.keys() if l_id not in borders.edge_lines] + ) + # here you can see that almost the line tensions act almost exclusively perpendicular to the cell borders. res_dict["average normal line tension"] = np.mean(np.abs(lt_normal)) res_dict["std normal line tension"] = np.std(np.abs(lt_normal)) # average normal component of the line tension - lt_shear = np.concatenate([lt[l_id]["t_shear"] for l_id in lt.keys() if l_id not in borders.edge_lines]) + lt_shear = np.concatenate( + [lt[l_id]["t_shear"] for l_id in lt.keys() if l_id not in borders.edge_lines] + ) res_dict["average shear line tension"] = np.mean(np.abs(lt_shear)) res_dict["std shear line tension"] = np.std(np.abs(lt_shear)) # plotting the line tension # norm of the line tension vector - line_tension_norm = mean_stress_vector_norm(lt, borders, norm_level="points", vtype="t_vecs", - exclude_colony_edge=True) + line_tension_norm = mean_stress_vector_norm( + lt, borders, norm_level="points", vtype="t_vecs", exclude_colony_edge=True + ) res_dict["average magnitude line tension"] = line_tension_norm[1] res_dict["std magnitude line tension"] = line_tension_norm[2] # normal component of the line tension vector - line_tension_n = mean_stress_vector_norm(lt, borders, norm_level="points", vtype="t_normal", - exclude_colony_edge=True) + line_tension_n = mean_stress_vector_norm( + lt, borders, norm_level="points", vtype="t_normal", exclude_colony_edge=True + ) res_dict["average normal line tension"] = line_tension_n[1] res_dict["std normal line tension"] = line_tension_n[2] # shear component of the line tension vector - line_tension_sh = mean_stress_vector_norm(lt, borders, norm_level="points", vtype="t_shear", - exclude_colony_edge=True) + line_tension_sh = mean_stress_vector_norm( + lt, borders, norm_level="points", vtype="t_shear", exclude_colony_edge=True + ) res_dict["average shear line tension"] = line_tension_sh[1] res_dict["std shear line tension"] = line_tension_sh[2] - avg_cell_force = mean_stress_vector_norm(lt, borders, norm_level="cells", vtype="t_vecs", - exclude_colony_edge=True) + avg_cell_force = mean_stress_vector_norm( + lt, borders, norm_level="cells", vtype="t_vecs", exclude_colony_edge=True + ) res_dict["average cell force"] = avg_cell_force[1] res_dict["std cell force"] = avg_cell_force[2] - avg_cell_pressure = mean_stress_vector_norm(lt, borders, norm_level="cells", vtype="t_normal", - exclude_colony_edge=True) + avg_cell_pressure = mean_stress_vector_norm( + lt, borders, norm_level="cells", vtype="t_normal", exclude_colony_edge=True + ) res_dict["average cell pressure"] = avg_cell_pressure[1] res_dict["std cell pressure"] = avg_cell_pressure[2] - avg_cell_shear = mean_stress_vector_norm(lt, borders, norm_level="cells", vtype="t_shear", - exclude_colony_edge=True) + avg_cell_shear = mean_stress_vector_norm( + lt, borders, norm_level="cells", vtype="t_shear", exclude_colony_edge=True + ) res_dict["average cell shear"] = avg_cell_shear[1] res_dict["std cell shear"] = avg_cell_shear[2] attributes = dict( - fx=fx, fy=fy, ms=ms, cv=cv, borders_edge_lines=borders.edge_lines, borders_inter_shape=borders.inter_shape, - lt=lt, min_v=min_v, max_v=max_v + fx=fx, + fy=fy, + ms=ms, + cv=cv, + borders_edge_lines=borders.edge_lines, + borders_inter_shape=borders.inter_shape, + lt=lt, + min_v=min_v, + max_v=max_v, ) return res_dict, attributes diff --git a/saenopy/pyTFM/correct_stage_drift.py b/saenopy/pyTFM/correct_stage_drift.py index dd8d15e..6066365 100644 --- a/saenopy/pyTFM/correct_stage_drift.py +++ b/saenopy/pyTFM/correct_stage_drift.py @@ -14,17 +14,17 @@ def normalizing(img): def cropping_after_shift(image, shift_x, shift_y): if shift_x <= 0: - image = image[:, int(np.ceil(-shift_x)):] + image = image[:, int(np.ceil(-shift_x)) :] else: - image = image[:, :-int(np.ceil(shift_x))] + image = image[:, : -int(np.ceil(shift_x))] if shift_y <= 0: - image = image[int(np.ceil(-shift_y)):, :] + image = image[int(np.ceil(-shift_y)) :, :] else: - image = image[:-int(np.ceil(shift_y)), :] + image = image[: -int(np.ceil(shift_y)), :] return np.array(image, dtype=float) -def correct_stage_drift(image1, image2, additional_images=[]): +def correct_stage_drift(image1, image2, additional_images=None): """ # correcting frame shift between images of beads before and after cell removal. @@ -39,6 +39,8 @@ def correct_stage_drift(image1, image2, additional_images=[]): :param additional_images: :return: """ + if additional_images is None: + additional_images = [] # find shift with image registration shift_values = phase_cross_correlation(image1, image2, upsample_factor=100) @@ -59,7 +61,9 @@ def correct_stage_drift(image1, image2, additional_images=[]): additional_images_save = [] for add_image in additional_images: add_image_shift = shift(add_image, shift=(-shift_y, -shift_x), order=5) - add_image_norm = normalizing(cropping_after_shift(add_image_shift, shift_x, shift_y)) + add_image_norm = normalizing( + cropping_after_shift(add_image_shift, shift_x, shift_y) + ) add_image_save = Image.fromarray(add_image_norm * 255) additional_images_save.append(add_image_save) diff --git a/saenopy/pyTFM/graph_theory_for_cell_boundaries.py b/saenopy/pyTFM/graph_theory_for_cell_boundaries.py index f2e14bd..51e2db7 100755 --- a/saenopy/pyTFM/graph_theory_for_cell_boundaries.py +++ b/saenopy/pyTFM/graph_theory_for_cell_boundaries.py @@ -7,8 +7,8 @@ import numpy as np from scipy.spatial import cKDTree -inf = float('inf') -Edge = namedtuple('Edge', 'start, end, cost') +inf = float("inf") +Edge = namedtuple("Edge", "start, end, cost") class FindingBorderError(Exception): @@ -17,7 +17,9 @@ class FindingBorderError(Exception): def graph_to_mask(graph, points, dims): m = np.zeros(dims) - ps = np.array([y for x in list(graph.values()) for y in x]) # flattening list of point ids + ps = np.array( + [y for x in list(graph.values()) for y in x] + ) # flattening list of point ids ps_coord = points[ps] # getting coordinates m[ps_coord[:, 0], ps_coord[:, 1]] = 1 # writing points return m @@ -36,20 +38,25 @@ def check_connectivity(graph, ep): """ checking if removing a node from a graph changes the connectivity of the neighbouring nodes. in other words: are the neighbouring nodes still connected to one another even if I remove the original node. - The neighouring points must be connected via maximally one other node. (more would represent ana actual hole in the - sceletonized mask) + The neighbouring points must be connected via maximally one other node. (more would represent ana actual hole in the + skeletonized mask) This classifies loose ends and points that can be removed. :return: """ - # this works for points 1,2 or nneighbours + # this works for points 1,2 or neighbours # identifying an endpoint by checking if removing this point changes the connectivity of its neighbouring nodes - # i.e. if i break a connection between two points by removing ep + # i.e. if I break a connection between two points by removing ep l1_ps = graph[ep] # first layer of points # check if this causes a line break - l2_ps = [[pi for pi in graph[p] if pi != ep] for p in l1_ps] # second layer of points if original point was removed + l2_ps = [ + [pi for pi in graph[p] if pi != ep] for p in l1_ps + ] # second layer of points if original point was removed # third layer of points // don't need to go deeper due to properties of skeletonize - # also adding points form secnod layer - l3_ps = [np.unique(list(chain.from_iterable([graph[p] for p in l2_p] + [l2_p]))) for l2_p in l2_ps] + # also adding points form second layer + l3_ps = [ + np.unique(list(chain.from_iterable([graph[p] for p in l2_p] + [l2_p]))) + for l2_p in l2_ps + ] # check if all points in l1_ps are connected even if ep is removed connectivity = all([all([p in sub_group for p in l1_ps]) for sub_group in l3_ps]) # check for connection between points in layer 1--> no connection means @@ -57,19 +64,25 @@ def check_connectivity(graph, ep): return connectivity -def remove_endpoints(graph, ep, removed=[]): +def remove_endpoints(graph, ep, removed=None): """ recursive function to remove dead ends in a graph starting from point ep. Ep has one neighbour. Function stops if it hits a point with 3 neighbours or the removal of a point would cause the appearance of two more loose lines. :param graph: graph as a dictionary :param ep: start point + :param removed: :return: """ + if removed is None: + removed = [] + connectivity = check_connectivity(graph, ep) if connectivity: nps = graph[ep] - remove_point_from_graph(graph, ep) # removes the point and all connections form the graph + remove_point_from_graph( + graph, ep + ) # removes the point and all connections form the graph removed.append(ep) else: return @@ -86,7 +99,7 @@ def remove_point_from_graph(graph, point): :param point: :return: """ - nps = graph[point] # neighbouring points/nddes + nps = graph[point] # neighbouring points/nodes graph.pop(point) # removing the node of the graph for p in nps: # removing all connections to this node graph[p].remove(point) @@ -114,7 +127,9 @@ def find_dead_end_lines(graph, non_dead_end_points, max_id): dead_end_lines = {} for ep in eps_id: lps = find_path(graph, ep, non_dead_end_points, path=[]) - non_dead_end_points.extend(lps) # adding points in the newly discoverd line segments to termination points + non_dead_end_points.extend( + lps + ) # adding points in the newly discovered line segments to termination points if len(lps) > 3: # filtering single points and very small bits max_id += 1 dead_end_lines[max_id] = lps @@ -127,21 +142,27 @@ def find_lines_simple(graph): lines_points = {} i = 0 while len(graph_cp.keys()) > 0: - # first ednpoint, if no endpoint the first point - new_endpoint = next((x for x in iter(graph_cp.keys()) if len(graph_cp[x]) == 1), next(iter(graph_cp.keys()))) - line = find_path_to_endpoint(graph_cp, new_endpoint, path=[], first=True) # only explores one direction + # first endpoint, if no endpoint the first point + new_endpoint = next( + (x for x in iter(graph_cp.keys()) if len(graph_cp[x]) == 1), + next(iter(graph_cp.keys())), + ) + line = find_path_to_endpoint( + graph_cp, new_endpoint, path=[], first=True + ) # only explores one direction for p in line: remove_point_from_graph(graph_cp, p) if len(line) > 2: lines_points[i] = line i += 1 if i > 10000: - raise FindingBorderError("found more than 100000 lines; something went wrong") + raise FindingBorderError( + "found more than 100000 lines; something went wrong" + ) return lines_points -def find_path(graph, start, end, path=[]): - +def find_path(graph, start, end, path=None): """ recursive function finds a path (not necessarily the shortest one) through a graph from start to an end node (not necessarily the @@ -153,6 +174,8 @@ def find_path(graph, start, end, path=[]): :param path: list, all nodes visited on the way from start to the first endpoint :return: """ + if path is None: + path = [] path = path + [start] if start in end: return path @@ -160,14 +183,13 @@ def find_path(graph, start, end, path=[]): return None for node in graph[start]: if node not in path: - newpath = find_path(graph, node, end, path) - if newpath: - return newpath + new_path = find_path(graph, node, end, path) + if new_path: + return new_path return None # only partial path -def find_path_to_endpoint(graph, start, path=[], first=False): - +def find_path_to_endpoint(graph, start, path=None, first=False): """ recursive function finds a path to a (not specific) point with only one neighbour @@ -175,26 +197,31 @@ def find_path_to_endpoint(graph, start, path=[], first=False): :param graph: dict, graph :param start: int, start point, must be a key in the graph :param path: list, all nodes visited on the way from start to the first endpoint + :param first: :return: """ + if path is None: + path = [] path = path + [start] - if len(graph[start]) < 2 and not first: # stop if we reached a point with only one neighbour + if ( + len(graph[start]) < 2 and not first + ): # stop if we reached a point with only one neighbour return path if start not in graph.keys(): return None for node in graph[start]: if node not in path: - newpath = find_path_to_endpoint(graph, node, path) - if newpath: - return newpath + new_path = find_path_to_endpoint(graph, node, path) + if new_path: + return new_path return path # only partial path -def find_line_segement_recursive(graph, start, path=[], left_right=0): +def find_line_segment_recursive(graph, start, path=None, left_right=0): """ ---> would sometimes cause stack overflow/recursion error recursive function - finds path from a point going from the first or second neigbour until + finds path from a point going from the first or second neighbour until it reaches an intersection (point with three neighbours) :param graph: graph as a dictionary :param start: start point @@ -202,28 +229,34 @@ def find_line_segement_recursive(graph, start, path=[], left_right=0): :param left_right: define which neighbour from start to explore (0 or 1 :return: """ - if len(graph[start]) > 2: # stop if intersection (point with 3 neighbours is reached + if path is None: + path = [] + if ( + len(graph[start]) > 2 + ): # stop if intersection (point with 3 neighbours is reached return path # returns the function before next recursion - # other wise there is some overlap + # otherwise there is some overlap path = path + [start] new_ps = np.array(graph[start]) # next points if len(path) == 1: new_p = new_ps[left_right] # just choose one point else: - new_p = new_ps[new_ps != path[-2]][0] # next point that wasn't the previous point + new_p = new_ps[new_ps != path[-2]][ + 0 + ] # next point that wasn't the previous point # recursive function - newpath = find_line_segement_recursive(graph, new_p, path) # next step + new_path = find_line_segment_recursive(graph, new_p, path) # next step - if newpath: - return newpath # return if recursion is completed + if new_path: + return new_path # return if recursion is completed -def find_line_segement(graph, start, path=None, left_right=0): +def find_line_segment(graph, start, path=None, left_right=0): """ ---> would sometimes cause stack overflow/recursion error recursive function - finds path from a point going from the first or second neigbour until + finds path from a point going from the first or second neighbour until it reaches an intersection (point with three neighbours) :param graph: graph as a dictionary :param start: start point @@ -233,8 +266,7 @@ def find_line_segement(graph, start, path=None, left_right=0): """ # first point - path = [] - path.append(start) + path = [start] new_ps = graph[start] new_p = new_ps[left_right] # just choose one point to the left or right # break of if we already hit another intersection @@ -248,7 +280,7 @@ def find_line_segement(graph, start, path=None, left_right=0): new_p = [p for p in new_ps if p not in path] new_p = new_p[0] # check if the point has more (or less) then 2 neighbours. - # should be more then 2 to indicate intersection + # should be more than 2 to indicate intersection if len(graph[new_p]) != 2: break else: @@ -266,7 +298,7 @@ def mask_to_graph(mask, d=np.sqrt(2)): """ graph = defaultdict(list) points = np.array(np.where(mask)).T - point_tree = cKDTree(points) # look up table for nearest neigbours ?? + point_tree = cKDTree(points) # look up table for nearest neighbours ?? for i, p in enumerate(points): neighbours = point_tree.query_ball_point(p, d) neighbours.remove(i) # removing the point itself from list of its neighbours @@ -276,37 +308,43 @@ def mask_to_graph(mask, d=np.sqrt(2)): def identify_line_segments(graph, points): # """ - function to identify all line segments (representing individual cell boundaries. Segments are returned as a dictionary - with an id as key and a list of points (referring to the points array) that are included in the line. The - points are in correct order already + function to identify all line segments (representing individual cell boundaries). Segments are returned as a + dictionary with an id as key and a list of points (referring to the points array) that are included in the line. + The points are in correct order already :param graph: :param points: - :return: dictionary with orderd points in the line + :return: dictionary with ordered points in the line """ lines_dict = {} n = 0 # counter in while loop all_points = list(graph.keys()) # all points in the graph - intersect_ps = [key for key, values in graph.items() if len(values) > 2] # finding intersection points + intersect_ps = [ + key for key, values in graph.items() if len(values) > 2 + ] # finding intersection points if len(intersect_ps) == 0: raise FindingBorderError("Can't identify internal cell borders.") remaining = set(all_points) - set(intersect_ps) # remaining point ids - while len(remaining) > 0: # stop when all points are assigned to intersect or line segment + while ( + len(remaining) > 0 + ): # stop when all points are assigned to intersect or line segment start = next(iter(remaining)) # first point of remaining points # finding a line segment - line_seg = list(reversed(find_line_segement(graph, start=start, left_right=0)))[:-1] + find_line_segement(graph, - start=start, - left_right=1) + line_seg = list(reversed(find_line_segment(graph, start=start, left_right=0)))[ + :-1 + ] + find_line_segment(graph, start=start, left_right=1) remaining -= set(line_seg) # updating remaining list - # avoid single point lines, which are not usefull (cant really get normal vectors from them) + # avoid single point lines, which are not useful (cant really get normal vectors from them) if len(line_seg) > 1: lines_dict[n] = line_seg n = n + 1 if n > 20000: # expectation if loop should get suck - raise FindingBorderError("found more than 20000 cell borders; something went wrong") + raise FindingBorderError( + "found more than 20000 cell borders; something went wrong" + ) # plot to confirm correct lines # plt.figure() @@ -319,7 +357,9 @@ def identify_line_segments(graph, points): # return lines_dict -def find_neighbor_lines(graph, start_ps, other_endpoint, own_points, end_points, visited=[], neighbours=[]): +def find_neighbor_lines( + graph, start_ps, other_endpoint, own_points, end_points, visited=None, neighbours=None +): """ recursive function to find neighbouring line. Explores the graph around the endpoint of a line. Notes the id of the line if it hits another line. Doesn't explore any points beyond the endpoints of lines. @@ -335,9 +375,13 @@ def find_neighbor_lines(graph, start_ps, other_endpoint, own_points, end_points, :return: visited: list of visited nodes :return: neighbours """ + if neighbours is None: + neighbours = [] + if visited is None: + visited = [] visited = visited + start_ps # update visited list ## start must already be a list - next_ps = [graph[s] for s in start_ps] ## next points + next_ps = [graph[s] for s in start_ps] # next points next_ps = [p for ps in next_ps for p in ps] # flatten next_ps = list(np.unique(next_ps)) # removing duplication @@ -347,12 +391,16 @@ def find_neighbor_lines(graph, start_ps, other_endpoint, own_points, end_points, next_ps.remove(other_endpoint) # remove if point is in visited list or in own line - for p in copy.deepcopy(next_ps): ##### change in the list while iterating is not a nice idea--> + for p in copy.deepcopy( + next_ps + ): # change in the list while iterating is not a nice idea--> if p in visited or p in own_points: next_ps.remove(p) # extract if point can be found in other line - for p in copy.deepcopy(next_ps): ##### change in the list while iterating is not a nice idea--> make a copy + for p in copy.deepcopy( + next_ps + ): # change in the list while iterating is not a nice idea--> make a copy if p in end_points: next_ps.remove(p) neighbours.append(p) @@ -361,8 +409,15 @@ def find_neighbor_lines(graph, start_ps, other_endpoint, own_points, end_points, if len(next_ps) == 0: # stop recursion if no more next points are left return visited, neighbours - visited, neighbours = find_neighbor_lines(graph, next_ps, other_endpoint, own_points, end_points, visited=visited, - neighbours=neighbours) + visited, neighbours = find_neighbor_lines( + graph, + next_ps, + other_endpoint, + own_points, + end_points, + visited=visited, + neighbours=neighbours, + ) # return when iteration is finished if visited: return visited, neighbours @@ -371,30 +426,40 @@ def find_neighbor_lines(graph, start_ps, other_endpoint, own_points, end_points, def find_exact_line_endpoints(lines_points, points, graph): """ function to find the exact meeting points of lines. - First find the next closes points on neighbouring lines by exploring the graph. Then calcualtes a - new endpoint as center of mass of these neighbouring points. Results are stored in a seperate dictionaryl to be - used in spline interpoltaion + First find the next closes points on neighbouring lines by exploring the graph. Then calculates a + new endpoint as center of mass of these neighbouring points. Results are stored in a separate dictionary to be + used in spline interpolation :param lines_points: dictionary with line_id: list of all points in correct order :param points: array of point coordinates :param graph: dictionary with connectivity of points :return: lines_endpoints_com: dictionary with the line_id:[new endpoint at start, new_endpoint at end] """ - end_points = [[ps[0], ps[-1]] for ps in lines_points.values()] # all end points in lines + end_points = [ + [ps[0], ps[-1]] for ps in lines_points.values() + ] # all end points in lines end_points = [p for ps in end_points for p in ps] # finding all neighbouring edpoints for one endpoint of a line lines_endpoints = {} for line, l_points in lines_points.items(): - # points on the line with out both endpoints, other wise the algorithm can connect two endpoints on the same line + # points on the line without both endpoints, + # otherwise the algorithm can connect two endpoints on the same line l_points_core = l_points[1:-1] end1 = l_points[0] end2 = l_points[-1] - v, neighbours1 = find_neighbor_lines(graph, [end1], end2, l_points_core, end_points, visited=[], neighbours=[]) - v, neighbours2 = find_neighbor_lines(graph, [end2], end1, l_points_core, end_points, visited=[], neighbours=[]) - lines_endpoints[line] = (neighbours1 + [end1], neighbours2 + [end2]) # also adding own endpoints here + v, neighbours1 = find_neighbor_lines( + graph, [end1], end2, l_points_core, end_points, visited=[], neighbours=[] + ) + v, neighbours2 = find_neighbor_lines( + graph, [end2], end1, l_points_core, end_points, visited=[], neighbours=[] + ) + lines_endpoints[line] = ( + neighbours1 + [end1], + neighbours2 + [end2], + ) # also adding own endpoints here # note adding endpoints after find_neighbour_lines is easiest - ## calcualte new endpoints: + # calculate new endpoints: # just center of mass of the endpoints # write to new dictionary and use this in splines calculation, without any points in between lines_endpoints_com = {} diff --git a/saenopy/pyTFM/grid_setup_solids_py.py b/saenopy/pyTFM/grid_setup_solids_py.py index 0d0866f..2c3b9a5 100644 --- a/saenopy/pyTFM/grid_setup_solids_py.py +++ b/saenopy/pyTFM/grid_setup_solids_py.py @@ -12,37 +12,46 @@ import solidspy.postprocesor as pos import solidspy.solutil as sol from .stress_functions import calculate_stress_tensor, normal_vectors_from_splines -from .graph_theory_for_cell_boundaries import (mask_to_graph, FindingBorderError, find_lines_simple, - remove_endpoints_wrapper, graph_to_mask, identify_line_segments, - find_dead_end_lines, find_exact_line_endpoints) +from .graph_theory_for_cell_boundaries import ( + mask_to_graph, + FindingBorderError, + find_lines_simple, + remove_endpoints_wrapper, + graph_to_mask, + identify_line_segments, + find_dead_end_lines, + find_exact_line_endpoints, +) from .utilities_TFM import make_random_discrete_color_range, join_dictionary from scipy.interpolate import splprep, splev -from scipy.ndimage import binary_fill_holes +from scipy.ndimage import binary_fill_holes, binary_closing from scipy.optimize import least_squares from scipy.sparse import csr_matrix from scipy.sparse.linalg import lsqr from skimage.measure import regionprops -from skimage.morphology import skeletonize, remove_small_holes, remove_small_objects, label, binary_dilation - - -def show_points(ps, mask): - plt.figure() - plt.imshow(mask) - plt.plot(ps[:, 1], ps[:, 0], "or") +from skimage.morphology import ( + skeletonize, + remove_small_holes, + remove_small_objects, + label, + binary_dilation, +) def identify_cells(mask_area, mask_boundaries, points): """ function to identify cells. Each cell is a dictionary entry with a list of ids, referring to points. - :param mask: - :param area: + :param mask_area: :param mask_boundaries: + :param points: :return: """ cells = {} # dictionary containg a list of point idsthat sourround each cell - cells_area = {} # dictionary containg a all pixels belonging to that cell as boolean aray + cells_area = ( + {} + ) # dictionary containg a all pixels belonging to that cell as boolean aray # points_to_flatt array map: # labeling each cell m = mask_area.astype(int) - mask_boundaries.astype(int) @@ -54,13 +63,19 @@ def identify_cells(mask_area, mask_boundaries, points): points_fl = points_fl[sort_ids] # for i, l in enumerate( - np.unique(ml)[1:]): # getting each cell border by binary dilation of one pixel; iterating over each cell + np.unique(ml)[1:] + ): # getting each cell border by binary dilation of one pixel; iterating over each cell m_part = (ml == l).astype(bool) # extracting a cell area - edge = np.logical_and(binary_dilation(m_part), ~m_part) # getting the boundary of a cell + edge = np.logical_and( + binary_dilation(m_part), ~m_part + ) # getting the boundary of a cell ps = np.array(np.where(edge)).T # finding coordinates - ps_fl = (ps[:, 0]) + mask_area.shape[0] * (ps[:, 1]) # convert coordinates to the one of a flat array - p_ids = sort_ids[np.searchsorted(points_fl, - ps_fl)] # find indices, where i would need to insert,supposed to be the fastest way + ps_fl = (ps[:, 0]) + mask_area.shape[0] * ( + ps[:, 1] + ) # convert coordinates to the one of a flat array + p_ids = sort_ids[ + np.searchsorted(points_fl, ps_fl) + ] # find indices, where i would need to insert,supposed to be the fastest way # and read index from unsorted list cells[i] = p_ids # save to dictionary cells_area[i] = m_part @@ -105,7 +120,9 @@ def spline_interpolation(line, points, k=3, endpoints=None): x = np.concatenate([x, [endpoints[1][1]]], axis=0) y = np.concatenate([y, [endpoints[1][0]]], axis=0) - k = len(x) - 1 if len(x) <= 3 else 3 # addapt spline order, according to number of points + k = ( + len(x) - 1 if len(x) <= 3 else 3 + ) # addapt spline order, according to number of points tck, u = splprep([x, y], s=10, k=k) ### parametric spline interpolation # fits essentially function: [x,y] =f(t) , t(paramter) is default np.linspace(0,1,len(x)) # tck is array with x_position of knot, y_position of knot, parameters of the plne, order of the spline @@ -139,9 +156,13 @@ def arrange_lines_from_endpoints(cells_lines, lines_endpoints_com): for cell_id, line_ids in cells_lines.items(): # extracting relevant endpoints - local_endpoints = {line_id: lines_endpoints_com[line_id] for line_id in line_ids} + local_endpoints = { + line_id: lines_endpoints_com[line_id] for line_id in line_ids + } # rearranging endpoints into an suitabel array: axis0: lines axis 1:[endpoint1, endpoint2], axis3: x,y coordinates - eps = np.array([np.array([value[0], value[1]]) for value in local_endpoints.values()]) + eps = np.array( + [np.array([value[0], value[1]]) for value in local_endpoints.values()] + ) new_line_ids = [] # newly_arranged line_ids p_ind1 = 0 # start ids @@ -152,10 +173,14 @@ def arrange_lines_from_endpoints(cells_lines, lines_endpoints_com): eps[p_ind1, p_ind2] = np.nan # remove the point from array # find second occurrence, by taking the norm of the diffrence between the current point and all other points # this should be zero - np_ind1, np_ind2 = np.array(np.where(np.linalg.norm(eps - point, axis=2) == 0)).T[0] + np_ind1, np_ind2 = np.array( + np.where(np.linalg.norm(eps - point, axis=2) == 0) + ).T[0] new_line_ids.append(line_ids[np_ind1]) # note corresponding line_ids eps[np_ind1, np_ind2] = np.nan # remove this point from array - p_ind1, p_ind2 = np_ind1, np.abs(np_ind2 - 1) # change to the other end of the line + p_ind1, p_ind2 = np_ind1, np.abs( + np_ind2 - 1 + ) # change to the other end of the line cells_lines_new[cell_id] = new_line_ids # update dictionary @@ -169,10 +194,13 @@ def find_edge_lines(cells_lines): :param cells_lines: dictionary with cell_id:[associated line_ids] :return: edge_lines: lsit of line ids at the edge of the cell colony """ - all_lines = np.array(list(chain.from_iterable(cells_lines.values()))) # unpacking all line ids + all_lines = np.array( + list(chain.from_iterable(cells_lines.values())) + ) # unpacking all line ids counts = Counter(all_lines) # counting occurences - edge_lines = [line for line, count in counts.items() if - count == 1] # select if line_id was only associated to one cell + edge_lines = [ + line for line, count in counts.items() if count == 1 + ] # select if line_id was only associated to one cell return edge_lines @@ -190,7 +218,9 @@ def center_of_mass_cells(cells_points, points): return cells_com -def remove_circular_line(allLines_points, lines_endpoints_com, lines_points, lines_endpoints): +def remove_circular_line( + allLines_points, lines_endpoints_com, lines_points, lines_endpoints +): """ finds lines that are circular by checking if the first and second endpoint are identical. The lines are delted from all input dictionaries @@ -200,8 +230,11 @@ def remove_circular_line(allLines_points, lines_endpoints_com, lines_points, lin :return: """ # finding all lines where first and second endpoint is identical - circular = [l_id for l_id, endpoints in lines_endpoints_com.items() if - np.linalg.norm(endpoints[0] - endpoints[1]) == 0] + circular = [ + l_id + for l_id, endpoints in lines_endpoints_com.items() + if np.linalg.norm(endpoints[0] - endpoints[1]) == 0 + ] # print(circular) # clearing these lines from the input dictionaries for l_id in circular: @@ -221,8 +254,10 @@ def interpolate_cell_area(cells_area, shape): def interpolate_points_dict(points_dict, shape_target, shape_orgin): points_dict_interp = {} for p_id, coords in points_dict.items(): - points_dict_interp[p_id] = (interpolation_single_point(coords[0], shape_target, shape_orgin), - interpolation_single_point(coords[1], shape_target, shape_orgin)) + points_dict_interp[p_id] = ( + interpolation_single_point(coords[0], shape_target, shape_orgin), + interpolation_single_point(coords[1], shape_target, shape_orgin), + ) return points_dict_interp @@ -237,12 +272,18 @@ def __init__(self, mask_boundaries, shape, graph, points): # any point id is the index in the points array (contains coordinate of these points self.graph_wp = graph self.points_wp = points - self.graph, self.points, removed = remove_endpoints_wrapper(self.graph_wp, self.points_wp) + self.graph, self.points, removed = remove_endpoints_wrapper( + self.graph_wp, self.points_wp + ) # masks, graph and points excluding dead-end lines - self.mask_boundaries = graph_to_mask(self.graph, self.points, mask_boundaries.shape) # rebuilding the mask + self.mask_boundaries = graph_to_mask( + self.graph, self.points, mask_boundaries.shape + ) # rebuilding the mask # interpolate points to the size of the future FEM_grid - self.points_interpol = interpolation_single_point(self.points, self.inter_shape, self.mask_boundaries.shape) + self.points_interpol = interpolation_single_point( + self.points, self.inter_shape, self.mask_boundaries.shape + ) # interpolation factors used in the fun # self.inerpol_factors=np.array([self.inter_shape[0] /self.mask_boundaries.shape[0], self.inter_shape[1] / self.mask_boundaries.shape[1]]) # points as dictionary with key=points id, values: points coordinates @@ -252,31 +293,43 @@ def __init__(self, mask_boundaries, shape, graph, points): self.lines_points = identify_line_segments(self.graph, self.points_interpol) # cells as a dictionary with key=cell id, values: ids of containing points (not ordered) - self.cells_points, self.cells_area = identify_cells(self.mask_boundaries, - binary_fill_holes(self.mask_boundaries), self.points) + self.cells_points, self.cells_area = identify_cells( + self.mask_boundaries, binary_fill_holes(self.mask_boundaries), self.points + ) # interpolate the area of individual cells to the size of deformation - self.cells_area_interpol = interpolate_cell_area(self.cells_area, self.inter_shape) + self.cells_area_interpol = interpolate_cell_area( + self.cells_area, self.inter_shape + ) # self.points_lines = invert_dictionary(self.lines_points) # point_id:line_id self.max_line_id = np.max(list(self.lines_points.keys())) - self.de_lines_points, self.max_line_id = find_dead_end_lines(self.graph_wp, list(self.graph.keys()), - self.max_line_id) - self.de_endpoints = {key: (self.points[value[0]], self.points[value[-1]]) for key, value in - self.de_lines_points.items()} # using exact endpoints for the dead end lines + self.de_lines_points, self.max_line_id = find_dead_end_lines( + self.graph_wp, list(self.graph.keys()), self.max_line_id + ) + self.de_endpoints = { + key: (self.points[value[0]], self.points[value[-1]]) + for key, value in self.de_lines_points.items() + } # using exact endpoints for the dead end lines self.allLines_points = join_dictionary(self.lines_points, self.de_lines_points) # dictionary with endpoints, needed to completely fill the gaps between all cell_lines - self.lines_endpoints_com, self.lines_endpoints = find_exact_line_endpoints(self.lines_points, self.points, - self.graph) + self.lines_endpoints_com, self.lines_endpoints = find_exact_line_endpoints( + self.lines_points, self.points, self.graph + ) - #self.simple_line_plotting(self.allLines_points, subset=np.inf) - #for p in self.lines_endpoints_com.values(): + # self.simple_line_plotting(self.allLines_points, subset=np.inf) + # for p in self.lines_endpoints_com.values(): # plt.plot(p[0][1],p[0][0],"o") # plt.plot(p[1][1], p[1][0], "o") # removing all lines that are predicted to be circular. Mostly a problem for very short lines - remove_circular_line(self.allLines_points, self.lines_endpoints_com, self.lines_points, self.lines_endpoints) + remove_circular_line( + self.allLines_points, + self.lines_endpoints_com, + self.lines_points, + self.lines_endpoints, + ) # center of mass of cells, calculated only from the hull points self.cells_com = center_of_mass_cells(self.cells_points, self.points) @@ -287,18 +340,25 @@ def __init__(self, mask_boundaries, shape, graph, points): self.lines_cells = defaultdict(list) for l_id, l in self.lines_points.items(): for c_id, c in self.cells_points.items(): - if l[int(len(l)/2)] in c: + if l[int(len(l) / 2)] in c: self.cells_lines[c_id].append(l_id) self.lines_cells[l_id].append(c_id) # using the new endpoints to arrange lines in the correct way - self.cells_lines = arrange_lines_from_endpoints(self.cells_lines, self.lines_endpoints_com) + self.cells_lines = arrange_lines_from_endpoints( + self.cells_lines, self.lines_endpoints_com + ) # adding dead end endpoints only now to avoid complications when identifying cells - self.de_endpoints = {key: (self.points[value[0]], self.points[value[-1]]) for key, value in - self.de_lines_points.items()} - self.lines_endpoints_com = join_dictionary(self.lines_endpoints_com, self.de_endpoints) - self.lines_endpoints_interpol = interpolate_points_dict(self.lines_endpoints_com, self.inter_shape, - self.mask_boundaries.shape) + self.de_endpoints = { + key: (self.points[value[0]], self.points[value[-1]]) + for key, value in self.de_lines_points.items() + } + self.lines_endpoints_com = join_dictionary( + self.lines_endpoints_com, self.de_endpoints + ) + self.lines_endpoints_interpol = interpolate_points_dict( + self.lines_endpoints_com, self.inter_shape, self.mask_boundaries.shape + ) # list of ids self.point_ids = list(self.points_dict.keys()) @@ -309,11 +369,16 @@ def __init__(self, mask_boundaries, shape, graph, points): # list of dead end lines self.dead_end_lines = list(self.de_lines_points.keys()) # list of central boundary (non dead-end, non-edge lines) - self.central_lines = [line_id for line_id in self.allLines_points.keys() if - line_id not in self.edge_lines and line_id not in self.dead_end_lines] + self.central_lines = [ + line_id + for line_id in self.allLines_points.keys() + if line_id not in self.edge_lines and line_id not in self.dead_end_lines + ] self.n_cells = len(self.cell_ids) # list with original line lengths--> later used for interpolation - self.line_lengths = {key: len(value) for key, value in self.allLines_points.items()} + self.line_lengths = { + key: len(value) for key, value in self.allLines_points.items() + } # dictionary containing the spline represetation of the points as a parametric function # [x,y]=f(u). u is always np.linspace(0,1,"number of points in the line). Use scipy.interpolate.splev @@ -330,14 +395,19 @@ def __init__(self, mask_boundaries, shape, graph, points): for line_id, line_ps in self.allLines_points.items(): endpoints = self.lines_endpoints_interpol[line_id] - k = len(line_ps) + 2 - 1 if len(line_ps) <= 3 else 3 # addapt spline order, according to number of points + k = ( + len(line_ps) + 2 - 1 if len(line_ps) <= 3 else 3 + ) # addapt spline order, according to number of points # spline order must be below nomber of points, so choose 2 if lne has one point + 2 endpoints - tck, u, points_new = spline_interpolation(line_ps, self.points_interpol, k=k, - endpoints=endpoints) # spline interpolation + tck, u, points_new = spline_interpolation( + line_ps, self.points_interpol, k=k, endpoints=endpoints + ) # spline interpolation self.lines_splines[line_id] = tck # saving the spline object # saving a few oints and n vectors for easy representations/ debugging, # these points will not be used further - n_vectors = normal_vectors_from_splines(u, tck) # calculating normal vectors as a list + n_vectors = normal_vectors_from_splines( + u, tck + ) # calculating normal vectors as a list self.lines_n_vectors[line_id] = n_vectors self.lines_spline_points[line_id] = points_new @@ -358,26 +428,42 @@ def cut_to_FEM_grid(self, FEM_mask): self.lines_endpoints_com.pop(l_id, None) self.de_lines_points.pop(l_id, None) self.lines_points.pop(l_id, None) - with suppress(ValueError): self.edge_lines.remove(l_id) - with suppress(ValueError): self.edge_lines.remove(l_id) - with suppress(ValueError): self.edge_lines.remove(l_id) + with suppress(ValueError): + self.edge_lines.remove(l_id) + with suppress(ValueError): + self.edge_lines.remove(l_id) + with suppress(ValueError): + self.edge_lines.remove(l_id) self.lines_outside.append(l_id) def filter_small_de_line(self, min_length): - for l_id in copy.deepcopy(self.dead_end_lines): # does not filter small line segments around cells - + for l_id in copy.deepcopy( + self.dead_end_lines + ): # does not filter small line segments around cells - if self.line_lengths[l_id] < min_length: - with suppress(AttributeError): self.allLines_points.pop(l_id, None) - with suppress(AttributeError): self.lines_endpoints_com.pop(l_id, None) - with suppress(AttributeError): self.lines_endpoints_interpol.pop(l_id, None) - with suppress(AttributeError): self.lines_splines.pop(l_id, None) - with suppress(AttributeError): self.lines_n_vectors.pop(l_id, None) - with suppress(AttributeError): self.lines_spline_points.pop(l_id, None) - with suppress(AttributeError): self.line_lengths.pop(l_id, None) - with suppress(AttributeError): self.de_endpoints.pop(l_id, None) - with suppress(AttributeError): self.de_lines_points.pop(l_id, None) - - with suppress(ValueError, AttributeError): self.line_ids.remove(l_id) # list - with suppress(ValueError, AttributeError): self.dead_end_lines.remove(l_id) + with suppress(AttributeError): + self.allLines_points.pop(l_id, None) + with suppress(AttributeError): + self.lines_endpoints_com.pop(l_id, None) + with suppress(AttributeError): + self.lines_endpoints_interpol.pop(l_id, None) + with suppress(AttributeError): + self.lines_splines.pop(l_id, None) + with suppress(AttributeError): + self.lines_n_vectors.pop(l_id, None) + with suppress(AttributeError): + self.lines_spline_points.pop(l_id, None) + with suppress(AttributeError): + self.line_lengths.pop(l_id, None) + with suppress(AttributeError): + self.de_endpoints.pop(l_id, None) + with suppress(AttributeError): + self.de_lines_points.pop(l_id, None) + + with suppress(ValueError, AttributeError): + self.line_ids.remove(l_id) # list + with suppress(ValueError, AttributeError): + self.dead_end_lines.remove(l_id) def return_n_array(self, fill_nan=True): """ @@ -385,10 +471,19 @@ def return_n_array(self, fill_nan=True): :return: """ if fill_nan: - n_array = np.zeros((self.mask_boundaries.shape[0], self.mask_boundaries.shape[1], 2)) + np.nan + n_array = ( + np.zeros( + (self.mask_boundaries.shape[0], self.mask_boundaries.shape[1], 2) + ) + + np.nan + ) else: - n_array = np.zeros((self.mask_boundaries.shape[0], self.mask_boundaries.shape[1], 2)) - for vecs, ps in zip(self.lines_n_vectors.values(), self.lines_spline_points.values()): + n_array = np.zeros( + (self.mask_boundaries.shape[0], self.mask_boundaries.shape[1], 2) + ) + for vecs, ps in zip( + self.lines_n_vectors.values(), self.lines_spline_points.values() + ): n_array[ps[:, 1], ps[:, 0]] = vecs return n_array @@ -415,21 +510,35 @@ def vizualize_lines_and_cells(self, sample_factor=1, plot_n_vectors=False): color = colors[np.array([l in line_ids for line_ids in line_classifier])][0] p_indices = np.array(list(range(len(ps)))) # all indicces # randomly choosing a few indices, without first and laast index - ps_select = np.random.choice(p_indices[1:-1], size=int((len(ps) - 2) * sample_factor), replace=False) - ps_select = np.append(ps_select, p_indices[np.array([0, -1])]) # adding back first and last index + ps_select = np.random.choice( + p_indices[1:-1], size=int((len(ps) - 2) * sample_factor), replace=False + ) + ps_select = np.append( + ps_select, p_indices[np.array([0, -1])] + ) # adding back first and last index plt.plot(self.points[ps][:, 1], self.points[ps][:, 0], "o", color=color) for p in ps[ps_select]: # labeling selected points - plt.text(self.points[p][1] + 1 * offset * l, self.points[p][0] + 1 * offset * l, s=str(l), - color="green") + plt.text( + self.points[p][1] + 1 * offset * l, + self.points[p][0] + 1 * offset * l, + s=str(l), + color="green", + ) # plotting cel id at center of mass of cell for cell_id, com in self.cells_com.items(): plt.text(com[1], com[0], str(cell_id), color="red", fontsize=13) if plot_n_vectors: - for points, n_vectors in zip(self.lines_spline_points.values(), self.lines_n_vectors.values()): - for n_vec, p in zip(n_vectors, - interpolation_single_point(points, self.mask_boundaries.shape, self.inter_shape)): + for points, n_vectors in zip( + self.lines_spline_points.values(), self.lines_n_vectors.values() + ): + for n_vec, p in zip( + n_vectors, + interpolation_single_point( + points, self.mask_boundaries.shape, self.inter_shape + ), + ): plt.arrow(p[0], p[1], n_vec[0], n_vec[1], head_width=0.15) plt.legend() return fig @@ -442,8 +551,12 @@ def vizualize_splines(self, sample_factor=1, subset=np.inf): for i, (l_id, points) in enumerate(self.lines_spline_points.items()): if i > subset: break - points = interpolation_single_point(points, self.mask_boundaries.shape, self.inter_shape) - plt.plot(points[::sample_factor, 0], points[::sample_factor, 1], color=colors[i]) + points = interpolation_single_point( + points, self.mask_boundaries.shape, self.inter_shape + ) + plt.plot( + points[::sample_factor, 0], points[::sample_factor, 1], color=colors[i] + ) def simple_line_plotting(self, lines, subset=np.inf): plt.figure() @@ -471,7 +584,9 @@ def __init__(self, mask_boundaries, shape, graph, points): self.points = points # finding line segments self.lines_points = find_lines_simple(self.graph) - self.points_interpol = interpolation_single_point(self.points, self.inter_shape, self.mask_boundaries.shape) + self.points_interpol = interpolation_single_point( + self.points, self.inter_shape, self.mask_boundaries.shape + ) self.lines_splines = defaultdict(list) self.lines_n_vectors = defaultdict(list) @@ -480,11 +595,15 @@ def __init__(self, mask_boundaries, shape, graph, points): # spline interpolation for line_id, line_ps in self.lines_points.items(): # spline order must be below nomber of points, so choose 2 if lne has one point + 2 endpoints - tck, u, points_new = spline_interpolation(line_ps, self.points_interpol) # spline interpolation + tck, u, points_new = spline_interpolation( + line_ps, self.points_interpol + ) # spline interpolation self.lines_splines[line_id] = tck # saving the spline object # saving a few oints and n vectors for easy representations/ debugging, # these points will not be used further - n_vectors = normal_vectors_from_splines(u, tck) # calculating normal vectors as a list + n_vectors = normal_vectors_from_splines( + u, tck + ) # calculating normal vectors as a list self.lines_n_vectors[line_id] = n_vectors self.lines_spline_points[line_id] = points_new @@ -494,10 +613,14 @@ def __init__(self, mask_boundaries, shape, graph, points): self.edge_lines = [] self.dead_end_lines = [] self.central_lines = self.line_ids - self.line_lengths = {key: len(value) for key, value in self.lines_points.items()} + self.line_lengths = { + key: len(value) for key, value in self.lines_points.items() + } # very rough estimate of cell number - label_mask, self.n_cells = label(~mask_boundaries, connectivity=1, return_num=True) + label_mask, self.n_cells = label( + ~mask_boundaries, connectivity=1, return_num=True + ) def find_borders(mask, shape, raise_error=True, type="colony", min_length=0): @@ -546,7 +669,9 @@ def interpolation(mask, dims, min_cell_size=100, dtype=bool): coords[1] = coords[1] * interpol_factors[1] # interpolating xy coordinates coords = np.round(coords).astype(int) - coords[0, coords[0] >= dims[0]] = dims[0] - 1 # fixing issue when interpolated object is just at the image border + coords[0, coords[0] >= dims[0]] = ( + dims[0] - 1 + ) # fixing issue when interpolated object is just at the image border coords[1, coords[1] >= dims[1]] = dims[1] - 1 mask_int = np.zeros(dims) @@ -554,15 +679,19 @@ def interpolation(mask, dims, min_cell_size=100, dtype=bool): mask_int = mask_int.astype(int) # filling gaps if we interpolate upwards if dims[0] * dims[1] >= mask.shape[0] * mask.shape[1]: - iter = int(np.ceil(np.max([mask.shape[0] / dims[0], mask.shape[0] / dims[0]])) * 5) # times 5 is safety factor - mask_int = binary_clo(mask_int, iterations=10) + iter = int( + np.ceil(np.max([mask.shape[0] / dims[0], mask.shape[0] / dims[0]])) * 5 + ) # times 5 is safety factor + mask_int = binary_closing(mask_int, iterations=10) print(iter) return mask_int.astype(bool) def interpolation_single_point(point, shape_target, shape_origin): # is also works with 2d arrays of shape(n,2) - interpol_factors = np.array([shape_target[0] / shape_origin[0], shape_target[1] / shape_origin[1]]) + interpol_factors = np.array( + [shape_target[0] / shape_origin[0], shape_target[1] / shape_origin[1]] + ) point_interp = point * interpol_factors return point_interp @@ -570,25 +699,34 @@ def interpolation_single_point(point, shape_target, shape_origin): def cut_mask_from_edge(mask, cut_factor, warn_flag=False, fill=True): sum_mask1 = np.sum(mask) dims = mask.shape - inds = [int(dims[0] * cut_factor), int(dims[0] - (dims[0] * cut_factor)), int(dims[1] * cut_factor), - int(dims[1] - (dims[1] * cut_factor))] + inds = [ + int(dims[0] * cut_factor), + int(dims[0] - (dims[0] * cut_factor)), + int(dims[1] * cut_factor), + int(dims[1] - (dims[1] * cut_factor)), + ] if fill: # filling to the original shape mask_cut = copy.deepcopy(mask) - mask_cut[:inds[0], :] = 0 - mask_cut[inds[1]:, :] = 0 - mask_cut[:, :inds[2]] = 0 - mask_cut[:, inds[3]:] = 0 + mask_cut[: inds[0], :] = 0 + mask_cut[inds[1] :, :] = 0 + mask_cut[:, : inds[2]] = 0 + mask_cut[:, inds[3] :] = 0 else: # new array with new shape mask_cut = np.zeros((inds[1] - inds[0], inds[3] - inds[2])) - mask_cut = mask[inds[0]:inds[1], inds[2]:inds[3]] + mask_cut = mask[inds[0] : inds[1], inds[2] : inds[3]] sum_mask2 = np.sum(mask_cut) - warn = "mask was cut close to image edge" if (sum_mask2 < sum_mask1 and warn_flag) else "" + warn = ( + "mask was cut close to image edge" + if (sum_mask2 < sum_mask1 and warn_flag) + else "" + ) return mask_cut, warn def FEM_simulation(nodes, elements, loads, mats, mask_area, verbose=False, **kwargs): from packaging import version + if version.parse(solidspy.__version__) > version.parse("1.1.0"): cond = nodes[:, -2:] nodes = nodes[:, :-2] @@ -607,11 +745,15 @@ def FEM_simulation(nodes, elements, loads, mats, mask_area, verbose=False, **kwa # System assembly KG = ass.assembler(elements, mats, nodes, neq, DME, sparse=True) - if isinstance(KG, tuple) or version.parse(solidspy.__version__) > version.parse("1.1.0"): + if isinstance(KG, tuple) or version.parse(solidspy.__version__) > version.parse( + "1.1.0" + ): KG = KG[0] RHSG = ass.loadasem(loads, IBC, neq) - if np.sum(IBC == -1) < 3: # 1 or zero fixed nodes/ pure neumann-boundary-condition system needs further constraints + if ( + np.sum(IBC == -1) < 3 + ): # 1 or zero fixed nodes/ pure neumann-boundary-condition system needs further constraints # System solution with custom conditions # solver with constraints to zero translation and zero rotation UG_sol, rx = custom_solver(KG, RHSG, mask_area, nodes, IBC, verbose=verbose) @@ -624,8 +766,12 @@ def FEM_simulation(nodes, elements, loads, mats, mask_area, verbose=False, **kwa # average shear and normal stress on the colony area UC = pos.complete_disp(IBC, nodes, UG_sol) # uc are x and y displacements - E_nodes, S_nodes = pos.strain_nodes(nodes, elements, mats, UC) # stresses and strains - stress_tensor = calculate_stress_tensor(S_nodes, nodes, dims=mask_area.shape) # assembling the stress tensor + E_nodes, S_nodes = pos.strain_nodes( + nodes, elements, mats, UC + ) # stresses and strains + stress_tensor = calculate_stress_tensor( + S_nodes, nodes, dims=mask_area.shape + ) # assembling the stress tensor return UG_sol, stress_tensor @@ -642,7 +788,9 @@ def grid_setup(mask_area, f_x, f_y, E=1, sigma=0.5, edge_factor=0): :return: """ - coords = np.array(np.where(mask_area)) # retrieving all coordintates from the points in the mask + coords = np.array( + np.where(mask_area) + ) # retrieving all coordintates from the points in the mask # setting up nodes list:[node_id,x_coordinate,y_coordinate,fixation_y,fixation_x] nodes = np.zeros((coords.shape[1], 5)) @@ -652,13 +800,19 @@ def grid_setup(mask_area, f_x, f_y, E=1, sigma=0.5, edge_factor=0): # creating an 2D array, with the node id of each pixel. Non assigned pixel is -1. ids = np.zeros(mask_area.shape).T - 1 - ids[coords[0], coords[1]] = np.arange(coords.shape[1], dtype=int) # filling with node ids + ids[coords[0], coords[1]] = np.arange( + coords.shape[1], dtype=int + ) # filling with node ids # fix all nodes that are exactely at the edge of the image (minus any regions close to the image edge that are # supposed to be ignored)in the movement direction perpendicular to the edge ids_cut, w = cut_mask_from_edge(ids, edge_factor, "", fill=False) - edge_nodes_horizontal = np.hstack([ids_cut[:, 0], ids_cut[:, -1]]).astype(int) # upper and lower image edge - edge_nodes_vertical = np.hstack([ids_cut[0, :], ids_cut[-1, :]]).astype(int) # left and right image edge + edge_nodes_horizontal = np.hstack([ids_cut[:, 0], ids_cut[:, -1]]).astype( + int + ) # upper and lower image edge + edge_nodes_vertical = np.hstack([ids_cut[0, :], ids_cut[-1, :]]).astype( + int + ) # left and right image edge edge_nodes_horizontal = edge_nodes_horizontal[edge_nodes_horizontal >= 0] edge_nodes_vertical = edge_nodes_vertical[edge_nodes_vertical >= 0] nodes[edge_nodes_vertical, 3] = -1 # fixed in x direction @@ -672,14 +826,21 @@ def grid_setup(mask_area, f_x, f_y, E=1, sigma=0.5, edge_factor=0): # list the square(node,node left,node left down, node down) for each node. These are all posiible square shaped # elements, with the coorect orientation - sqr = [(coords[0], coords[1] - 1), (coords[0] - 1, coords[1] - 1), (coords[0] - 1, coords[1]), - (coords[0], coords[1])] + sqr = [ + (coords[0], coords[1] - 1), + (coords[0] - 1, coords[1] - 1), + (coords[0] - 1, coords[1]), + (coords[0], coords[1]), + ] # this produce negative indices, when at the edge of the mask # filtering these values - filter = np.sum(np.array([(s[0] < 0) + (s[1] < 0) for s in sqr]), - axis=0) > 0 # logical to find any square with negative coordinates + filter = ( + np.sum(np.array([(s[0] < 0) + (s[1] < 0) for s in sqr]), axis=0) > 0 + ) # logical to find any square with negative coordinates sqr = [(s[0][~filter], s[1][~filter]) for s in sqr] # applying filter - elements = elements[~filter] # shortening length of elements list according to the same filter + elements = elements[ + ~filter + ] # shortening length of elements list according to the same filter # enter node ids in elements, needs counter clockwise arangement # check by calling pyTFM.functions_for_cell_colonie.plot_grid(nodes,elements,inverted_axis=False,symbol_size=4,arrows=True,image=0) @@ -708,11 +869,13 @@ def grid_setup(mask_area, f_x, f_y, E=1, sigma=0.5, edge_factor=0): def prepare_forces(tx, ty, ps, mask): - f_x = tx * ((ps * (10 ** -6)) ** 2) # point force for each node from tractions - f_y = ty * ((ps * (10 ** -6)) ** 2) + f_x = tx * ((ps * (10**-6)) ** 2) # point force for each node from tractions + f_y = ty * ((ps * (10**-6)) ** 2) f_x[~mask] = np.nan # setting all values outside of mask area to zero f_y[~mask] = np.nan - f_x_c1 = f_x - np.nanmean(f_x) # normalizing traction force to sum up to zero (no displacement) + f_x_c1 = f_x - np.nanmean( + f_x + ) # normalizing traction force to sum up to zero (no displacement) f_y_c1 = f_y - np.nanmean(f_y) f_x_c2, f_y_c2, p = correct_torque(f_x_c1, f_y_c1, mask) return f_x_c2, f_y_c2 @@ -722,22 +885,31 @@ def correct_torque(fx, fy, mask_area): com = regionprops(mask_area.astype(int))[0].centroid # finding center of mass com = (com[1], com[0]) # as x y coordinate - c_x, c_y = np.meshgrid(range(fx.shape[1]), range(fx.shape[0])) # arrays with all x and y coordinates + c_x, c_y = np.meshgrid( + range(fx.shape[1]), range(fx.shape[0]) + ) # arrays with all x and y coordinates r = np.zeros((fx.shape[0], fx.shape[1], 2)) # array with all positional vectors r[:, :, 0] = c_x # note maybe its also enough to chose any point as refernece point r[:, :, 1] = c_y r = r - np.array(com) - f = np.zeros((fx.shape[0], fx.shape[1], 2), dtype="float64") # array with all force vectors + f = np.zeros( + (fx.shape[0], fx.shape[1], 2), dtype="float64" + ) # array with all force vectors f[:, :, 0] = fx f[:, :, 1] = fy - q = np.zeros((fx.shape[0], fx.shape[1], 2), dtype="float64") # rotated positional vectors + q = np.zeros( + (fx.shape[0], fx.shape[1], 2), dtype="float64" + ) # rotated positional vectors def get_torque_angle(p): - q[:, :, 0] = + np.cos(p) * (f[:, :, 0]) - np.sin(p) * (f[:, :, 1]) # whats the mathematics behind this?? - q[:, :, 1] = + np.sin(p) * (f[:, :, 0]) + np.cos(p) * (f[:, :, 1]) + q[:, :, 0] = +np.cos(p) * (f[:, :, 0]) - np.sin(p) * ( + f[:, :, 1] + ) # whats the mathematics behind this?? + q[:, :, 1] = +np.sin(p) * (f[:, :, 0]) + np.cos(p) * (f[:, :, 1]) torque = np.abs( - np.nansum(np.cross(r, q, axisa=2, axisb=2))) ## using nna sum to only look at force values in mask + np.nansum(np.cross(r, q, axisa=2, axisb=2)) + ) ## using nna sum to only look at force values in mask return torque.astype("float64") # plotting torque angle relation ship @@ -763,17 +935,33 @@ def get_torque_angle(p): # there seems to be a bug when using very small tolerances close to the machine precision limit (eps) # in rare cases there is an error. see also https://github.com/scipy/scipy/issues/11572 try: - p = least_squares(fun=get_torque_angle, x0=pstart, method="lm", - max_nfev=100000000, xtol=eps, ftol=eps, gtol=eps, args=())["x"] + p = least_squares( + fun=get_torque_angle, + x0=pstart, + method="lm", + max_nfev=100000000, + xtol=eps, + ftol=eps, + gtol=eps, + args=(), + )["x"] except KeyError: eps *= 5 - p = least_squares(fun=get_torque_angle, x0=pstart, method="lm", - max_nfev=100000000, xtol=eps, ftol=eps, gtol=eps, args=())["x"] - - - - q[:, :, 0] = + np.cos(p) * (f[:, :, 0]) - np.sin(p) * (f[:, :, 1]) # corrected forces - q[:, :, 1] = + np.sin(p) * (f[:, :, 0]) + np.cos(p) * (f[:, :, 1]) + p = least_squares( + fun=get_torque_angle, + x0=pstart, + method="lm", + max_nfev=100000000, + xtol=eps, + ftol=eps, + gtol=eps, + args=(), + )["x"] + + q[:, :, 0] = +np.cos(p) * (f[:, :, 0]) - np.sin(p) * ( + f[:, :, 1] + ) # corrected forces + q[:, :, 1] = +np.sin(p) * (f[:, :, 0]) + np.cos(p) * (f[:, :, 1]) return q[:, :, 0], q[:, :, 1], p # returns corrected forces and rotation angle @@ -783,8 +971,12 @@ def find_eq_position(nodes, IBC, neq): nloads = IBC.shape[0] RHSG = np.zeros((neq, 2)) - x_points = np.zeros((neq)).astype(bool) # mask showing which point has x deformation - y_points = np.zeros((neq)).astype(bool) # mask showing which point has y deformation + x_points = np.zeros((neq)).astype( + bool + ) # mask showing which point has x deformation + y_points = np.zeros((neq)).astype( + bool + ) # mask showing which point has y deformation for i in range(nloads): il = int(nodes[i, 0]) # index of the node ilx = IBC[il, 0] # indices in RHSG or fixed nodes, if -1 @@ -830,13 +1022,21 @@ def custom_solver(mat, rhs, mask_area, nodes, IBC, verbose=False): com = regionprops(mask_area.astype(int))[0].centroid # finding center of mass com = (com[1], com[0]) # as x y coordinate - c_x, c_y = np.meshgrid(range(mask_area.shape[1]), range(mask_area.shape[0])) # arrays with all x and y coordinates - r = np.zeros((mask_area.shape[0], mask_area.shape[1], 2)) # array with all positional vectors - r[:, :, 0] = c_x # Note: maybe its also enough to chose any point as reference point + c_x, c_y = np.meshgrid( + range(mask_area.shape[1]), range(mask_area.shape[0]) + ) # arrays with all x and y coordinates + r = np.zeros( + (mask_area.shape[0], mask_area.shape[1], 2) + ) # array with all positional vectors + r[:, :, 0] = ( + c_x # Note: maybe its also enough to chose any point as reference point + ) r[:, :, 1] = c_y # solidspy function that is used to construct the loads vector (rhs) nodes_xy_ordered, x_points, y_points = find_eq_position(nodes, IBC, len_disp) - r = r[nodes_xy_ordered[:, 1], nodes_xy_ordered[:, 0], :] # ordering r in the same order as rhs + r = r[ + nodes_xy_ordered[:, 1], nodes_xy_ordered[:, 0], : + ] # ordering r in the same order as rhs r = r - np.array(com) zero_disp_x[x_points] = 1 @@ -851,11 +1051,22 @@ def custom_solver(mat, rhs, mask_area, nodes, IBC, verbose=False): if type(mat) is csr_matrix: import scipy.sparse + # convert additional conditions to sparse matrix mat = scipy.sparse.vstack([mat, csr_matrix(add_matrix)], format="csr") - u_sol, error = \ - np.array(lsqr(mat, rhs, atol=10 ** -12, btol=10 ** -12, iter_lim=200000, show=verbose, conlim=10 ** 12))[ - [0, 3]] # sparse least squares solver + u_sol, error = np.array( + lsqr( + mat, + rhs, + atol=10**-12, + btol=10**-12, + iter_lim=200000, + show=verbose, + conlim=10**12, + ) + )[ + [0, 3] + ] # sparse least squares solver elif type(mat) is np.ndarray: # adding to matrix mat = np.append(mat, add_matrix, axis=0) diff --git a/saenopy/pyTFM/plotting.py b/saenopy/pyTFM/plotting.py index ce7554a..9aa6eb8 100644 --- a/saenopy/pyTFM/plotting.py +++ b/saenopy/pyTFM/plotting.py @@ -1,47 +1,44 @@ from contextlib import suppress -import copy import numpy as np import matplotlib import matplotlib.pyplot as plt -from matplotlib.ticker import MultipleLocator, ScalarFormatter from mpl_toolkits.axes_grid1.inset_locator import inset_axes -from scipy.ndimage import binary_erosion from tqdm import tqdm -def plot_continuous_boundary_stresses(plot_values, mask_boundaries=None, plot_t_vecs=False, plot_n_arrows=False, - figsize=(10, 7), - scale_ratio=0.2, border_arrow_filter=1, cbar_str="line tension in N/m", vmin=None, - vmax=None, - cbar_width="2%", cbar_height="50%", cbar_axes_fraction=0.2, - cbar_tick_label_size=20, - background_color="white", cbar_borderpad=0.1, linewidth=4, cmap="jet", - plot_cbar=True, cbar_style="clickpoints", - boundary_resolution=3, cbar_title_pad=1, outer_cb_color="grey", - outer_cb_style="-", - **kwargs): +def plot_continuous_boundary_stresses( + plot_values, + mask_boundaries=None, + plot_t_vecs=False, + plot_n_arrows=False, + figsize=(10, 7), + scale_ratio=0.2, + border_arrow_filter=1, + cbar_str="line tension in N/m", + vmin=None, + vmax=None, + cbar_width="2%", + cbar_height="50%", + cbar_axes_fraction=0.2, + # cbar_tick_label_size=20, + background_color="white", + cbar_borderpad=0.1, + linewidth=4, + cmap="jet", + plot_cbar=True, + cbar_style="clickpoints", + boundary_resolution=3, + cbar_title_pad=1, + outer_cb_color="grey", + outer_cb_style="-", +): """ plotting the line stresses (total transmitted force of cell boundaries), colored by their absolute values - as continous lines. - :param shape: - :param edge_lines: - :param lines_interpol: - :param min_v: - :param max_v: - :param mask_boundaries: - :param plot_t_vecs: - :param plot_n_arrows: - :param figsize: - :param scale_ratio: - :param arrow_filter: - :param cbar_str: - :param vmin: overwrites max_v and min_v if provided - :param vmax: overwrites max_v and min_v if provided - :return: + as continuous lines. """ - if not isinstance(plot_values[0], (list)): + if not isinstance(plot_values[0], list): plot_values = [plot_values] min_v = np.min([pv[3] for pv in plot_values]) # minimum over all objects @@ -50,22 +47,40 @@ def plot_continuous_boundary_stresses(plot_values, mask_boundaries=None, plot_t_ print("plotting cell border stresses") min_v = vmin if isinstance(vmin, (float, int)) else min_v max_v = vmax if isinstance(vmax, (float, int)) else max_v - mask_boundaries = np.zeros(shape) if not isinstance(mask_boundaries, np.ndarray) else mask_boundaries + mask_boundaries = ( + np.zeros(shape) + if not isinstance(mask_boundaries, np.ndarray) + else mask_boundaries + ) fig = plt.figure(figsize=figsize) - ax = plt.Axes(fig, [0., 0., 1., 1.]) - background_color = matplotlib.cm.get_cmap(cmap)(0) if background_color == "cmap_0" else background_color + ax = plt.Axes(fig, [0.0, 0.0, 1.0, 1.0]) + background_color = ( + matplotlib.cm.get_cmap(cmap)(0) + if background_color == "cmap_0" + else background_color + ) fig.set_facecolor(background_color) ax.set_facecolor(background_color) ax.set_axis_off() fig.add_axes(ax) + scale = 1 for shape, edge_lines, lines_interpol, *rest in plot_values: - all_t_vecs = np.vstack([subdict["t_vecs"] for subdict in lines_interpol.values()]) + all_t_vecs = np.vstack( + [subdict["t_vecs"] for subdict in lines_interpol.values()] + ) if plot_t_vecs: - scale = scale_for_quiver(all_t_vecs[:, 0], all_t_vecs[:, 1], dims=mask_boundaries.shape, - scale_ratio=scale_ratio, return_scale=True) - for line_id, interp in tqdm(lines_interpol.items(), total=len(lines_interpol.values())): + scale = scale_for_quiver( + all_t_vecs[:, 0], + all_t_vecs[:, 1], + dims=mask_boundaries.shape, + scale_ratio=scale_ratio, + return_scale=True, + ) + for line_id, interp in tqdm( + lines_interpol.items(), total=len(lines_interpol.values()) + ): p_new = interp["points_new"] x_new = p_new[:, 0] y_new = p_new[:, 1] @@ -75,14 +90,27 @@ def plot_continuous_boundary_stresses(plot_values, mask_boundaries=None, plot_t_ # plotting line segments c = matplotlib.colormaps.get_cmap(cmap)( - (t_norm - min_v) / (max_v - min_v)) # normalization and creating a color range - ## see how well that works + (t_norm - min_v) / (max_v - min_v) + ) # normalization and creating a color range + # see how well that works if line_id in edge_lines: # plot lines at the edge - plt.plot(x_new, y_new, outer_cb_style, color=outer_cb_color, linewidth=linewidth) + plt.plot( + x_new, + y_new, + outer_cb_style, + color=outer_cb_color, + linewidth=linewidth, + ) else: - for i in range(0, len(x_new) - boundary_resolution, boundary_resolution): - plt.plot([x_new[i], x_new[i + boundary_resolution]], [y_new[i], y_new[i + boundary_resolution]], - color=c[i], linewidth=linewidth) + for i in range( + 0, len(x_new) - boundary_resolution, boundary_resolution + ): + plt.plot( + [x_new[i], x_new[i + boundary_resolution]], + [y_new[i], y_new[i + boundary_resolution]], + color=c[i], + linewidth=linewidth, + ) # plotting stress vectors if plot_t_vecs: @@ -94,7 +122,13 @@ def plot_continuous_boundary_stresses(plot_values, mask_boundaries=None, plot_t_ if plot_n_arrows: for i in range(len(x_new) - 1): if i % border_arrow_filter == 0: - plt.arrow(x_new[i], y_new[i], n_vecs[i][0], n_vecs[i][1], head_width=0.5) + plt.arrow( + x_new[i], + y_new[i], + n_vecs[i][0], + n_vecs[i][1], + head_width=0.5, + ) plt.gca().invert_yaxis() # to get the typicall imshow orientation plt.xlim(0, shape[1]) @@ -102,27 +136,63 @@ def plot_continuous_boundary_stresses(plot_values, mask_boundaries=None, plot_t_ # background_color=matplotlib.cm.get_cmap(cmap)(0) if background_color=="cmap_0" else background_color # ax.set_facecolor(background_color) if plot_cbar: - add_colorbar(min_v, max_v, cmap, ax=ax, cbar_style=cbar_style, cbar_width=cbar_width, cbar_height=cbar_height, - cbar_borderpad=cbar_borderpad, v=cbar_tick_label_size, cbar_str=cbar_str, - cbar_axes_fraction=cbar_axes_fraction, cbar_title_pad=cbar_title_pad) + add_colorbar( + min_v, + max_v, + cmap, + ax=ax, + cbar_style=cbar_style, + cbar_width=cbar_width, + cbar_height=cbar_height, + cbar_borderpad=cbar_borderpad, + # v=cbar_tick_label_size, + cbar_str=cbar_str, + cbar_axes_fraction=cbar_axes_fraction, + cbar_title_pad=cbar_title_pad, + ) return fig, ax -def add_colorbar(vmin, vmax, cmap="rainbow", ax=None, cbar_style="not-clickpoints", cbar_width="2%", - cbar_height="50%", cbar_borderpad=0.1, cbar_tick_label_size=15, cbar_str="", - cbar_axes_fraction=0.2, shrink=0.8, aspect=20, cbar_title_pad=1, **kwargs): +def add_colorbar( + vmin, + vmax, + cmap="rainbow", + ax=None, + cbar_style="not-clickpoints", + cbar_width="2%", + cbar_height="50%", + cbar_borderpad=0.1, + cbar_tick_label_size=15, + cbar_str="", + cbar_axes_fraction=0.2, + shrink=0.8, + aspect=20, + cbar_title_pad=1, +): norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax) sm = plt.cm.ScalarMappable(cmap=matplotlib.colormaps.get_cmap(cmap), norm=norm) sm.set_array([]) # bug fix for lower matplotlib version if cbar_style == "clickpoints": # colorbar inside of the plot - cbaxes = inset_axes(ax, width=cbar_width, height=cbar_height, loc=5, borderpad=cbar_borderpad * 30) + cbaxes = inset_axes( + ax, + width=cbar_width, + height=cbar_height, + loc=5, + borderpad=cbar_borderpad * 30, + ) cb0 = plt.colorbar(sm, cax=cbaxes) with suppress(TypeError, AttributeError): cbaxes.set_title(cbar_str, color="white", pad=cbar_title_pad) cbaxes.tick_params(colors="white", labelsize=cbar_tick_label_size) else: # colorbar outide of the plot - cb0 = plt.colorbar(sm, aspect=aspect, shrink=shrink, fraction=cbar_axes_fraction, - pad=cbar_borderpad, ax=plt.gca()) # just exploiting the axis generation by a plt.colorbar + cb0 = plt.colorbar( + sm, + aspect=aspect, + shrink=shrink, + fraction=cbar_axes_fraction, + pad=cbar_borderpad, + ax=plt.gca(), + ) # just exploiting the axis generation by a plt.colorbar cb0.outline.set_visible(False) cb0.ax.tick_params(labelsize=cbar_tick_label_size) with suppress(TypeError, AttributeError): @@ -130,17 +200,51 @@ def add_colorbar(vmin, vmax, cmap="rainbow", ax=None, cbar_style="not-clickpoint return cb0 -def show_quiver(fx, fy, filter=[0, 1], scale_ratio=0.2, headwidth=None, headlength=None, headaxislength=None, - width=None, cmap="rainbow", - figsize=None, cbar_str="", ax=None, fig=None - , vmin=None, vmax=None, cbar_axes_fraction=0.2, cbar_tick_label_size=15 - , cbar_width="2%", cbar_height="50%", cbar_borderpad=0.1, - cbar_style="not-clickpoints", plot_style="not-clickpoints", cbar_title_pad=1, plot_cbar=True, alpha=1, - ax_origin="upper", filter_method="regular", filter_radius=5, **kwargs): +def show_quiver( + fx, + fy, + filter=None, + scale_ratio=0.2, + headwidth=None, + headlength=None, + headaxislength=None, + width=None, + cmap="rainbow", + figsize=None, + cbar_str="", + ax=None, + fig=None, + vmin=None, + vmax=None, + cbar_axes_fraction=0.2, + # cbar_tick_label_size=15, + cbar_width="2%", + cbar_height="50%", + cbar_borderpad=0.1, + cbar_style="not-clickpoints", + plot_style="not-clickpoints", + cbar_title_pad=1, + plot_cbar=True, + alpha=1, + ax_origin="upper", + filter_method="regular", + filter_radius=5, +): # list of all necessary quiver parameters - quiver_parameters = {"headwidth": headwidth, "headlength": headlength, "headaxislength": headaxislength, - "width": width, "scale_units": "xy", "angles": "xy", "scale": None} - quiver_parameters = {key: value for key, value in quiver_parameters.items() if not value is None} + if filter is None: + filter = [0, 1] + quiver_parameters = { + "headwidth": headwidth, + "headlength": headlength, + "headaxislength": headaxislength, + "width": width, + "scale_units": "xy", + "angles": "xy", + "scale": None, + } + quiver_parameters = { + key: value for key, value in quiver_parameters.items() if value is not None + } fx = fx.astype("float64") fy = fy.astype("float64") @@ -148,23 +252,43 @@ def show_quiver(fx, fy, filter=[0, 1], scale_ratio=0.2, headwidth=None, headleng if not isinstance(ax, matplotlib.axes.Axes): fig = plt.figure(figsize=figsize) ax = plt.axes() - map_values = np.sqrt(fx ** 2 + fy ** 2) + map_values = np.sqrt(fx**2 + fy**2) vmin, vmax = set_vmin_vmax(map_values, vmin, vmax) - im = plt.imshow(map_values, cmap=cmap, vmin=vmin, vmax=vmax, alpha=alpha, origin=ax_origin) # imshowing + plt.imshow( + map_values, cmap=cmap, vmin=vmin, vmax=vmax, alpha=alpha, origin=ax_origin + ) if plot_style == "clickpoints": ax.set_position([0, 0, 1, 1]) ax.set_axis_off() # plotting arrows # filtering every n-th value and every value smaller then x - fx, fy, xs, ys = filter_values(fx, fy, abs_filter=filter[0], f_dist=filter[1],filter_method=filter_method, radius=filter_radius) + fx, fy, xs, ys = filter_values( + fx, + fy, + abs_filter=filter[0], + f_dist=filter[1], + filter_method=filter_method, + radius=filter_radius, + ) if scale_ratio: # optional custom scaling with the image axis lenght fx, fy = scale_for_quiver(fx, fy, dims=dims, scale_ratio=scale_ratio) - quiver_parameters["scale"] = 1 # disabeling the auto scaling behavior of quiver + quiver_parameters["scale"] = 1 # disabeling the auto-scaling behavior of quiver plt.quiver(xs, ys, fx, fy, **quiver_parameters) # plotting the arrows if plot_cbar: - add_colorbar(vmin, vmax, cmap, ax=ax, cbar_style=cbar_style, cbar_width=cbar_width, cbar_height=cbar_height, - cbar_borderpad=cbar_borderpad, v=cbar_tick_label_size, cbar_str=cbar_str, - cbar_axes_fraction=cbar_axes_fraction, cbar_title_pad=cbar_title_pad) + add_colorbar( + vmin, + vmax, + cmap, + ax=ax, + cbar_style=cbar_style, + cbar_width=cbar_width, + cbar_height=cbar_height, + cbar_borderpad=cbar_borderpad, + # v=cbar_tick_label_size, + cbar_str=cbar_str, + cbar_axes_fraction=cbar_axes_fraction, + cbar_title_pad=cbar_title_pad, + ) return fig, ax @@ -178,28 +302,35 @@ def set_vmin_vmax(x, vmin, vmax): return vmin, vmax -def find_maxima(ar1,ar2,radius=5,shape="circle"): +def find_maxima(ar1, ar2, radius=5, shape="circle"): # generating circle - ys,xs = np.indices((radius*2 + 1,radius*2+1)) + ys, xs = np.indices((radius * 2 + 1, radius * 2 + 1)) xs = (xs - radius).astype(float) ys = (ys - radius).astype(float) - if shape=="circle": - out = np.sqrt(xs ** 2 + ys ** 2) <= radius + if shape == "circle": + out = np.sqrt(xs**2 + ys**2) <= radius xs[~out] = np.nan ys[~out] = np.nan - abs = np.sqrt(ar1**2+ar2**2) - lmax = np.unravel_index(np.nanargmax(abs),shape=abs.shape) - maxis = [lmax] + vector_abs = np.sqrt(ar1**2 + ar2**2) + lmax = np.unravel_index(np.nanargmax(vector_abs), shape=vector_abs.shape) + maxis = [lmax] while True: x_exclude = (lmax[1] + xs).flatten() y_exclude = (lmax[0] + ys).flatten() - outside_image = (x_exclude>=abs.shape[1]) | (x_exclude<0) | (y_exclude>=abs.shape[0]) | (y_exclude<0) | (np.isnan(x_exclude)) | (np.isnan(y_exclude)) + outside_image = ( + (x_exclude >= vector_abs.shape[1]) + | (x_exclude < 0) + | (y_exclude >= vector_abs.shape[0]) + | (y_exclude < 0) + | (np.isnan(x_exclude)) + | (np.isnan(y_exclude)) + ) x_exclude = x_exclude[~outside_image] y_exclude = y_exclude[~outside_image] - abs[y_exclude.astype(int),x_exclude.astype(int)] = np.nan + vector_abs[y_exclude.astype(int), x_exclude.astype(int)] = np.nan try: - lmax = np.unravel_index(np.nanargmax(abs), shape=abs.shape) + lmax = np.unravel_index(np.nanargmax(vector_abs), shape=vector_abs.shape) except ValueError: break maxis.append(lmax) @@ -212,20 +343,14 @@ def find_maxima(ar1,ar2,radius=5,shape="circle"): def filter_values(ar1, ar2, abs_filter=0, f_dist=3, filter_method="regular", radius=5): """ function to filter out values from an array for better display - :param ar1: - :param ar2: - :param ar: - :param f_dist: distance betweeen filtered values - :return: """ - if filter_method == "regular": pixx = np.arange(np.shape(ar1)[0]) pixy = np.arange(np.shape(ar1)[1]) xv, yv = np.meshgrid(pixy, pixx) - def_abs = np.sqrt((ar1 ** 2 + ar2 ** 2)) + def_abs = np.sqrt((ar1**2 + ar2**2)) select_x = ((xv - 1) % f_dist) == 0 select_y = ((yv - 1) % f_dist) == 0 select_size = def_abs > abs_filter @@ -234,19 +359,21 @@ def filter_values(ar1, ar2, abs_filter=0, f_dist=3, filter_method="regular", rad s2 = ar2[select] x_ind = xv[select] y_ind = yv[select] - if filter_method == "local_maxima": - y_ind,x_ind = find_maxima(ar1, ar2, radius=radius,shape="circle") + elif filter_method == "local_maxima": + y_ind, x_ind = find_maxima(ar1, ar2, radius=radius, shape="circle") s1 = ar1[y_ind, x_ind] s2 = ar2[y_ind, x_ind] - if filter_method == "local_maxima_square": - y_ind,x_ind = find_maxima(ar1, ar2, radius=radius,shape="square") + elif filter_method == "local_maxima_square": + y_ind, x_ind = find_maxima(ar1, ar2, radius=radius, shape="square") s1 = ar1[y_ind, x_ind] s2 = ar2[y_ind, x_ind] + else: + raise ValueError("Filter method unknown", filter_method) return s1, s2, x_ind, y_ind def scale_for_quiver(ar1, ar2, dims, scale_ratio=0.2, return_scale=False): - scale = scale_ratio * np.max(dims) / np.nanmax(np.sqrt((ar1) ** 2 + (ar2) ** 2)) + scale = scale_ratio * np.max(dims) / np.nanmax(np.sqrt(ar1 ** 2 + ar2 ** 2)) if return_scale: return scale return ar1 * scale, ar2 * scale diff --git a/saenopy/pyTFM/stress_functions.py b/saenopy/pyTFM/stress_functions.py index 42f93f8..9a30fcc 100644 --- a/saenopy/pyTFM/stress_functions.py +++ b/saenopy/pyTFM/stress_functions.py @@ -7,15 +7,14 @@ def exclude_by_key(d, ex_list): - ex_dict = {key: values for key, values in d.items() if - key not in ex_list} + ex_dict = {key: values for key, values in d.items() if key not in ex_list} return ex_dict def normal_vectors_from_splines(u, tck): """ computes the normal vector from a spline representation of a curve at the points defined by parameters u. - The normal vector is the derivative, with [dx/du,dy/du]-->[dy/du,-dx/du]). + The normal vector is the derivative, with [dx/du,dy/du]-->[dy/du,-dx/du]. (Derivative is parallel vector to spline curve) :param u: parameters defining on which points to evaluate the splines :param tck: spline @@ -23,15 +22,20 @@ def normal_vectors_from_splines(u, tck): """ n_vectors = np.array( - splev(u, tck, der=1)).T # evaluation of the spline derivative, only at the points inx and y + splev(u, tck, der=1) + ).T # evaluation of the spline derivative, only at the points inx and y # (returns result as x,y value array - n_vectors[:, [1, 0]] = n_vectors[:, [0, 1]] # normal vectors by switching x andy and tacking negative x values + n_vectors[:, [1, 0]] = n_vectors[ + :, [0, 1] + ] # normal vectors by switching x andy and tacking negative x values n_vectors[:, 1] *= -1 n_vectors = n_vectors / np.linalg.norm(n_vectors, axis=1)[:, None] # normalizing return n_vectors -def lineTension(lines_splines, line_lengths, stress_tensor, pixel_length, interpol_factor=1): +def line_tension( + lines_splines, line_lengths, stress_tensor, pixel_length, interpol_factor=1 +): """ function to perform interpolation on lines, to get new x,y coordinates, calculate the normal vectors on these points and calculate the stress vector and the norm of the stress vectors, across the lines at the interpolated points. @@ -46,13 +50,19 @@ def lineTension(lines_splines, line_lengths, stress_tensor, pixel_length, interp norm of the stress vectors from the interpolation min_v,max_v: maximal and minimal values of the stress_vector norm. This is used to get a uniform color bar when plotting later. + + Parameters + ---------- + line_lengths """ # interpolating the stress vector: - pixx = np.linspace(0, 1, stress_tensor.shape[1]) # coordinate space on which to perform the interpolation + pixx = np.linspace( + 0, 1, stress_tensor.shape[1] + ) # coordinate space on which to perform the interpolation pixy = np.linspace(0, 1, stress_tensor.shape[0]) - # using 2 dimensional interpolation on each component + # using 2-dimensional interpolation on each component # this returns a continuous function f(x,y)=z, that can be evaluated at any point x,y sig_xx_inter = RectBivariateSpline(pixx, pixy, stress_tensor[:, :, 0, 0].T) sig_yx_inter = RectBivariateSpline(pixx, pixy, stress_tensor[:, :, 1, 0].T) @@ -73,33 +83,42 @@ def lineTension(lines_splines, line_lengths, stress_tensor, pixel_length, interp x_new, y_new = splev(u_new, tck, der=0) # interpolating new x,y points p_new = np.vstack([x_new, y_new]).T # convenient form for new points - n_vecs = normal_vectors_from_splines(u_new, - tck) # new normal_vectors, using the derivative at interpolation points + n_vecs = normal_vectors_from_splines( + u_new, tck + ) # new normal_vectors, using the derivative at interpolation points # traction vectors using interpolation of the stress tensor - t_vecs, t_norm, t_vecs_n, t_vecs_shear = stress_vector_from_tensor_interpolation(p_new, n_vecs, sig_xx_inter, - sig_yx_inter, - sig_yy_inter, inter_ranges) + t_vecs, t_norm, t_vecs_n, t_vecs_shear = ( + stress_vector_from_tensor_interpolation( + p_new, n_vecs, sig_xx_inter, sig_yx_inter, sig_yy_inter, inter_ranges + ) + ) # conversion to N/m (height not included) - t_vecs, t_norm, t_vecs_n, t_vecs_shear = [x / (pixel_length * 10 ** -6) for x in - [t_vecs, t_norm, t_vecs_n, t_vecs_shear]] + t_vecs, t_norm, t_vecs_n, t_vecs_shear = [ + x / (pixel_length * 10**-6) + for x in [t_vecs, t_norm, t_vecs_n, t_vecs_shear] + ] # updating minimal value to find global minimum eventually min_v = np.min([min_v, np.nanmin(t_norm)]) max_v = np.max([max_v, np.nanmax(t_norm)]) - lines_interpol[i] = {"points_new": p_new, # array of new points from interpolation - "t_vecs": t_vecs, # stress vectors at the new points - "t_norm": t_norm, # norm of stress vectors at the new points - "n_vecs": n_vecs, # normal vectors at the new points - "t_normal": t_vecs_n, # normal component of stress vectors at the new points - "t_shear": t_vecs_shear} # shear component of stress vectors at the new points + lines_interpol[i] = { + "points_new": p_new, # array of new points from interpolation + "t_vecs": t_vecs, # stress vectors at the new points + "t_norm": t_norm, # norm of stress vectors at the new points + "n_vecs": n_vecs, # normal vectors at the new points + "t_normal": t_vecs_n, # normal component of stress vectors at the new points + "t_shear": t_vecs_shear, + } # shear component of stress vectors at the new points return lines_interpol, min_v, max_v -def stress_vector_from_tensor_interpolation(ps, n_vecs, sig_xx_inter, sig_yx_inter, sig_yy_inter, inter_ranges): +def stress_vector_from_tensor_interpolation( + ps, n_vecs, sig_xx_inter, sig_yx_inter, sig_yy_inter, inter_ranges +): """ calculates the stress vector with t=sigma*n ; sigma: stress tensor, n: normal vector to the cut over which the stress tensor is calculated. This function uses 2-interpolation functions for the individual stress components, @@ -117,32 +136,50 @@ def stress_vector_from_tensor_interpolation(ps, n_vecs, sig_xx_inter, sig_yx_int t_norm: norm of these vectors """ - u = copy.deepcopy(ps) # points represented on the (0,1) interval used for the interpolation function. + u = copy.deepcopy( + ps + ) # points represented on the (0,1) interval used for the interpolation function. u[:, 0] = u[:, 0] / inter_ranges[0] u[:, 1] = u[:, 1] / inter_ranges[1] t_vecs = [] for i in range(len(u)): # stress vectors according to cauchy theorem - t_vec = [sig_xx_inter(u[i][0], u[i][1]) * n_vecs[i][0] + sig_yx_inter(u[i][0], u[i][1]) * n_vecs[i][1], - sig_yx_inter(u[i][0], u[i][1]) * n_vecs[i][0] + sig_yy_inter(u[i][0], u[i][1]) * n_vecs[i][1]] + t_vec = [ + sig_xx_inter(u[i][0], u[i][1]) * n_vecs[i][0] + + sig_yx_inter(u[i][0], u[i][1]) * n_vecs[i][1], + sig_yx_inter(u[i][0], u[i][1]) * n_vecs[i][0] + + sig_yy_inter(u[i][0], u[i][1]) * n_vecs[i][1], + ] t_vecs.append(t_vec) t_vecs = np.array(t_vecs).squeeze() # removing unused array dimension t_norm = np.linalg.norm(t_vecs, axis=1) # calculating norm of the stress vector - t_vec_n = np.abs(np.sum(t_vecs * n_vecs, axis=1)) # length of the normal component of the line tension - t_vec_shear = np.sqrt(t_norm ** 2 - t_vec_n ** 2) # shear component of the line tension + t_vec_n = np.abs( + np.sum(t_vecs * n_vecs, axis=1) + ) # length of the normal component of the line tension + t_vec_shear = np.sqrt(t_norm**2 - t_vec_n**2) # shear component of the line tension return t_vecs, t_norm, t_vec_n, t_vec_shear def calculate_stress_tensor(s_nodes, nodes, dims=None): if not dims: - stress_tensor = np.zeros((int(np.sqrt(len(s_nodes))), int(np.sqrt(len(s_nodes))), 2, 2)) + stress_tensor = np.zeros( + (int(np.sqrt(len(s_nodes))), int(np.sqrt(len(s_nodes))), 2, 2) + ) else: stress_tensor = np.zeros((dims[0], dims[1], 2, 2)) - stress_tensor[nodes[:, 2].astype(int), nodes[:, 1].astype(int), 0, 0] = s_nodes[:, 0] # sigma_x - stress_tensor[nodes[:, 2].astype(int), nodes[:, 1].astype(int), 1, 1] = s_nodes[:, 1] # sigma_y - stress_tensor[nodes[:, 2].astype(int), nodes[:, 1].astype(int), 1, 0] = s_nodes[:, 2] # sigma_yx - stress_tensor[nodes[:, 2].astype(int), nodes[:, 1].astype(int), 0, 1] = s_nodes[:, 2] # sigma_xy + stress_tensor[nodes[:, 2].astype(int), nodes[:, 1].astype(int), 0, 0] = s_nodes[ + :, 0 + ] # sigma_x + stress_tensor[nodes[:, 2].astype(int), nodes[:, 1].astype(int), 1, 1] = s_nodes[ + :, 1 + ] # sigma_y + stress_tensor[nodes[:, 2].astype(int), nodes[:, 1].astype(int), 1, 0] = s_nodes[ + :, 2 + ] # sigma_yx + stress_tensor[nodes[:, 2].astype(int), nodes[:, 1].astype(int), 0, 1] = s_nodes[ + :, 2 + ] # sigma_xy return stress_tensor @@ -152,7 +189,8 @@ def coefficient_of_variation(mask, x, border_pad=0): if border_pad > 0: mask_cp = binary_erosion(mask_cp, iterations=border_pad) return np.nanstd(x[mask_cp]) / np.abs( - np.nanmean(x[mask_cp])) # absolute value of mean is an alternative definition + np.nanmean(x[mask_cp]) + ) # absolute value of mean is an alternative definition def all_stress_measures(st, px_size=1): @@ -161,11 +199,11 @@ def all_stress_measures(st, px_size=1): tau_xy = st[:, :, 0, 1] # shear stress # principal (normal stresses) - sigma_max = (sig_x + sig_y) / 2 + np.sqrt(((sig_x - sig_y) / 2) ** 2 + tau_xy ** 2) - sigma_min = (sig_x + sig_y) / 2 - np.sqrt(((sig_x - sig_y) / 2) ** 2 + tau_xy ** 2) + sigma_max = (sig_x + sig_y) / 2 + np.sqrt(((sig_x - sig_y) / 2) ** 2 + tau_xy**2) + sigma_min = (sig_x + sig_y) / 2 - np.sqrt(((sig_x - sig_y) / 2) ** 2 + tau_xy**2) sigma_max_abs = np.maximum(np.abs(sigma_max), np.abs(sigma_min)) # maximum shear stress - tau_max = np.sqrt(((sig_x - sig_y) / 2) ** 2 + tau_xy ** 2) + tau_max = np.sqrt(((sig_x - sig_y) / 2) ** 2 + tau_xy**2) # angle of maximal principal stress phi_n = np.arctan(2 * tau_xy / (sig_x - sig_y)) / 2 # angel of maximal shear stress @@ -173,11 +211,20 @@ def all_stress_measures(st, px_size=1): # side note: (phi_n-phi_shear) = pi/4 should always hold sigma_mean = (sigma_max + sigma_min) / 2 # mean normal stress - return sigma_max / px_size, sigma_min / px_size, sigma_max_abs / px_size, tau_max / px_size, \ - phi_n, phi_shear, sigma_mean / px_size + return ( + sigma_max / px_size, + sigma_min / px_size, + sigma_max_abs / px_size, + tau_max / px_size, + phi_n, + phi_shear, + sigma_mean / px_size, + ) -def reorder_vectors_inward(borders, lines_interpol, cell_id, line_ids, plot_n_vecs=False, mask_boundaries=None): +def reorder_vectors_inward( + borders, lines_interpol, cell_id, line_ids, plot_n_vecs=False, mask_boundaries=None +): """ reorientation of normal and traction force vectors, so that the normal vectors of one cell all point inwards. :param borders: @@ -189,21 +236,30 @@ def reorder_vectors_inward(borders, lines_interpol, cell_id, line_ids, plot_n_ve :return: """ cell_area = borders.cells_area[cell_id] - n_vectors = {line_id: lines_interpol[line_id]["n_vecs"] for line_id in - line_ids} # extracting relevant normal vectors - t_vectors = {line_id: lines_interpol[line_id]["t_vecs"] for line_id in - line_ids} # extracting relevant normal vectors - points_dict = {line_id: lines_interpol[line_id]["points_new"] for line_id in - line_ids} # extracting relevant normal vectors + n_vectors = { + line_id: lines_interpol[line_id]["n_vecs"] for line_id in line_ids + } # extracting relevant normal vectors + t_vectors = { + line_id: lines_interpol[line_id]["t_vecs"] for line_id in line_ids + } # extracting relevant normal vectors + points_dict = { + line_id: lines_interpol[line_id]["points_new"] for line_id in line_ids + } # extracting relevant normal vectors for l_id in line_ids: - nps1 = np.round(points_dict[l_id] + n_vectors[l_id]).astype(int) # predict points original orientation - nps2 = np.round(points_dict[l_id] - n_vectors[l_id]).astype(int) # reversed orientation + nps1 = np.round(points_dict[l_id] + n_vectors[l_id]).astype( + int + ) # predict points original orientation + nps2 = np.round(points_dict[l_id] - n_vectors[l_id]).astype( + int + ) # reversed orientation s1 = np.sum(cell_area[nps1[:, 1], nps1[:, 0]]) s2 = np.sum(cell_area[nps2[:, 1], nps2[:, 0]]) change_orientation = s1 < s2 - n_vectors[l_id] *= (-(change_orientation * 2 - 1)) # changes orientation inwards or keep it - t_vectors[l_id] *= (-(change_orientation * 2 - 1)) # change t vector accordingly + n_vectors[l_id] *= -( + change_orientation * 2 - 1 + ) # changes orientation inwards or keep it + t_vectors[l_id] *= -(change_orientation * 2 - 1) # change t vector accordingly # confirmation of correct vector orientation if plot_n_vecs: @@ -222,12 +278,17 @@ def reorder_vectors_inward(borders, lines_interpol, cell_id, line_ids, plot_n_ve def normal_and_shear(t_vecs, n_vecs): tn = np.sum(t_vecs * n_vecs, axis=1) # dot product of t vec and n vec # remaining component by triangle equation (Pythagorean theorem) - ts = np.sqrt(np.sum(t_vecs * t_vecs, axis=1) - tn ** 2) + ts = np.sqrt(np.sum(t_vecs * t_vecs, axis=1) - tn**2) return tn, ts -def mean_stress_vector_norm(lines_interpolation, borders, exclude_colony_edge=True, norm_level="points", - vtype="t_norm"): +def mean_stress_vector_norm( + lines_interpolation, + borders, + exclude_colony_edge=True, + norm_level="points", + vtype="t_norm", +): """ average norm of the stress vector.First some vectors are averaged. Averaging can be eft out (level="points"), for each line(level="lines") or for each cell (level="cells") @@ -246,15 +307,27 @@ def mean_stress_vector_norm(lines_interpolation, borders, exclude_colony_edge=Tr all_values = [] if norm_level == "cells": # only makes sense if normal vectors are orientated all in the same manner - for cell_id, line_ids in borders.cells_lines.items(): # only uses lines associated to cells + for ( + cell_id, + line_ids, + ) in borders.cells_lines.items(): # only uses lines associated to cells # orientation correction: - n_vectors, t_vectors = reorder_vectors_inward(borders, lines_interpolation, cell_id, line_ids, - plot_n_vecs=False) + n_vectors, t_vectors = reorder_vectors_inward( + borders, lines_interpolation, cell_id, line_ids, plot_n_vecs=False + ) # optional exclude borders at the colony edge - exclude_cond = borders.edge_lines if exclude_colony_edge else list(lines_interpolation.keys()) + exclude_cond = ( + borders.edge_lines + if exclude_colony_edge + else list(lines_interpolation.keys()) + ) # excluding some lines (or none) and joining to single array - t_vectors = np.vstack(list(exclude_by_key(t_vectors, exclude_cond).values())) - n_vectors = np.vstack(list(exclude_by_key(n_vectors, exclude_cond).values())) + t_vectors = np.vstack( + list(exclude_by_key(t_vectors, exclude_cond).values()) + ) + n_vectors = np.vstack( + list(exclude_by_key(n_vectors, exclude_cond).values()) + ) # calculating normal and shear components of traction vector after reorientation tns, tss = normal_and_shear(t_vectors, n_vectors) single_cell_force = None @@ -273,9 +346,16 @@ def mean_stress_vector_norm(lines_interpolation, borders, exclude_colony_edge=Tr if exclude_colony_edge: lines_interpolation = exclude_by_key(lines_interpolation, borders.edge_lines) if norm_level == "lines": # mean over a line - all_values = np.vstack([np.mean(sub_dict[vtype], axis=0) for sub_dict in lines_interpolation.values()]) + all_values = np.vstack( + [ + np.mean(sub_dict[vtype], axis=0) + for sub_dict in lines_interpolation.values() + ] + ) if norm_level == "points": # each point individually - all_values = np.concatenate([sub_dict[vtype] for sub_dict in lines_interpolation.values()]) + all_values = np.concatenate( + [sub_dict[vtype] for sub_dict in lines_interpolation.values()] + ) # returning the norm of the mean t_vector all_values = np.linalg.norm(all_values, axis=1) if vtype == "t_vecs" else all_values diff --git a/saenopy/pyTFM/suppress_warnings.py b/saenopy/pyTFM/suppress_warnings.py index 20d4825..0610654 100644 --- a/saenopy/pyTFM/suppress_warnings.py +++ b/saenopy/pyTFM/suppress_warnings.py @@ -1,13 +1,12 @@ import warnings -class suppress_warnings(): +class suppress_warnings: def __init__(self, warning_type): self.warning_type = warning_type def __enter__(self): warnings.filterwarnings("ignore", category=self.warning_type) - def __exit__(self, type, value, traceback): + def __exit__(self, warn_type, value, traceback): warnings.filterwarnings("default", category=self.warning_type) - diff --git a/saenopy/pyTFM/utilities_TFM.py b/saenopy/pyTFM/utilities_TFM.py index 68d8041..a02b9b7 100755 --- a/saenopy/pyTFM/utilities_TFM.py +++ b/saenopy/pyTFM/utilities_TFM.py @@ -1,17 +1,28 @@ import numpy as np +import copy -def join_dictionary(d1, d2, update_keys=False): - if update_keys: +def update_keys(d1, d2): + # recounts keys of d2 by starting with last key of d1. keys must all be integers + d3 = copy.deepcopy(d2) + max_key = np.max(list(d1.keys())) + for key, value in d2.items(): + max_key += 1 + d3[max_key] = value + return d3 + + +def join_dictionary(d1, d2, do_update_keys=False): + if do_update_keys: d3 = update_keys(d1, d2) else: d3 = d2 return {**d1, **d3} - # note:z = {**x, **y} and "update" are nices tricks here + # note:z = {**x, **y} and "update" are nice tricks here def make_random_discrete_color_range(size): colors = [] for i in range(size): - colors.append('#%06X' % np.random.randint(0, 0xFFFFFF)) + colors.append("#%06X" % np.random.randint(0, 0xFFFFFF)) return colors