diff --git a/Data/empty.txt b/Data/empty.txt new file mode 100644 index 0000000..89e5656 --- /dev/null +++ b/Data/empty.txt @@ -0,0 +1 @@ +emptyfile diff --git a/FiberFit_toshare_3_4_2021.ipynb b/FiberFit_toshare_3_4_2021.ipynb new file mode 100644 index 0000000..fbdb82e --- /dev/null +++ b/FiberFit_toshare_3_4_2021.ipynb @@ -0,0 +1,398 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Fiber Analysis; Modified from FiberFit Software\n", + "## software is from [this paper](https://link.springer.com/article/10.1007/s10237-016-0776-3) and Github with original program can be viewed [here.](https://github.com/NTMatBoiseState/FiberFit/tree/master/src/fiberfit_model)\n", + "#### this software below is same as above, but it allows for bulk proccessing and provides export of raw data\n", + "\n", + "### Steps: First, Upload your Images for analysis." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/srv/conda/envs/notebook/lib/python3.7/site-packages/IPython/html.py:14: ShimWarning: The `IPython.html` package has been deprecated since IPython 4.0. You should import from `notebook` instead. `IPython.html.widgets` has moved to `ipywidgets`.\n", + " \"`IPython.html.widgets` has moved to `ipywidgets`.\", ShimWarning)\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "448a008d479d43e7acb8b816338c4bfe", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "FileUpload(value={}, description='Upload', multiple=True, style=ButtonStyle(button_color='lightblue'))" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "from widgetsmp import s,run,run2,save,clear\n", + "from src.foruploading5 import ff,upload,prepimages\n", + "import ipywidgets as widgets\n", + "from IPython.display import display, Markdown, clear_output\n", + "from ipywidgets import interact, interact_manual\n", + "import os\n", + "import glob\n", + "%matplotlib inline\n", + "## Remove blank file from output folder\n", + "files = glob.glob('Data/*',recursive=True)\n", + "for f in files:\n", + " try:\n", + " os.remove(f)\n", + " except OSError as e:\n", + " print(\"Error: %s : %s\" % (f, e.strerror))\n", + "s = widgets.FileUpload(accept='', multiple=True)\n", + "upload(s)\n", + "s" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Next, Run Code. Input maximum and minimum fiber thickness (unit - pixels) for fibers. Suggested values are shown in boxes, but still need to be typed in the box. Output will be in \"output\" folder on home page. CSV output provides alignment data in order of images uploaded." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "22a99f3dcad74e48b9f794dca567e39a", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "VBox(children=(Text(value='', description='Max fiber Thickness (pixels)', layout=Layout(width='300px'), placeh…" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "17d380b8fb844c878356cacc729c1a94", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Output()" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "run(_)### Next, Run Code. Input maximum and minimum fiber thickness for fibers. Suggested values are shown in boxes, but still need to be typed in the box. Output will be in \"output\" folder on home page. Can Download images and .csv from there (working on including .zip file download)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Next, Create zip file to download." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "output/data_export.csv\r\n", + "output/FFT Norm Test Image_90_0.7_90.03_0.7_0.5729723.png\r\n", + "output/GroupedNorm Test Image_90_0.7_90.03_0.7_0.5729723.png\r\n", + "output/Normdist Norm Test Image_90_0.7_90.03_0.7_0.5729723.png\r\n", + "output/polarpositiveNorm Test Image_90_0.7_90.03_0.7_0.5729723.png\r\n" + ] + } + ], + "source": [ + "!tar chvfz notebook.tar.gz output/*;" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Finally, download output file by clicking on download link. Files include images for analysis and bulk .csv file. Delete files in Data and output file by clicking on clear data button prior to running a new set of files." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "0326fc2d463f4487ae43dc6b83407cb0", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "HTML(value='\\n \\n \\n …" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "f0cacc9ed7d44e49a62ead2ffea602bd", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "VBox(children=(Button(description='Clear Data', layout=Layout(width='300px'), style=ButtonStyle()),))" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "5ab496d2336e4ca0b3de20750368f736", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Output()" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "save(_)\n", + "clear(_)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Next, if desired, run code below for color maps of fiber orientation and diameter thickness (code modified from [this paper](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0215137). Again, output will be in \"output\" folder. \n", + "###### Note: These color maps are from a program called Quanfima. This is a method that utilizes alternative calculation method than above. It takes longer to calculate, but can be used qualitatively against csv + image outputs in above program to confirm if similar results. PNG or JPEG files are suggested for speedier calculations." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Steps: First, Upload a single image for analysis" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "02033485ccb64585b7172f35bed4c8b0", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "FileUpload(value={}, description='Upload', multiple=True, style=ButtonStyle(button_color='lightblue'))" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "s = widgets.FileUpload(accept='', multiple=True)\n", + "upload(s)\n", + "s" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Steps: Next, run analysis" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "90747e3bec0f4fd98a67ebd80f4be938", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "VBox(children=(Button(description='Run Fiber Orientation Visualization', layout=Layout(width='300px'), style=B…" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "eb5c6f7237d1431a84ebbf21a8d8c89c", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Output()" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "run2(_)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Create output file to download" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "output/DIAMETER_estimate_diameter_map.png\r\n", + "output/orientation_map.png\r\n" + ] + } + ], + "source": [ + "!tar chvfz notebook.tar.gz output/*;" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Download output file by clicking on download link. Files include images shown above. Click clear data to clear input and output data to run again." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "15b97f9033a74a96acdca2e904f7c1a6", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "HTML(value='\\n \\n \\n …" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "ff16db18b43242a591a73b8c9611b86b", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "VBox(children=(Button(description='Clear Data', layout=Layout(width='300px'), style=ButtonStyle()),))" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "18a3abad11f4459999349b409fbaf7e7", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Output()" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "save(_)\n", + "clear(_)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.9" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/apt.txt b/apt.txt new file mode 100644 index 0000000..f359f07 --- /dev/null +++ b/apt.txt @@ -0,0 +1 @@ +libgl1-mesa-glx diff --git a/environment.yml b/environment.yml new file mode 100644 index 0000000..e352463 --- /dev/null +++ b/environment.yml @@ -0,0 +1,44 @@ +name: test-environment +channels: + - conda-forge + - defaults +dependencies: + - Cycler + - decorator + - greenlet + - GitPython + - statsmodels + - imageio + - ipykernel + - ipython_genutils + - ipython + - ipywidgets + - jedi + - matplotlib==3.1.3 + - more-itertools + - numpy + - pandas + - pathlib + - pip + - python-dateutil>=2.7 + - scikit-image + - scipy + - seaborn + - vigra=1.11.1 + - pip: + - ffmpy + - isapi_wsgi + - jupyterlab + - opencv-python + - mpl-utils + - mpl4qt + - mpl-tools + - PyQt5 + - py-dateutil + - runtime + - cv_utils + - cv2studio + - smop + - syspath + + diff --git a/output/blank.txt b/output/blank.txt new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/output/blank.txt @@ -0,0 +1 @@ + diff --git a/quanfima/.gitignore b/quanfima/.gitignore new file mode 100644 index 0000000..0d327a6 --- /dev/null +++ b/quanfima/.gitignore @@ -0,0 +1,9 @@ +# Compiled python modules. +*.pyc + +# Setuptools distribution folder. +/dist/ + +# Python egg metadata, regenerated from source files by setuptools. +/*.egg-info +/*.egg diff --git a/quanfima/.travis.yml b/quanfima/.travis.yml new file mode 100644 index 0000000..3a3d687 --- /dev/null +++ b/quanfima/.travis.yml @@ -0,0 +1,11 @@ +language: python + +python: + - "2.7" + +install: + - pip install -r requirements.txt + - pip install . + +script: + - nosetests diff --git a/quanfima/LICENSE b/quanfima/LICENSE new file mode 100644 index 0000000..9825b19 --- /dev/null +++ b/quanfima/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2018 Roman Shkarin, Andrei Shkarin + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/quanfima/MANIFEST.in b/quanfima/MANIFEST.in new file mode 100644 index 0000000..9561fb1 --- /dev/null +++ b/quanfima/MANIFEST.in @@ -0,0 +1 @@ +include README.rst diff --git a/quanfima/README.rst b/quanfima/README.rst new file mode 100644 index 0000000..083e361 --- /dev/null +++ b/quanfima/README.rst @@ -0,0 +1,95 @@ +.. image:: docs/source/_static/logo.png + :align: left + +----------- + +.. image:: https://travis-ci.org/rshkarin/quanfima.svg?branch=master + :target: https://travis-ci.org/rshkarin/quanfima + +.. image:: https://readthedocs.org/projects/quanfima/badge/?version=latest + :target: http://quanfima.readthedocs.io/en/latest/?badge=latest + :alt: Documentation Status + +.. image:: https://zenodo.org/badge/127795855.svg + :target: https://zenodo.org/badge/latestdoi/127795855 + + + + +*Quanfima* (**qu**\ antitative **an**\ alysis of **fi**\ brous **ma**\ terials) +is a collection of useful functions for morphological analysis and visualization +of 2D/3D data from various areas of material science. The aim is to simplify +the analysis process by providing functionality for frequently required tasks +in the same place. + +More examples of usage you can find in the documentation. + +- Analysis of fibrous structures by tensor-based method in 2D / 3D datasets. +- Estimation of structure diameters in 2D / 3D by a ray-casting method. +- Counting of particles in 2D / 3D datasets and providing a detailed report in + pandas.DataFrame format. +- Calculation of porosity measure for each material in 2D / 3D datasets. +- Visualization in 2D / 3D using matplotlib, visvis packages. + +Installation +------------ + +The easiest way to install the latest version is by using pip:: + + $ pip install quanfima + +You may also use Git to clone the repository and install it manually:: + + $ git clone https://github.com/rshkarin/quanfima.git + $ cd quanfima + $ python setup.py install + +Usage +----- +Open a grayscale image, perform segmentation, estimate porosity, analyze fiber +orientation and diameters, and plot the results. + +.. code-block:: python + + import numpy as np + from skimage import io, filters + from quanfima import morphology as mrph + from quanfima import visualization as vis + from quanfima import utils + + img = io.imread('../data/polymer_slice.tif') + + th_val = filters.threshold_otsu(img) + img_seg = (img > th_val).astype(np.uint8) + + # estimate porosity + pr = mrph.calc_porosity(img_seg) + for k,v in pr.items(): + print 'Porosity ({}): {}'.format(k, v) + + # prepare data and analyze fibers + data, skeleton, skeleton_thick = utils.prepare_data(img_seg) + cskel, fskel, omap, dmap, ovals, dvals = \ + mrph.estimate_fiber_properties(data, skeleton) + + # plot results + vis.plot_orientation_map(omap, fskel, min_label=u'0°', max_label=u'180°', + figsize=(10,10), + name='2d_polymer', + output_dir='/path/to/output/dir') + vis.plot_diameter_map(dmap, cskel, figsize=(10,10), cmap='gist_rainbow', + name='2d_polymer', + output_dir='/path/to/output/dir') + +.. code-block:: python + + >> Porosity (Material 1): 0.845488888889 + +.. image:: docs/source/_static/2d_polymer_data.png + :align: center + +.. image:: docs/source/_static/2d_polymer_orientation_map_600px.png + :align: center + +.. image:: docs/source/_static/2d_polymer_diameter_map_600px.png + :align: center diff --git a/quanfima/data/polymer3d_8bit_128x128x128.raw b/quanfima/data/polymer3d_8bit_128x128x128.raw new file mode 100644 index 0000000..05560e6 Binary files /dev/null and b/quanfima/data/polymer3d_8bit_128x128x128.raw differ diff --git a/quanfima/data/polymer_slice.tif b/quanfima/data/polymer_slice.tif new file mode 100644 index 0000000..f4dfc08 Binary files /dev/null and b/quanfima/data/polymer_slice.tif differ diff --git a/quanfima/data/test.txt b/quanfima/data/test.txt new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/quanfima/data/test.txt @@ -0,0 +1 @@ + diff --git a/quanfima/docs/Makefile b/quanfima/docs/Makefile new file mode 100644 index 0000000..2c30804 --- /dev/null +++ b/quanfima/docs/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line. +SPHINXOPTS = +SPHINXBUILD = sphinx-build +SPHINXPROJ = quanfima +SOURCEDIR = source +BUILDDIR = build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/quanfima/docs/make.bat b/quanfima/docs/make.bat new file mode 100644 index 0000000..41fc86d --- /dev/null +++ b/quanfima/docs/make.bat @@ -0,0 +1,36 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=source +set BUILDDIR=build +set SPHINXPROJ=quanfima + +if "%1" == "" goto help + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.http://sphinx-doc.org/ + exit /b 1 +) + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% + +:end +popd diff --git a/quanfima/docs/source/change_log.rst b/quanfima/docs/source/change_log.rst new file mode 100644 index 0000000..27a61f8 --- /dev/null +++ b/quanfima/docs/source/change_log.rst @@ -0,0 +1,32 @@ +========= +ChangeLog +========= + +Version 0.1a1 +============== + +Features: + +- Orientation and diameter estimation for 2D / 3D data +- Counting of particles with subsequent estimation of properties +- Porosity calculation for ND data +- Visualization of 2D orientation and diameter with MPL package +- Visualization of 3D orientation and diameter with VisVis package +- Simulation of fibers and particles with various scenarios + + +Version 0.1a2 +============== + +Fix: + +- Incorrect display of measure units at a legend of a pie chart + + +Version 0.1a3 +============== + +New features: + +- Block-wise orientation estimation with 2D FFT and PCA-approach [estimate_fourier_orientation] +- Visualization of block-wise orientation with MPL package [plot_fourier_orientation] diff --git a/quanfima/docs/source/conf.py b/quanfima/docs/source/conf.py new file mode 100644 index 0000000..3483c70 --- /dev/null +++ b/quanfima/docs/source/conf.py @@ -0,0 +1,170 @@ +# -*- coding: utf-8 -*- +# +# Configuration file for the Sphinx documentation builder. +# +# This file does only contain a selection of the most common options. For a +# full list see the documentation: +# http://www.sphinx-doc.org/en/stable/config + +# -- Path setup -------------------------------------------------------------- + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +# +import os +import sys + +sys.path.insert(0, os.path.abspath('../..')) +print os.path.abspath('../..') + + +# -- Project information ----------------------------------------------------- + +project = u'quanfima' +copyright = u'2018, Roman Shkarin, Andrei Shkarin' +author = u'Roman Shkarin, Andrei Shkarin' + +# The short X.Y version +version = u'' +# The full version, including alpha/beta/rc tags +release = u'0.1a1' + + +# -- General configuration --------------------------------------------------- + +# If your documentation needs a minimal Sphinx version, state it here. +# +# needs_sphinx = '1.0' + +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom +# ones. +extensions = [ + 'sphinx.ext.autodoc', + 'sphinx.ext.doctest', + 'sphinx.ext.coverage', + 'sphinx.ext.mathjax', + 'sphinx.ext.viewcode', + 'sphinx.ext.githubpages', + # 'autoapi.extension' +] + +# Add any paths that contain templates here, relative to this directory. +templates_path = ['_templates'] + +# The suffix(es) of source filenames. +# You can specify multiple suffix as a list of string: +# +# source_suffix = ['.rst', '.md'] +source_suffix = '.rst' + +# The master toctree document. +master_doc = 'index' + +# The language for content autogenerated by Sphinx. Refer to documentation +# for a list of supported languages. +# +# This is also used if you do content translation via gettext catalogs. +# Usually you set "language" from the command line for these cases. +language = None + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +# This pattern also affects html_static_path and html_extra_path . +exclude_patterns = [] + +# The name of the Pygments (syntax highlighting) style to use. +pygments_style = 'sphinx' + +# autoapi_type = 'python' +# autoapi_dirs = [os.path.abspath('../../quanfima')] + + +# -- Options for HTML output ------------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +# +html_theme = 'sphinx_rtd_theme' + +# Theme options are theme-specific and customize the look and feel of a theme +# further. For a list of options available for each theme, see the +# documentation. +# +# html_theme_options = {} + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ['_static'] + +# Custom sidebar templates, must be a dictionary that maps document names +# to template names. +# +# The default sidebars (for documents that don't match any pattern) are +# defined by theme itself. Builtin themes are using these templates by +# default: ``['localtoc.html', 'relations.html', 'sourcelink.html', +# 'searchbox.html']``. +# +# html_sidebars = {} + + +# -- Options for HTMLHelp output --------------------------------------------- + +# Output file base name for HTML help builder. +htmlhelp_basename = 'quanfimadoc' + + +# -- Options for LaTeX output ------------------------------------------------ + +latex_elements = { + # The paper size ('letterpaper' or 'a4paper'). + # + # 'papersize': 'letterpaper', + + # The font size ('10pt', '11pt' or '12pt'). + # + # 'pointsize': '10pt', + + # Additional stuff for the LaTeX preamble. + # + # 'preamble': '', + + # Latex figure (float) alignment + # + # 'figure_align': 'htbp', +} + +# Grouping the document tree into LaTeX files. List of tuples +# (source start file, target name, title, +# author, documentclass [howto, manual, or own class]). +latex_documents = [ + (master_doc, 'quanfima.tex', u'quanfima Documentation', + u'Roman Shkarin, Andrei Shkarin', 'manual'), +] + + +# -- Options for manual page output ------------------------------------------ + +# One entry per manual page. List of tuples +# (source start file, name, description, authors, manual section). +man_pages = [ + (master_doc, 'quanfima', u'quanfima Documentation', + [author], 1) +] + + +# -- Options for Texinfo output ---------------------------------------------- + +# Grouping the document tree into Texinfo files. List of tuples +# (source start file, target name, title, author, +# dir menu entry, description, category) +texinfo_documents = [ + (master_doc, 'quanfima', u'quanfima Documentation', + author, 'quanfima', 'One line description of project.', + 'Miscellaneous'), +] + + +# -- Extension configuration ------------------------------------------------- diff --git a/quanfima/docs/source/index.rst b/quanfima/docs/source/index.rst new file mode 100644 index 0000000..9d4626b --- /dev/null +++ b/quanfima/docs/source/index.rst @@ -0,0 +1,29 @@ +.. title:: Quanfima + +====================== +Quanfima Documentation +====================== + +*Quanfima* (**qu**\antitative **an**\alysis of **fi**\brous **ma**\terials) +is a collection of useful routines for morphological analysis and visualization +of 2D/3D data from various areas of material science. The aim is to simplify +the analysis process by providing functionality for frequently required tasks +in the same place. + +- Analysis of fibrous structures by tensor-based method in 2D / 3D datasets. +- Estimation of structure diameters in 2D / 3D by a ray-casting method. +- Counting of particles in 2D / 3D datasets and providing a detailed report in + pandas.DataFrame format. +- Calculation of porosity measure for each material in 2D / 3D datasets. +- Visualization in 2D / 3D using matplotlib, visvis packages. + +User Guide +========== + +.. toctree:: + :maxdepth: 2 + + overview + quickstart + quanfima + change_log diff --git a/quanfima/docs/source/overview.rst b/quanfima/docs/source/overview.rst new file mode 100644 index 0000000..d17cde9 --- /dev/null +++ b/quanfima/docs/source/overview.rst @@ -0,0 +1,28 @@ +======== +Overview +======== + +Requirements +============ + +- Python 2.7 +- PyCUDA >= 2017.1.1 +- Numpy >= 1.13 +- SciPy >= 0.19 +- scikit-image >= 0.12 +- scikit-learn >= 0.18 +- Pandas >= 0.19 +- Matplotlib >= 2.0 + +Installation +============ + +The easiest way to install the latest version is by using pip:: + + $ pip install quanfima + +You may also use Git to clone the repository and install it manually:: + + $ git clone https://github.com/rshkarin/quanfima.git + $ cd quanfima + $ python setup.py install diff --git a/quanfima/docs/source/quanfima.rst b/quanfima/docs/source/quanfima.rst new file mode 100644 index 0000000..08e2089 --- /dev/null +++ b/quanfima/docs/source/quanfima.rst @@ -0,0 +1,35 @@ +============= +API Reference +============= + +quanfima\.morphology module +--------------------------- + +.. automodule:: quanfima.morphology + :members: + :undoc-members: + :show-inheritance: + +quanfima\.simulation module +--------------------------- + +.. automodule:: quanfima.simulation + :members: + :undoc-members: + :show-inheritance: + +quanfima\.utils module +---------------------- + +.. automodule:: quanfima.utils + :members: + :undoc-members: + :show-inheritance: + +quanfima\.visualization module +------------------------------ + +.. automodule:: quanfima.visualization + :members: + :undoc-members: + :show-inheritance: diff --git a/quanfima/docs/source/quickstart.rst b/quanfima/docs/source/quickstart.rst new file mode 100644 index 0000000..e4e61cf --- /dev/null +++ b/quanfima/docs/source/quickstart.rst @@ -0,0 +1,266 @@ +========== +Quickstart +========== + +Analysis of fibrous 2D data +^^^^^^^^^^^^^^^^^^^^^^^^^^^ +Open a grayscale image, perform segmentation, estimate porosity, analyze fiber +orientation and diameters, and plot the results. + +.. code-block:: python + + import numpy as np + from skimage import io, filters + from quanfima import morphology as mrph + from quanfima import visualization as vis + from quanfima import utils + + img = io.imread('../../data/polymer_slice.tif') + + th_val = filters.threshold_otsu(img) + img_seg = (img > th_val).astype(np.uint8) + + # estiamte porosity + pr = mrph.calc_porosity(img_seg) + for k,v in pr.items(): + print 'Porosity ({}): {}'.format(k, v) + + # prepare data and analyze fibers + data, skeleton, skeleton_thick = utils.prepare_data(img_seg) + cskel, fskel, omap, dmap, ovals, dvals = \ + mrph.estimate_fiber_properties(data, skeleton) + + # plot results + vis.plot_orientation_map(omap, fskel, min_label=u'0°', max_label=u'180°', + figsize=(10,10), + name='2d_polymer', + output_dir='../../data/results') + vis.plot_diameter_map(dmap, cskel, figsize=(10,10), cmap='gist_rainbow', + name='2d_polymer', + output_dir='../../data/results') + +.. code-block:: python + + >> Porosity (Material 1): 0.845488888889 + +.. image:: _static\2d_polymer_data.png + :width: 90 % + :align: center +.. image:: _static\2d_polymer_orientation_map.png + :width: 90 % + :align: center +.. image:: _static\2d_polymer_diameter_map.png + :width: 90 % + :align: center + +Analysis of 3D data of fibrous material +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +Open a micro-CT dataset, perform slice-wise segmentation, estimate porosity, +analyze 3D fiber orientation and diameters, and visualize the results. + +.. code:: python + + import numpy as np + from skimage import filters + from quanfima import morphology as mrph + from quanfima import visualization as vis + from quanfima import utils + + data = np.memmap('../../data/polymer3d_8bit_128x128x128.raw', + shape=(128,128,128), dtype=np.uint8, mode='r') + + data_seg = np.zeros_like(data, dtype=np.uint8) + for i in xrange(data_seg.shape[0]): + th_val = filters.threshold_otsu(data[i]) + data_seg[i] = (data[i] > th_val).astype(np.uint8) + + # estimate porosity + pr = mrph.calc_porosity(data_seg) + for k,v in pr.items(): + print 'Porosity ({}): {}'.format(k, v) + + # prepare data and analyze fibers + pdata, pskel, pskel_thick = utils.prepare_data(data_seg) + oprops = mrph.estimate_tensor_parallel('polymer_orientation_w32', pskel, + pskel_thick, 32, + '../../data/results') + + odata = np.load(oprops['output_path']).item() + lat, azth, skel = odata['lat'], odata['azth'], odata['skeleton'] + + dprops = mrph.estimate_diameter_single_run('polymer_diameter', + '../../data/results', + pdata, skel, lat, azth) + dmtr = np.load(dprops['output_path']).item()['diameter'] + + # plot results + vis.plot_3d_orientation_map('polymer_w32', lat, azth, + output_dir='../../data/results', + camera_azth=40.47, + camera_elev=32.5, + camera_fov=35.0, + camera_loc=(40.85, 46.32, 28.85), + camera_zoom=0.005124) + + vis.plot_3d_diameter_map('polymer_w32', dmtr, + output_dir='../../data/results', + measure_quantity='vox', + camera_azth=40.47, + camera_elev=32.5, + camera_fov=35.0, + camera_loc=(40.85, 46.32, 28.85), + camera_zoom=0.005124, + cb_x_offset=5, + width=620) + +.. code:: python + + >> Porosity (Material 1): 0.855631351471 + +.. image:: _static\polymer_w32_3d_orientation.png + :width: 78 % + :align: center +.. image:: _static\polymer_w32_3d_diameter.png + :width: 80 % + :align: center + +Estimation of p-values +^^^^^^^^^^^^^^^^^^^^^^ +Estimate p-values between several groups of samples with corresponding +measurements of some material's property. + +.. code:: python + + import pandas as pd + from quanfima import utils + + prop_vals = [174.93, 182.42, 194.61, 234.6, 229.73, 242.6, 38.78, 37.79, + 32.06, 14.81, 15.23, 13.84] + mat_groups = ['PCL_cl', 'PCL_cl', 'PCL_cl', 'PCL_wa', 'PCL_wa', 'PCL_wa', + 'PCL_SiHA_cl', 'PCL_SiHA_cl', 'PCL_SiHA_cl', 'PCL_SiHA_wa', + 'PCL_SiHA_wa', 'PCL_SiHA_wa'] + df_elongation = pd.DataFrame({'elongation': prop_vals, 'type': mat_groups}) + + _, _ = utils.calculate_tukey_posthoc(df_elongation, 'elongation', + name='samples_elongation', + write=True, + output_dir='../../data/results') + +.. code:: python + + >> Tukey post-hoc (elongation) + >> Multiple Comparison of Means - Tukey HSD,FWER=0.05 + >> ========================================================= + >> group1 group2 meandiff lower upper reject + >> --------------------------------------------------------- + >> PCL_SiHA_cl PCL_SiHA_wa -21.5833 -37.8384 -5.3282 True + >> PCL_SiHA_cl PCL_cl 147.7767 131.5216 164.0318 True + >> PCL_SiHA_cl PCL_wa 199.4333 183.1782 215.6884 True + >> PCL_SiHA_wa PCL_cl 169.36 153.1049 185.6151 True + >> PCL_SiHA_wa PCL_wa 221.0167 204.7616 237.2718 True + >> PCL_cl PCL_wa 51.6567 35.4016 67.9118 True + >> --------------------------------------------------------- + >> ['PCL_SiHA_cl' 'PCL_SiHA_wa' 'PCL_cl' 'PCL_wa'] + >> PCL_SiHA_cl-PCL_SiHA_wa : 0.011919282004 + >> PCL_SiHA_wa-PCL_cl : 0.001 + >> PCL_SiHA_wa-PCL_wa : 0.001 + >> PCL_cl-PCL_wa : 0.001 + >> PCL_SiHA_cl-PCL_wa : 0.001 + >> PCL_SiHA_cl-PCL_cl : 0.001 + +Simulate and count particles in 3D data +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +Count and estimate properties of particles in a generated dataset comprised +of spheres of varying radius. + +.. code:: python + + from quanfima import simulation + from quanfima import morphology as mrph + + volume, diameter, _, _ = simulation.simulate_particles((512,512,512), + n_particles=1000) + + stats, labeled_volume = mrph.object_counter(volume) + +.. code:: python + + >> Label Area Perimeter Sphericity + >> 0 1.0 20479.0 2896.958728 1.249533 + >> 1 2.0 5575.0 1184.028571 1.284158 + >> 2 3.0 57777.0 5816.142853 1.242660 + >> 3 4.0 17077.0 2545.194226 1.260001 + >> 4 5.0 5575.0 1184.028571 1.284158 + >> 5 6.0 65267.0 6348.926691 1.234752 + >> .. ... ... ... ... + >> 791 792.0 2109.0 605.185858 1.314154 + >> 792 793.0 257.0 134.225397 1.456369 + >> 793 794.0 257.0 134.225397 1.456369 + >> 794 795.0 123.0 78.627417 1.521179 + + >> [795 rows x 4 columns] + +Simulate fibers and estimate properties +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +Simulate a 3D dataset containing some number of fibers, estimate their +properties and visualize. + +.. code:: python + + import numpy as np + from scipy import ndimage as ndi + from skimage import morphology + from quanfima import simulation + from quanfima import morphology as mrph + from quanfima import utils + from quanfima import visualization as vis + + volume, lat_ref, azth_ref, diameter, _, _ = \ + simulation.simulate_fibers((128,128,128), n_fibers=30, max_fails=100, + radius_lim=(2, 3), gap_lim=(3,5)) + volume = volume.astype(np.uint8) + volume = ndi.binary_fill_holes(volume) + volume = ndi.median_filter(volume, footprint=morphology.ball(2)) + lat_ref = ndi.median_filter(lat_ref, footprint=morphology.ball(2)) + azth_ref = ndi.median_filter(azth_ref, footprint=morphology.ball(2)) + + # prepare data and analyze fibers + pdata, pskel, pskel_thick = utils.prepare_data(volume) + oprops = mrph.estimate_tensor_parallel('dataset_orientation_w36', + pskel, pskel_thick, 36, + '../../data/results') + + odata = np.load(oprops['output_path']).item() + lat, azth, skel = odata['lat'], odata['azth'], odata['skeleton'] + + dprops = mrph.estimate_diameter_single_run('dataset_diameter', + '../../data/results', + pdata, skel, lat, azth) + dmtr = np.load(dprops['output_path']).item()['diameter'] + + # plot results + vis.plot_3d_orientation_map('dataset_w36', lat, azth, + output_dir='../../data/results', + camera_azth=40.47, + camera_elev=32.5, + camera_fov=35.0, + camera_loc=(40.85, 46.32, 28.85), + camera_zoom=0.005124) + + vis.plot_3d_diameter_map('dataset_w36', dmtr, + output_dir='../../data/results', + measure_quantity='vox', + camera_azth=40.47, + camera_elev=32.5, + camera_fov=35.0, + camera_loc=(40.85, 46.32, 28.85), + camera_zoom=0.005124, + cb_x_offset=5, + width=620) + +.. image:: _static\dataset_w36_3d_orientation.png + :width: 78 % + :align: center +.. image:: _static\dataset_w36_3d_diameter.png + :width: 80 % + :align: center diff --git a/quanfima/docs/source/static/2d_polymer_data.png b/quanfima/docs/source/static/2d_polymer_data.png new file mode 100644 index 0000000..f4883fc Binary files /dev/null and b/quanfima/docs/source/static/2d_polymer_data.png differ diff --git a/quanfima/docs/source/static/2d_polymer_diameter_map.png b/quanfima/docs/source/static/2d_polymer_diameter_map.png new file mode 100644 index 0000000..2a814a5 Binary files /dev/null and b/quanfima/docs/source/static/2d_polymer_diameter_map.png differ diff --git a/quanfima/docs/source/static/2d_polymer_diameter_map_600px.png b/quanfima/docs/source/static/2d_polymer_diameter_map_600px.png new file mode 100644 index 0000000..4401f3b Binary files /dev/null and b/quanfima/docs/source/static/2d_polymer_diameter_map_600px.png differ diff --git a/quanfima/docs/source/static/2d_polymer_orientation_map.png b/quanfima/docs/source/static/2d_polymer_orientation_map.png new file mode 100644 index 0000000..9e01e85 Binary files /dev/null and b/quanfima/docs/source/static/2d_polymer_orientation_map.png differ diff --git a/quanfima/docs/source/static/2d_polymer_orientation_map_600px.png b/quanfima/docs/source/static/2d_polymer_orientation_map_600px.png new file mode 100644 index 0000000..7175eb2 Binary files /dev/null and b/quanfima/docs/source/static/2d_polymer_orientation_map_600px.png differ diff --git a/quanfima/docs/source/static/dataset_w36_3d_diameter.png b/quanfima/docs/source/static/dataset_w36_3d_diameter.png new file mode 100644 index 0000000..5ffadd2 Binary files /dev/null and b/quanfima/docs/source/static/dataset_w36_3d_diameter.png differ diff --git a/quanfima/docs/source/static/dataset_w36_3d_orientation.png b/quanfima/docs/source/static/dataset_w36_3d_orientation.png new file mode 100644 index 0000000..4f3f688 Binary files /dev/null and b/quanfima/docs/source/static/dataset_w36_3d_orientation.png differ diff --git a/quanfima/docs/source/static/dataset_w36_color_bar.png b/quanfima/docs/source/static/dataset_w36_color_bar.png new file mode 100644 index 0000000..9aa094c Binary files /dev/null and b/quanfima/docs/source/static/dataset_w36_color_bar.png differ diff --git a/quanfima/docs/source/static/logo.png b/quanfima/docs/source/static/logo.png new file mode 100644 index 0000000..7339957 Binary files /dev/null and b/quanfima/docs/source/static/logo.png differ diff --git a/quanfima/docs/source/static/polymer_w32_3d_diameter.png b/quanfima/docs/source/static/polymer_w32_3d_diameter.png new file mode 100644 index 0000000..8513db5 Binary files /dev/null and b/quanfima/docs/source/static/polymer_w32_3d_diameter.png differ diff --git a/quanfima/docs/source/static/polymer_w32_3d_orientation.png b/quanfima/docs/source/static/polymer_w32_3d_orientation.png new file mode 100644 index 0000000..ea12607 Binary files /dev/null and b/quanfima/docs/source/static/polymer_w32_3d_orientation.png differ diff --git a/quanfima/docs/source/static/polymer_w32_color_bar.png b/quanfima/docs/source/static/polymer_w32_color_bar.png new file mode 100644 index 0000000..9aa094c Binary files /dev/null and b/quanfima/docs/source/static/polymer_w32_color_bar.png differ diff --git a/quanfima/docs/source/static/test.txt b/quanfima/docs/source/static/test.txt new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/quanfima/docs/source/static/test.txt @@ -0,0 +1 @@ + diff --git a/quanfima/docs/test.txt b/quanfima/docs/test.txt new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/quanfima/docs/test.txt @@ -0,0 +1 @@ + diff --git a/quanfima/quanfima/__init__.py b/quanfima/quanfima/__init__.py new file mode 100644 index 0000000..7a3aba5 --- /dev/null +++ b/quanfima/quanfima/__init__.py @@ -0,0 +1,19 @@ +"""Quanfima initialization.""" + +import logging + +__version__ = '0.1a' +__log_name__ = 'matscipy' +log = logging.getLogger(__log_name__) + +cuda_available = True +try: + import pycuda +except ImportError: + cuda_available = False + +visvis_available = True +try: + import visvis +except ImportError: + visvis_available = False diff --git a/quanfima/quanfima/morphology.py b/quanfima/quanfima/morphology.py new file mode 100644 index 0000000..9d383b2 --- /dev/null +++ b/quanfima/quanfima/morphology.py @@ -0,0 +1,1159 @@ +from __future__ import print_function +import math +import time +import os +import itertools +import numpy as np +from skimage import feature, measure, filters +from skimage.util.shape import view_as_blocks +from scipy import ndimage as ndi +from scipy.spatial import distance +import vigra +import pandas as pd +from multiprocessing import Pool +from quanfima.quanfima import cuda_available + +if cuda_available: + import pycuda.autoinit + import pycuda.driver as cuda + from pycuda.compiler import SourceModule + from pycuda import gpuarray + +# Constants of object counter +_MEASUREMENTS = { + 'Label': 'label', + 'Area': 'area', + 'Perimeter': 'perimeter' +} + + +def cast_ray(theta, y0, x0, fiber_mask, ray_len=100): + """Computes a distance between two detected points at `fiber_mask`. + + Casts a ray from the point (x0, y0) at `fiber_mask` under angle `theta` of maximum + length of `ray_len`. The ray is casted towards to opposite directions, thus two points + (x1,y1) and (x2,y2) are detected at opposite borders of fiber. Finally, the Euclidean + distance between these points is calculated and returned. + + Parameters + ---------- + theta : float + Indicates the angle under which the ray is casted. + + y0 : integer + Indicates the y-axis component of the origin of ray emission. + + x0 : integer + Indicates the x-axis component of the origin of ray emission. + + fiber_mask : 2D array + Indicates the binary array of fibers. + + ray_len : integer + Indicates the maximum length of the ray. + + Returns + ------- + distance : float + The distance between (x1,y0) and (x2,y2) detected points. + """ + m = math.tan(theta) + y1, x1, y2, x2 = 0, 0, 0, 0 + + for i in range(-ray_len, 0)[::-1]: + if np.abs(m) <= 1.: + y = int(y0 - i * m) + x = int(x0 + i) + else: + y = int(y0 - i) + x = int(x0 + i * 1./m) + + if x < 0 or x >= fiber_mask.shape[1] or \ + y < 0 or y >= fiber_mask.shape[0]: + continue + + if fiber_mask[y, x] == 0: + break + + y1, x1 = y, x + + for i in range(ray_len): + if np.abs(m) <= 1.: + y = int(y0 - i * m) + x = int(x0 + i) + else: + y = int(y0 - i) + x = int(x0 + i * 1./m) + + if x < 0 or x >= fiber_mask.shape[1] or \ + y < 0 or y >= fiber_mask.shape[0]: + continue + + if fiber_mask[y, x] == 0: + break + + y2, x2 = y, x + + return distance.euclidean((y1, x1), (y2, x2)) + + +def scan_fiber_thickness(angle, patch, angular_step=1., tilt_range=range(-5, 6), ray_len=100): + """Computes an average diameter of a structure within a patch. + + Estimates an average distance from sequence of distances calculated with `cast_ray` + function for a range of angles `tilt_range`. + + Parameters + ---------- + orientation : float + Indicates the orientation angle of a structure centered within a patch. + + patch : 2D array + Indicates the patch with a centered structure having orientation `angle`. + + angular_step : float + Indicates the angular step of the scanning range in degrees. + + tilt_range : integer + Indicates the scanning range of steps to estimate the distance. + + ray_len : integer + Indicates the maximum length of the ray. + + Returns + ------- + average_diameter : float + The average diameter of a structure centered within a patch. + """ + thickness_dist = [] + + y0, x0 = patch.shape[0]/2, patch.shape[1]/2 + pi2 = np.pi/2. + rad_p_deg = np.deg2rad(angular_step) + + for tilt_amount in tilt_range: + theta = angle + pi2 + tilt_amount * rad_p_deg + thickness_dist.append(cast_ray(theta, y0, x0, patch, ray_len=ray_len)) + + average_distance = np.mean(thickness_dist) + return average_distance + + +def estimate_fiber_properties(fiber_mask, fiber_skel, paddding=25, window_radius=12, + orient_type='tensor', diameter_window_radius=12): + """Computes orientation and diameter of fibers at every point of a skeleton. + + Estimates orientation and diameter of fibers at every point of a skeleton. The orientation + is estimated using either tensor-based or PCA-based approach. The distance is evaluated + by scanning a mask of fibers in a speficied angular range. + + Parameters + ---------- + fiber_mask : 2D array + Indicates the binary data of fibers produced by the segmentation process. + + fiber_skel : 2D array + Indicates the skeleton produced by thinning of the binary data of fibers. + + paddding : integer + Indicates the amount of padding at the corners to prevent an estimation error. + + window_radius : integer + Indicates the radius of the local window of orientation calculation, which leads to + patches of size (`window_radius`*2+1) x (`window_radius`*2+1). + + orient_type : str + Indicates the type of algorithm for orientation estimation ('tensor' or 'pca'). + + diameter_window_radius : integer + Indicates the radius of the local window of diameter estimation, which leads to + patches of size (`diameter_window_radius`*2+1) x (`diameter_window_radius`*2+1). + + Returns + ------- + (clear_fiber_skel, fiber_skel, output_orientation_map, output_diameter_map, + orientation_vals, diameter_vals) : tuple of arrays + The skeleton with removed intersections, the skeleton, the orientation map, the + diameter map, the arrays of orientation and diameter values. + """ + orientation_vals, diameter_vals = [], [] + + padded_fiber_skel = np.pad(fiber_skel, pad_width=(paddding,), mode='constant', constant_values=0) + padded_fiber_mask = np.pad(fiber_mask, pad_width=(paddding,), mode='constant', constant_values=0) + + output_orientation_map, output_diameter_map = \ + np.zeros_like(padded_fiber_skel, dtype=np.float32), \ + np.zeros_like(padded_fiber_skel, dtype=np.float32) + + pcy, pcx = (window_radius*2+1)/2., (window_radius*2+1)/2. + + ycords, xcords = np.nonzero(padded_fiber_skel) + + clear_padded_fiber_skel = padded_fiber_skel.copy() + + method = feature.corner_harris(clear_padded_fiber_skel, sigma=1.5) + corner_points = feature.corner_peaks(method, min_distance=3) + + for i, (yy, xx) in enumerate(corner_points): + clear_padded_fiber_skel[yy-1:yy+2, xx-1:xx+2] = np.zeros((3, 3), + dtype=clear_padded_fiber_skel.dtype) + + # PCA-based approach + def calc_orientation_pca(patch): + lbls = ndi.label(patch, structure=ndi.generate_binary_structure(2, 2))[0] + + rgns = measure.regionprops(lbls) + dists = [(r.label, distance.euclidean((pcy, pcx), (r.centroid[0], r.centroid[1]))) for r in rgns] + dists = sorted(dists, key=lambda x: x[1]) + orintations = [r.orientation for r in rgns if r.label == dists[0][0]] + centroids = [r.centroid for r in rgns if r.label == dists[0][0]] + + orientation = orintations[0] + y0, x0 = centroids[0] + + return orientation, y0, x0 + + # Tensor-based approach + def calc_orientation_tensor(patch): + Axx, Axy, Ayy = feature.structure_tensor(patch, sigma=0.1) + tensor_vals = np.array([[np.mean(Axx), np.mean(Axy)], [np.mean(Axy), np.mean(Ayy)]]) + + w, v = np.linalg.eig(tensor_vals) + orientation = math.atan2(*v[:,np.argmax(w)]) + y0, x0 = patch.shape[0]/2, patch.shape[1]/2 + + return orientation, y0, x0 + + for i, (yy, xx) in enumerate(zip(ycords, xcords)): + patch_skel = padded_fiber_skel[(yy-window_radius):(yy+window_radius+1), + (xx-window_radius):(xx+window_radius+1)] + + patch_mask = padded_fiber_mask[(yy-window_radius):(yy+window_radius+1), + (xx-window_radius):(xx+window_radius+1)] + + if orient_type == 'tensor': + orientation, y0, x0 = calc_orientation_tensor(patch_skel) + elif orient_type == 'pca': + orientation, y0, x0 = calc_orientation_pca(patch_skel) + + final_orientation = orientation if orientation >= 0.0 else (np.pi - np.abs(orientation)) + orientation_vals.append(final_orientation) + output_orientation_map[yy, xx] = final_orientation + + if clear_padded_fiber_skel[yy, xx] != 0: + gy, gx = yy - window_radius + y0, xx - window_radius + x0 + + fiber_tpatch = padded_fiber_mask[yy-diameter_window_radius:yy+diameter_window_radius+1, + xx-diameter_window_radius:xx+diameter_window_radius+1] + + final_thickness = scan_fiber_thickness(orientation, patch_mask) + output_diameter_map[yy, xx] = final_thickness + diameter_vals.append(final_thickness) + + output_orientation_map = output_orientation_map[paddding:-paddding, paddding:-paddding] + output_diameter_map = output_diameter_map[paddding:-paddding, paddding:-paddding] + clear_fiber_skel = clear_padded_fiber_skel[paddding:-paddding, paddding:-paddding] + + return (clear_fiber_skel, fiber_skel, + output_orientation_map, output_diameter_map, + orientation_vals, diameter_vals) + +def estimate_fourier_orientation(data, grid_shape=(2,2), sigma=2., zoom=1., order=3): + """Computes orientation at every block of the subdivided image. + + Subdivides the image into the grid of blocks `grid_shape`, it computes + 2D FFT for every block, then the real components are segmented by + the Otsu thresholding and the PCA-based approach calculates the orientation + of structures within each block. + + Parameters + ---------- + data : ndarray + Indicates the grayscale 2D image. + + grid_shape : tuple + Indicates the number of blocks to subdivide the image. + + sigma : float + Indicates the sigma value of the Gaussian filter to smooth the real part + of 2D Fourier spectrum. + + zoom : float + Indicates the upscaling factor of each block before applying 2D FFT. + + order : str + Indicates the order of interpolation used in the upscaling procedure. + + Returns + ------- + (orient_blocks, block_shape) : tuple of arrays + The 2D array of the orientation angle within each block, + the shape of a block. + """ + block_shape = tuple([int(math.floor(d/float(gs))) for d, gs in + zip(data.shape, grid_shape)]) + data_blocks = view_as_blocks(data, block_shape=block_shape) + orient_blocks = np.zeros(grid_shape, dtype=np.float32) + + for i in xrange(grid_shape[0]): + for j in xrange(grid_shape[1]): + y0, x0 = block_shape[0] * j + block_shape[0]/2, \ + block_shape[1] * i + block_shape[1]/2 + dblock = data_blocks[i, j] + + if zoom != 1.: + dblock = ndi.interpolation.zoom(dblock, zoom, order=order) + + dblock_freq = np.abs(np.fft.fftshift(np.fft.fft2(dblock)).real) + dblock_freq = ndi.gaussian_filter(dblock_freq, sigma=sigma) + dblock_mask = (dblock_freq > filters.threshold_otsu(dblock_freq)).astype(np.uint8) + + lbls = ndi.label(dblock_mask, + structure=ndi.generate_binary_structure(2,2))[0] + rgns = measure.regionprops(lbls) + + dists = [(r.label, r.area) for r in rgns] + dists = sorted(dists, key=lambda x: x[1]) + + orient_blocks[i, j] = [r.orientation for r in rgns + if r.label == dists[0][0]][0] + + return (orient_blocks, block_shape) + + +def numpy3d_to_array(np_array, allow_surface_bind=True): + """Converts 3D numpy array to 3D device array. + """ + # numpy3d_to_array + # taken from pycuda mailing list (striped for C ordering only) + + d, h, w = np_array.shape + + descr = cuda.ArrayDescriptor3D() + descr.width = w + descr.height = h + descr.depth = d + descr.format = cuda.dtype_to_array_format(np_array.dtype) + descr.num_channels = 1 + descr.flags = 0 + + if allow_surface_bind: + descr.flags = cuda.array3d_flags.SURFACE_LDST + + device_array = cuda.Array(descr) + + copy = cuda.Memcpy3D() + copy.set_src_host(np_array) + copy.set_dst_array(device_array) + copy.width_in_bytes = copy.src_pitch = np_array.strides[1] + copy.src_height = copy.height = h + copy.depth = d + + copy() + + return device_array + + +def _filter_coords(data, ws, th_nval=0.8): + """Filters coordinates of 3D local windows by number of non-zero values speficied by `th_nval`. + """ + ws2 = ws / 2 + + def check_point(data, pt, ws2, th): + size = max(data.shape) + lim0 = pt - ws2 + lim1 = pt + ws2 + + if any(np.array(lim0) < 0) or any(np.array(lim1) > (size - 1)): + return False + + z0, y0, x0 = lim0 + z1, y1, x1 = lim1 + + vol = data[z0:z1, y0:y1, x0:x1] + nz_vals = np.count_nonzero(vol) + + if np.count_nonzero(vol) < ws2*2.*th: + return False + + return True + + filtered_coords = [pt for pt in np.transpose(np.nonzero(data)) if check_point(data, pt, ws2, th_nval)] + + return np.transpose(filtered_coords) + + +def extract_patch(data, pt, ws2): + """Extract a path from `data` of radius `ws2`. + + Parameters + ---------- + data : 3D array + Indicates the 3D data array. + + pt : array or tuple + Indicates the coordinates of a point within `data`. + + ws2 : integer + Indicates a half of the window size of the patch to be extracted. + + Returns + ------- + patch : 3D array + The 3D patch extracted from the local window around the point `pt` in the `data` + with radius of `ws2`. + """ + lim0 = pt - ws2 + lim1 = pt + ws2 + + if any(np.array(lim0) < 0) or any(np.array(lim1) > (data.shape[0] - 1)): + return None + + z0, y0, x0 = lim0 + z1, y1, x1 = lim1 + patch = data[z0:z1, y0:y1, x0:x1] + return patch + + +def orientation_3d_tensor_vigra(data, sigma=0.1): + """Computes 3D orientation from a 3D structure tensor of `data`. + + Parameters + ---------- + data : ndarray + Indicates the N-dimensional array. + + sigma : float + Indicates the sigma value of the Gaussian filter. + + Returns + ------- + (lat, azth) : tuple of floats + The latitude / elevation and azimuth component of 3D orientation of structures within + the patch `data`. + """ + img = vigra.filters.structureTensor(data, 1, 1, sigma_d=sigma) + Axx = img[:, :, :, 0] + Axy = img[:, :, :, 1] + Axz = img[:, :, :, 2] + Ayy = img[:, :, :, 3] + Ayz = img[:, :, :, 4] + Azz = img[:, :, :, 5] + + tensor_vals = np.array([[np.mean(Azz), np.mean(Ayz), np.mean(Axz)], + [np.mean(Ayz), np.mean(Ayy), np.mean(Axy)], + [np.mean(Axz), np.mean(Axy), np.mean(Axx)]]) + + tensor_vals = tensor_vals[::-1, ::-1] + + eps = 1e-8 + w, v = np.linalg.eig(tensor_vals) + + mv = v[:, np.argmin(w)] # z, y, x + mv[np.abs(mv) < eps] = 0 + + G = np.sqrt(mv[2]**2 + mv[1]**2) + lat = np.arcsin(np.around(G, decimals=3)) + azth = np.arctan(mv[1] / mv[2]) if mv[2] else np.pi/2. + + return (lat, azth) + + +def estimate_tensor(name, skel, data, window_size, output_dir, sigma=0.025, make_output=True): + """Computes 3D orientation at every point of a skeleton of data. + + Estimates 3D orientation at every point of the skeleton `skel` extracted from the binary + data `data` within a 3D local window of size `window_size` using the tensor-based approach + with a Gaussian smoothing of `sigma`. + + Parameters + ---------- + name : str + Indicates the name of the output npy file. + + skel : 3D array + Indicates the skeleton of the binary data. + + data : 3D array + Indicates the 3D binary data. + + window_size : integer + Indicates the size of the 3D local window. + + output_dir : str + Indicates the path to the output folder where the data will be stored. + + output_fmt : str + Indicates the format of + + sigma : float + Indicates the sigma value of the Gaussian filter. + + make_output : boolean + Specifies if the estimated data should be stored. + + Returns + ------- + output_props : dict + The dictionary of properties specifying the sample name, the algorithm name, the number + of processes, and the execution time. + """ + output_props = dict() + + Z, Y, X = skel.nonzero() + tens_lat_arr = np.zeros_like(skel, dtype=np.float32) + tens_azth_arr = np.zeros_like(skel, dtype=np.float32) + skel_est = np.zeros_like(skel, dtype=np.int32) + + ws = np.uint32(window_size) + ws2 = ws/2 + skel_shape = skel.shape + + output_props['sample_name'] = name + output_props['type'] = 'tensor' + output_props['n_processes'] = 1 + + ts = time.time() + + for idx, pt in enumerate(zip(Z, Y, X)): + lim0 = pt - ws2 + lim1 = pt + ws2 + + if any(np.array(lim0) < 0) or any(np.array(lim1) > (skel_shape[0] - 1)): + skel_est[pt] = -1 + continue + + z0, y0, x0 = lim0 + z1, y1, x1 = lim1 + + area = data[z0:z1, y0:y1, x0:x1] + + lat, azth = orientation_3d_tensor_vigra(area, sigma) + + tens_lat_arr[pt] = lat + tens_azth_arr[pt] = azth + skel_est[pt] = 255 + + te = time.time() + output_props['time'] = te-ts + + print("Tensor time: {}s" % (output_props['time'])) + + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + if make_output: + opath = os.path.join(output_dir, '{}.npy').format(name) + output_props['output_path'] = opath + output = {'lat': tens_lat_arr, 'azth': tens_azth_arr, 'skeleton': skel_est, 'indices': skel_est > 0} + output['props'] = output_props + np.save(opath, output) + + return output_props + + +def execute_tensor(patch, sigma): + """Executes the tensor-based approach for a `patch` and `sigma` for Gaussian smoothing. + + Parameters + ---------- + patch : 3D aray + Indicates the data patch of some size. + + sigma : float + Indicates the sigma value of the Gaussian filter. + + Returns + ------- + (lat, azth, 255) : tuple + The latitude / elevation and azimuth components of 3D orientation with valid value + for the point at the skeleton. If the patch is None, then return (0, 0, -1) where + -1 indicates that this point of the skeleton is invalid. + """ + if patch is None: + return (0, 0, -1) + + lat, azth = orientation_3d_tensor_vigra(patch, sigma) + + return (lat, azth, 255) + + +def unpack_execute_tensor(args): + """Unpack input arguments and return result of `execute_tensor` function + """ + return execute_tensor(*args) + + +def estimate_tensor_parallel(name, skel, data, window_size, output_dir, sigma=0.025, + make_output=True, n_processes=4): + """Computes 3D orientation at every point of a skeleton of data in parallel processes. + + Estimates 3D orientation at every point of the skeleton `skel` extracted from the binary + data `data` within a 3D local window of size `window_size` using the tensor-based approach + with a Gaussian smoothing of `sigma`. The orientation is estimated simultaneously at + `n_processes` parallel processes. + + Parameters + ---------- + name : str + Indicates the name of the output npy file. + + skel : 3D array + Indicates the skeleton of the binary data. + + data : 3D array + Indicates the 3D binary data. + + window_size : integer + Indicates the size of the 3D local window. + + output_dir : str + Indicates the path to the output folder where the data will be stored. + + sigma : float + Indicates the sigma value of the Gaussian filter. + + make_output : boolean + Specifies if the estimated data should be stored. + + n_processes : integer + Indicates the number of the parallel processes. + + Returns + ------- + output_props : dict + The dictionary of properties specifying the sample name, the algorithm name, the number + of processes, and the execution time. + """ + output_props = dict() + + Z, Y, X = skel.nonzero() + tens_lat_arr = np.zeros_like(skel, dtype=np.float32) + tens_azth_arr = np.zeros_like(skel, dtype=np.float32) + skel_est = np.zeros_like(skel, dtype=np.int32) + + ws = np.uint32(window_size) + ws2 = ws/2 + + output_props['sample_name'] = name + output_props['type'] = 'tensor_parallel' + output_props['n_processes'] = n_processes + + ts = time.time() + + pts = zip(Z, Y, X) + data_patches = [extract_patch(data, pt, ws2) for pt in pts] + print(len(data_patches)) + args = zip(data_patches, itertools.repeat(sigma)) + + proc_pool = Pool(processes=n_processes) + results = np.array(proc_pool.map(unpack_execute_tensor, args)) + proc_pool.close() + proc_pool.join() + proc_pool.terminate() + + te = time.time() + output_props['time'] = te-ts + + lat_arr, azth_arr, skel_arr = results.T + tens_lat_arr[Z, Y, X] = lat_arr + tens_azth_arr[Z, Y, X] = azth_arr + skel_est[Z, Y, X] = skel_arr + + print("Tensor parallel time: %fs" % (output_props['time'])) + + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + if make_output: + opath = os.path.join(output_dir, '{}.npy').format(name) + output_props['output_path'] = opath + output = {'lat': tens_lat_arr, 'azth': tens_azth_arr, 'skeleton': skel_est, 'indices': skel_est > 0} + output['props'] = output_props + np.save(opath, output) + + return output_props + + +def _diameter_kernel(): + """Returns the CUDA kernel to estimate 3D diameter of structures in a 3D local window. + """ + diameter_kernel_src = """ + + texture tex_data; + + __global__ void diameter3d (unsigned int width, + unsigned int height, + unsigned int depth, + const int n_points, + const float norm_factor, + const int max_iters, + const int n_scan_angles, + const int *X, + const int *Y, + const int *Z, + const float *scan_angl_arr, + const float *azth_data, + const float *lat_data, + float *radius_arr) + { + unsigned long blockId, idx; + blockId = blockIdx.x + blockIdx.y * gridDim.x; + idx = blockId * (blockDim.x * blockDim.y) + (threadIdx.y * blockDim.x) + threadIdx.x; + + if (idx > n_points) { + return; + } + + float _x, _y, _z; + _x = (float)X[idx] + 0.5; + _y = (float)Y[idx] + 0.5; + _z = (float)Z[idx] + 0.5; + + // ----------------------------------------- + // Find the diameter + + float azth = azth_data[idx]; + azth += M_PI_2; + float lat = lat_data[idx]; + + float cy = cosf(lat), sy = sinf(lat); + float cz = cosf(azth), sz = sinf(azth); + //float cz = -sinf(azth), sz = cosf(azth); // taking into account azth + pi/2 + + //vector along fiber + //float dx = -sy * sz; + //float dy = cy* sz; + //float dz = cy; + float uvec[3] = {0.0, 0.0, 1.0}; + float fiber_vector_x[3] = {0.0, + uvec[2]*cy + uvec[3]*sy, + -uvec[2]*sy + uvec[3]*cy}; + + float fiber_vector_z[3] = {fiber_vector_x[0]*cz - fiber_vector_x[1]*sz, + fiber_vector_x[0]*sz + fiber_vector_x[1]*cz, + fiber_vector_x[2]}; + + float dx = fiber_vector_z[0], dy = fiber_vector_z[1], dz = fiber_vector_z[2]; + + //scan vector perpendicular to a fiber vector (rotation X -> Z) + float scan_vec[3] = {0, 1, 0}; // unit vector perpendicular to default (0,0,1) direction + float rot_scan_vec_x[3] = {0.0, + scan_vec[1]*cy + scan_vec[2]*sy, + -scan_vec[1]*sy + scan_vec[2]*cy}; + + + float rot_scan_vec_z[3] = {rot_scan_vec_x[0]*cz - rot_scan_vec_x[1]*sz, + rot_scan_vec_x[0]*sz + rot_scan_vec_x[1]*cz, + rot_scan_vec_x[2]}; + + float out_radius = 0; + + for (int scan_angl_idx = 0; scan_angl_idx < n_scan_angles; scan_angl_idx++) { + float theta = scan_angl_arr[scan_angl_idx]; + + float ct = cosf(theta), st = sinf(theta); + float x = rot_scan_vec_z[0], y = rot_scan_vec_z[1], z = rot_scan_vec_z[2]; + float u = dx, v = dy, w = dz; + + //rotation of point (x,y,z) around axis (u,v,w) + float scan_vec_coords[3] = + {u*(u*x + v*y + w*z)*(1.0f - ct) + x*ct + (-w*y + v*z)*st, + v*(u*x + v*y + w*z)*(1.0f - ct) + y*ct + (w*x - u*z)*st, + w*(u*x + v*y + w*z)*(1.0f - ct) + z*ct + (-v*x + u*y)*st}; + + float nc[3] = {_x, _y, _z}; + float p[3]; + + for (int i = 0; i < max_iters; i++) { + nc[0] += scan_vec_coords[0]; + nc[1] += scan_vec_coords[1]; + nc[2] += scan_vec_coords[2]; + + if (tex3D(tex_data, nc[0], nc[1], nc[2]) == 0) { + p[0] = nc[0]; + p[1] = nc[1]; + p[2] = nc[2]; + break; + } + } + + out_radius += norm3df(p[0] - _x, p[1] - _y, p[2] - _z); + } + + radius_arr[idx] = out_radius * norm_factor; + } + """ + + dm_program = SourceModule(diameter_kernel_src) + diameter3d = dm_program.get_function("diameter3d") + + return dm_program, diameter3d + + +def estimate_diameter_gpu(skel, data, lat_data, azth_data, n_scan_angles, max_iters=150, do_reshape=True): + """Computes 3D diameter at every point of a skeleton of data using GPU. + + Estimates 3D diameter at every point of the skeleton `skel` extracted from the binary + data `data` with help of orientation information provided by `lat_data` and `azth_data` + arrays. The diameter is evaluated with a ray casting approach `cast_ray` adapted for + a 3D case. + + Parameters + ---------- + skel : 3D array + Indicates the skeleton of the binary data. + + data : 3D array + Indicates the 3D binary data. + + lat_data : 3D array + Indicates the 3D array containing latitude / elevation angle at every point of + the skeleton in radians. + + azth_data : 3D array + Indicates the 3D array containing azimuth angle at every point of the skeleton + in radians. + + n_scan_angles : int + Indicates the number of scanning angles on a range [0, 360] degrees. + + max_iters : int + Indicates the maximum length of a ray in each direction. + + do_reshape : boolean + Specifies if the output array should be reshaped immediately after estimation. + + Returns + ------- + out : dict + The dictionary of the 3D array of estimated diameter and the execution time. + """ + if not cuda_available: + print('The pycuda package is not found. The diameter estimation cannot be done.') + return None + + program, diameter3d = _diameter_kernel() + + Z, Y, X = np.int32(skel.nonzero()) + depth, height, width = np.uint32(skel.shape) + + scan_angl_arr = np.deg2rad(np.float32(np.linspace(0, 360, num=n_scan_angles, endpoint=False))) + radius_arr = np.zeros_like(Z, dtype=np.float32) + + lat_data_1d = lat_data[skel.nonzero()] + azth_data_1d = azth_data[skel.nonzero()] + + gpu_X = gpuarray.to_gpu(X) + gpu_Y = gpuarray.to_gpu(Y) + gpu_Z = gpuarray.to_gpu(Z) + + gpu_lat_data_1d = gpuarray.to_gpu(lat_data_1d) + gpu_azth_data_1d = gpuarray.to_gpu(azth_data_1d) + + gpu_radius_arr = gpuarray.to_gpu(radius_arr) + gpu_scan_angl_arr = gpuarray.to_gpu(scan_angl_arr) + + gpu_rad_tex = program.get_texref('tex_data') + gpu_data = numpy3d_to_array(data) + gpu_rad_tex.set_array(gpu_data) + + n_scan_angles = np.uint32(n_scan_angles) + n_points = np.uint32(len(Z)) + max_iters = np.uint32(max_iters) + norm_factor = np.float32(1. / n_scan_angles) + + block = (16, 16, 1) + n_blocks = np.ceil(float(n_points)/(block[0] * block[1])) + g_cols = 2 + g_rows = np.int(np.ceil(n_blocks / g_cols)) + grid = (g_rows, g_cols, 1) + + start = cuda.Event() + end = cuda.Event() + + start.record() # start timing + diameter3d(width, height, depth, + n_points, + norm_factor, + max_iters, + n_scan_angles, + gpu_X, + gpu_Y, + gpu_Z, + gpu_scan_angl_arr, + gpu_azth_data_1d, + gpu_lat_data_1d, + gpu_radius_arr, + block=block, + grid=grid) + end.record() # end timing + end.synchronize() + + dm_time = start.time_till(end)*1e-3 + print("Diameter estimation time: %fs" % (dm_time)) + + radius_arr = gpu_radius_arr.get() + + if do_reshape: + radius_arr = np.reshape(radius_arr, data.shape, order='C') + + out = {'diameter': radius_arr * 2., 'time': dm_time} + return out + + +def estimate_diameter_single_run(name, output_dir, data, skel, lat_data, azth_data, + n_scan_angles=32, make_output=True): + """Computes 3D diameter using GPU and stores result in a npy file. + + Parameters + ---------- + name : str + Indicates the name of the output npy file. + + output_dir : str + Indicates the path to the output folder where the data will be stored. + + data : 3D array + Indicates the 3D binary data. + + skel : 3D array + Indicates the skeleton of the binary data. + + lat_data : 3D array + Indicates the 3D array containing latitude / elevation angle at every point of + the skeleton in radians. + + azth_data : 3D array + Indicates the 3D array containing azimuth angle at every point of the skeleton + in radians. + + n_scan_angles : int + Indicates the number of scanning angles on a range [0, 360] degrees. + + make_output : boolean + Specifies if the estimated data should be stored. + + Returns + ------- + output_props : dict + The dictionary of properties specifying the sample name, the algorithm name, + the number of processes, and the execution time. + """ + if not cuda_available: + print('The pycuda package is not found. The diameter estimation cannot be done.') + return None + + output = dict() + output_props = dict() + + output['diameter'] = np.zeros_like(data, dtype=np.float32) + + output_props['time'] = 0. + output_props['sample_name'] = name + output_props['type'] = 'diameter_gpu' + output_props['n_processes'] = 1 + + result = estimate_diameter_gpu(skel, data, lat_data, azth_data, n_scan_angles, do_reshape=False) + + output_props['time'] += result['time'] + output['diameter'][skel.nonzero()] = result['diameter'] + + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + print('Total diameter execution time: {}'%format(output_props['time'])) + + if make_output: + opath = os.path.join(output_dir, '{}.npy').format(name) + output_props['output_path'] = opath + output['indices'] = skel > 0 + output['props'] = output_props + np.save(opath, output) + + return output_props + + +def estimate_diameter_batches(name, output_dir, data, skel, lat_data, azth_data, border_gap, + n_scan_angles=32, out_arr_names=['diameter'], make_output=True, + slices_per_batch=100): + """Computes 3D diameter using GPU in batches and stores result in a npy file. + + Parameters + ---------- + name : str + Indicates the name of the output npy file. + + output_dir : str + Indicates the path to the output folder where the data will be stored. + + data : 3D array + Indicates the 3D binary data. + + skel : 3D array + Indicates the skeleton of the binary data. + + lat_data : 3D array + Indicates the 3D array containing latitude / elevation angle at every point of + the skeleton in radians. + + azth_data : 3D array + Indicates the 3D array containing azimuth angle at every point of the skeleton + in radians. + + border_gap : integer + Indicates the number of overlapping slices along z-axis, usually it should be more or + equal to a half of size of the 3D local window. + + n_scan_angles : int + Indicates the number of scanning angles on a range [0, 360] degrees. + + out_arr_names : array of str + Indicates the array of keys of the output dictionary. + + make_output : boolean + Specifies if the estimated data should be stored. + + slices_per_batch : integer + The number of slices along z-axis in a batch. + + Returns + ------- + output_props : dict + The dictionary of properties specifying the sample name, the algorithm name, + the number of processes, and the execution time. + """ + if not cuda_available: + print('The pycuda package is not found. The diameter estimation cannot be done.') + return None + + output = dict() + output_props = dict() + + for arr_name in out_arr_names: + output[arr_name] = np.zeros_like(data, dtype=np.float32) + + output_props['time'] = 0. + output_props['sample_name'] = name + output_props['type'] = 'diameter_gpu' + output_props['n_processes'] = 1 + + depth, height, width = data.shape + batches_idxs = np.array_split(np.arange(depth), np.ceil(depth / float(slices_per_batch))) + print(batches_idxs) + + for batch_idxs in batches_idxs: + batch_len = len(batch_idxs) + + if batch_idxs[0] == 0: + arr1, arr2 = None, np.arange(batch_idxs[-1]+1, batch_idxs[-1]+border_gap+1) + gaped_batch_idxs = np.concatenate((batch_idxs, arr2)) + elif batch_idxs[-1] == (depth - 1): + arr1, arr2 = np.arange(batch_idxs[0]-border_gap, batch_idxs[0]), None + gaped_batch_idxs = np.concatenate((arr1, batch_idxs)) + else: + arr1, arr2 = np.arange(batch_idxs[0]-border_gap, batch_idxs[0]), \ + np.arange(batch_idxs[-1], batch_idxs[-1]+border_gap+1) + gaped_batch_idxs = np.concatenate((arr1, batch_idxs, arr2)) + + batched_data = data[gaped_batch_idxs] + batched_skel = skel[gaped_batch_idxs] + batched_lat_data = lat_data[gaped_batch_idxs] + batched_azth_data = azth_data[gaped_batch_idxs] + + gaped_out_dict = estimate_diameter_gpu(batched_skel, + batched_data, + batched_lat_data, + batched_azth_data, + n_scan_angles, + do_reshape=False) + + output_props['time'] += gaped_out_dict['time'] + + for data_name in out_arr_names: + gaped_arr_values = gaped_out_dict[data_name] + gaped_arr = np.zeros_like(batched_data, dtype=np.float32) + gaped_arr[batched_skel.nonzero()] = gaped_arr_values + + if batch_idxs[0] == 0: + output[data_name][batch_idxs] = gaped_arr[:batch_len] + elif batch_idxs[-1] == (depth - 1): + output[data_name][batch_idxs] = gaped_arr[border_gap:] + else: + output[data_name][batch_idxs] = gaped_arr[border_gap:border_gap+batch_len] + + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + print('Total diameter execution time: {}'%(output_props['time'])) + + if make_output: + opath = os.path.join(output_dir, '{}.npy').format(name) + output_props['output_path'] = opath + output['indices'] = output['diameter'] > 0 + output['props'] = output_props + np.save(opath, output) + + return output_props + + +def _calc_sphericity(volume, perimeter): + """Computes sphericity from volume and perimeter of an object. + """ + r = ((3.0 * volume) / (4.0 * np.pi)) ** (1.0/3.0) + return (4.0 * np.pi * (r*r)) / perimeter + + +def object_counter(stack_binary_data): + """Label and counts particles in a binary data. + + Parameters + ---------- + stack_binary_data : 3D array + Indicates the 3D binary data. + + Returns + ------- + (objects_stats, labeled_stack) : tuple + The tuple of a DataFrame object of counted partilces and the labeled 3D data. + """ + measurements_vals = _MEASUREMENTS.values() + + print('Object counting - Labeling...') + labeled_stack, num_labels = ndi.measurements.label(stack_binary_data) + objects_stats = pd.DataFrame(columns=measurements_vals) + + print('Object counting - Stats gathering...') + for slice_idx in np.arange(labeled_stack.shape[0]): + for region in measure.regionprops(labeled_stack[slice_idx]): + objects_stats = objects_stats.append({mname: region[mval] + for mname, mval in _MEASUREMENTS.items()}, + ignore_index=True) + + print('Object counting - Stats grouping...') + objects_stats = objects_stats.groupby('Label', as_index=False).sum() + objects_stats['Sphericity'] = _calc_sphericity(objects_stats['Area'], + objects_stats['Perimeter']) + + return (objects_stats, labeled_stack) + + +def calc_porosity(data): + """Computes porosity. + + Parameters + ---------- + data : 3D array + Indicates the labeled 3D data. + + Returns + ------- + out : dict + The dictionary of materials and corresponding porosity values. + """ + total_volume = data.size + mats = np.unique(data) + out = {'Material {}'.format(m): \ + 1. - (float(data[data == m].size) / total_volume) \ + for m in mats[mats > 0]} + return out diff --git a/quanfima/quanfima/simulation.py b/quanfima/quanfima/simulation.py new file mode 100644 index 0000000..1386896 --- /dev/null +++ b/quanfima/quanfima/simulation.py @@ -0,0 +1,676 @@ +from __future__ import print_function +import time +import os +import itertools +import numpy as np +from multiprocessing import Pool +from scipy import ndimage as ndi +from sklearn import metrics +from skimage import filters, morphology, data as skidata, exposure, draw + + +def random_in(rng, number=1): + """Returns a random value within a given range. + """ + start, end = rng + values = np.random.random_sample(number) * (end - start) + start + return values[0] if number == 1 else values + + +def mkfiber(dims_size, length, radius, azth, lat, offset_xyz): + """Computes fiber coordinates and its length. + + Computes a fiber of speficied `length`, `radius`, oriented under azimuth `azth` and + latitude / elevation `lat` angles shifted to `offset_xyz` from the center of a volume of + size `dims_size`. + + Parameters + ---------- + dims_size : tuple + Indicates the size of the volume. + + length : integer + Indicates the length of the simulated fiber. + + radius : integer + Indicates the radius of the simulated fiber. + + azth : float + Indicates the azimuth component of the orientation angles of the fiber in radians. + + lat : float + Indicates the latitude / elevation component of the orientation angles of the fiber + in radians. + + offset_xyz : tuple + Indicates the offset from the center of the volume where the fiber will be generated. + + Returns + ------- + fiber_pts, fiber_len : tuple of array and number + The array of fiber coordinates and the length. + """ + dims_size = np.array(dims_size) + + half_pi = np.pi / 2. + mx = np.array([[1., 0., 0], + [0., np.cos(lat), np.sin(lat)], + [0., -np.sin(lat), np.cos(lat)]]) + + azth += half_pi + mz = np.array([[np.cos(azth), -np.sin(azth), 0], + [np.sin(azth), np.cos(azth), 0], + [0., 0., 1.]]) + + # Directional vector + dl = 1 + dir_vec = np.array([0, 0, 1]) + dir_vec = np.dot(mx, dir_vec) + dir_vec = np.dot(mz, dir_vec) + dx, dy, dz = dir_vec[0], dir_vec[1], dir_vec[2] + + # Compute length + n_steps = np.round(length / dl) + half_steps = int(np.ceil(n_steps / 2.)) + steps = range(half_steps - int(n_steps), half_steps) + + # Draw circle perpedicular to the directional vector + X, Y = draw.circle(0, 0, radius) + Z = np.repeat(0, len(Y)) + circle_pts = np.array([X, Y, Z]) + circle_pts = np.dot(mx, circle_pts) + circle_pts = np.dot(mz, circle_pts) + + # Propogate the circle profile along the directional vector + slice_pts = np.array(zip(*circle_pts)) + dxyz = np.array([dx, dy, dz]) + step_shifts = np.array([step * dxyz for step in steps]) # [(dx,dy,dz), ...] for each step + center_shift = dims_size * 0.5 + offset_xyz # (x, y ,z) + slices_pts = np.round(np.array([slice_pts + step_shift + center_shift + for step_shift in step_shifts])) + + # Filter all the points which are outside the boundary + pt_filter = lambda pt: np.all(np.greater_equal(pt, (0, 0, 0))) and \ + np.all(np.less(np.array(pt), dims_size)) + n_slices = 0 + fiber_pts = None + for pts in slices_pts: + slice_pts_mask = [pt_filter(pt) for pt in pts] + slice_pts = pts[slice_pts_mask].astype(np.int32) + if len(slice_pts) > 0: + n_slices += 1 + fiber_pts = slice_pts if fiber_pts is None else \ + np.concatenate((fiber_pts, slice_pts)) + # number of slices, e.g. steps. + fiber_len = np.round(n_slices * dl).astype(np.int32) + fiber_pts = fiber_pts.astype(np.int32) + return fiber_pts, fiber_len + + +def simulate_fibers(volume_shape, n_fibers=1, radius_lim=(4, 10), length_lim=(0.2, 0.8), + lat_lim=(0, np.pi), azth_lim=(0, np.pi), gap_lim=(3, 10), + max_fails=10, max_len_loss=0.5, intersect=False): + """Simulates fibers in a volume. + + Simulates `n_fibers` of the radii and lengths in ranges `radius_lim` and `length_lim`, + oriented in a range of azimuth `azth_lim` and latitude \ elevation 'lat_lim' angles, + separated with a gap in a range of `gap_lim`. The simulation process stops if the number + of attempts to generate a fiber exceeds `max_fails`. + + Parameters + ---------- + volume_shape : tuple + Indicates the size of the volume. + + n_fibers : integer + Indicates the number of fibers to be generated. + + radius_lim : tuple + Indicates the range of radii for fibers to be generated. + + length_lim : tuple + Indicates the range of lengths for fibers to be generated. + + lat_lim : tuple + Indicates the range of latitude / elevation component of the orientation angles of + the fibers to be generated. + + azth_lim : tuple + Indicates the range of azimuth component of the orientation angles of the fibers to + be generated. + + gap_lim : tuple + Indicates the range of gaps separating the fibers from each other. + + max_fails : integer + Indicates the maximum number of failures during the simulation process. + + max_len_loss : float + Indicates the maximum fraction of the generated fiber placed out of volume, exceeding + which the fiber is counted as failed. + + intersect : boolean + Specifies if generated fibers can intersect. + + Returns + ------- + (volume, lat_ref, azth_ref, diameter, n_generated, elapsed_time) : tuple of arrays and numbers + The binary volume of generated fibers, the volumes of latitude / elevation and + azimuth angles at every fiber point, the volume of diameters at every fibers point, + the number of generated fibers and the simulation time. + """ + ts = time.time() + + volume = np.zeros(volume_shape, dtype=np.int32) + lat_ref = np.zeros_like(volume, dtype=np.float32) + azth_ref = np.zeros_like(volume, dtype=np.float32) + diameter = np.zeros_like(volume, dtype=np.float32) + + dims = np.array(volume.shape)[::-1] + offset_lim = zip(-dims / 2, dims / 2) + + n_generated = 0 + n_fails = 0 + while n_generated < n_fibers and n_fails < max_fails: + length = min(volume_shape) + length = np.floor(length * random_in(length_lim, number=1)).astype(np.int32) + + offset = [random_in(olim, number=1) for olim in offset_lim] + offset = np.round(offset).astype(np.int32) + + azth = random_in(azth_lim, number=1) + lat = random_in(lat_lim, number=1) + radius = random_in(radius_lim, number=1) + + gap = random_in(gap_lim, number=1) + gap = np.round(gap).astype(np.int32) + + fiber_pts, fiber_len = mkfiber(dims, length, radius, azth, lat, offset) + gap_fiber_pts, gap_fiber_len = mkfiber(dims, length, radius + gap, azth, lat, offset) + + # Length loss + if (1. - float(gap_fiber_len) / length) > max_len_loss: + n_fails = n_fails + 1 + continue + + # Intersection + if gap_fiber_pts.size: + X_gap, Y_gap, Z_gap = gap_fiber_pts[:, 0], gap_fiber_pts[:, 1], gap_fiber_pts[:, 2] + X, Y, Z = fiber_pts[:, 0], fiber_pts[:, 1], fiber_pts[:, 2] + + if not intersect: + if np.any(volume[Z_gap, Y_gap, X_gap]): + n_fails = n_fails + 1 + + if n_fails == max_fails: + print("The number of fails exceeded. Generated {} fibers".\ + format(n_generated)) + + continue + + # Fill the volume + volume[Z, Y, X] = 1 + lat_ref[Z, Y, X] = lat + azth_ref[Z, Y, X] = azth + diameter[Z, Y, X] = radius * 2 + n_generated = n_generated + 1 + n_fails = 0 + + te = time.time() + elapsed_time = te - ts + + return (volume, lat_ref, azth_ref, diameter, n_generated, elapsed_time) + + +def generate_datasets(volume_size=(512, 512, 512), n_fibers=50, radius_lim=(4, 10), + length_lim=(0.2, 0.8), gap_lim=(3, 10), max_fails=100, + median_rad=3, intersect=False, output_dir=None, params=None): + """Simulates speficied configurations of fibers and stores in a npy file. + + Simulates a number of fiber configurations speficied in `params` with `n_fibers` of the + radii and lengths in ranges `radius_lim` and `length_lim`, separated with gaps in a range + of `gap_lim`. The simulation process stops if the number of attempts to generate a fiber + exceeds `max_fails`. + + Parameters + ---------- + volume_size : tuple + Indicates the size of the volume. + + n_fibers : integer + Indicates the number of fibers to be generated. + + radius_lim : tuple + Indicates the range of radii for fibers to be generated. + + length_lim : tuple + Indicates the range of lengths for fibers to be generated. + + gap_lim : tuple + Indicates the range of gaps separating the fibers from each other. + + max_fails : integer + Indicates the maximum number of failures during the simulation process. + + median_rad : integer + Indicates the radius of median filter to fill holes occured due to rounding of + coordinates of the generated fibers. + + intersect : boolean + Specifies if generated fibers can intersect. + + output_dir : str + Indicates the path to the output folder where the data will be stored. + + params : dict + Indicates the configurations of orientation of fibers to be generated. + + Returns + ------- + out : dict + The dictionary of generated datasets of specified configurations. + """ + if params is None: + params = {'aligned': {'lat_rng': (15, 15), 'azth_rng': (27, 27)}, + 'medium': {'lat_rng': (0, 45), 'azth_rng': (-45, 45)}, + 'disordered': {'lat_rng': (0, 90), 'azth_rng': (-89, 90)}} + + out = {} + for name, config in params.items(): + data, lat_data, azth_data, diameter_data, n_gen_fibers, elapsed_time = \ + simulate_fibers(volume_size, + lat_lim=tuple([np.deg2rad(v) for v in config['lat_rng']]), + azth_lim=tuple([np.deg2rad(v) for v in config['azth_rng']]), + radius_lim=radius_lim, + n_fibers=n_fibers, + max_fails=max_fails, + gap_lim=gap_lim, + length_lim=length_lim, + intersect=intersect) + + data_8bit = data.astype(np.uint8) + data_8bit = ndi.binary_fill_holes(data_8bit) + data_8bit = ndi.median_filter(data_8bit, footprint=morphology.ball(median_rad)) + lat_data = ndi.median_filter(lat_data, footprint=morphology.ball(median_rad)) + azth_data = ndi.median_filter(azth_data, footprint=morphology.ball(median_rad)) + diameter_data = ndi.median_filter(diameter_data, footprint=morphology.ball(median_rad)) + + out[name] = {'data': data_8bit, + 'lat': lat_data, + 'azth': azth_data, + 'diameter': diameter_data, + 'skeleton': morphology.skeletonize_3d(data_8bit).astype(np.float32), + 'props': {'n_gen_fibers': n_gen_fibers, + 'time': elapsed_time, + 'intersect': intersect}} + + if output_dir is not None and not os.path.exists(output_dir): + os.makedirs(output_dir) + + np.save(os.path.join(output_dir, 'dataset_fibers_n{}_r{}_{}_g{}_{}_mr{}_i{}.npy'. + format(n_fibers, + radius_lim[0], radius_lim[-1], + gap_lim[0], gap_lim[-1], + median_rad, int(intersect))), out) + + return out + + +def simulate_particles(volume_shape, n_particles=1, radius_lim=(3, 30), max_fails=10): + """Simulates particles in a volume. + + Simulates `n_particles` of the radii in a range `radius_lim`. The simulation process + stops if the number of attempts to generate a particle exceeds `max_fails`. + + Parameters + ---------- + volume_shape : tuple + Indicates the size of the volume. + + n_particles : integer + Indicates the number of particles to be generated. + + radius_lim : tuple + Indicates the range of radii for particles to be generated. + + max_fails : integer + Indicates the maximum number of failures during the simulation process. + + Returns + ------- + (volume, diameter, n_generated, elapsed_time) : tuple of arrays and numbers + The binary volume of generated particles, the volume of diameters at every point of + particles, the number of generated particles and the simulation time. + """ + ts = time.time() + + volume = np.zeros(volume_shape, dtype=np.uint8) + diameter = np.zeros_like(volume, dtype=np.int32) + + dims = np.array(volume.shape) + offset_lim = zip(itertools.repeat(0), dims) + + n_generated = 0 + n_fails = 0 + while n_generated < n_particles and n_fails < max_fails: + if (n_generated % 100 == 0) or (n_generated == n_particles): + print('n_generated = {}/{}, n_fails = {}/{}'.format(n_generated, n_particles, + n_fails, max_fails)) + + offset = [random_in(olim, number=1) for olim in offset_lim] + offset = np.round(offset).astype(np.int32) + + radius = np.round(random_in(radius_lim, number=1)) + + gen_ball = morphology.ball(radius, dtype=np.int32) + Z, Y, X = gen_ball.nonzero() + + Z += offset[0] + Y += offset[1] + X += offset[2] + + if np.max(X) >= dims[0] or np.min(X) < 0 or \ + np.max(Y) >= dims[1] or np.min(Y) < 0 or \ + np.max(Z) >= dims[2] or np.min(Z) < 0: + + n_fails = n_fails + 1 + continue + + if np.any(volume[Z, Y, X]): + n_fails = n_fails + 1 + + if n_fails == max_fails: + print("The number of fails exceeded. Generated {} particles".\ + format(n_generated)) + + continue + + # Fill the volume + volume[Z, Y, X] = 1 + diameter[Z, Y, X] = radius * 2 + n_generated = n_generated + 1 + n_fails = 0 + + te = time.time() + elapsed_time = te - ts + + return (volume, diameter, n_generated, elapsed_time) + + +def generate_particle_dataset(volume_size=(512, 512, 512), n_particles=500, + radius_lim=(4, 10), max_fails=100, output_dir=None): + """Simulates a speficied number of particles and stores complete dataset in a npy file. + + Parameters + ---------- + volume_size : tuple + Indicates the size of the volume. + + n_particles : integer + Indicates the number of particles to be generated. + + radius_lim : tuple + Indicates the range of radii for particles to be generated. + + max_fails : integer + Indicates the maximum number of failures during the simulation process. + + output_dir : str + Indicates the path to the output folder where the data will be stored. + + Returns + ------- + out : dict + The dictionary of generated dataset. + """ + out = {} + + data, diameter_data, n_gen_particle, elapsed_time = \ + simulate_particles(volume_size, + n_particles=n_particles, + radius_lim=radius_lim, + max_fails=max_fails) + + out['normal'] = {'data': data, + 'diameter': diameter_data, + 'props': {'n_gen_particles': n_gen_particle, + 'time': elapsed_time}} + + if output_dir is not None and not os.path.exists(output_dir): + os.makedirs(output_dir) + + np.save(os.path.join(output_dir, + 'dataset_particles_n{}_r{}_{}.npy'.format(n_particles, + radius_lim[0], + radius_lim[-1])), out) + return out + + +def generate_blobs(volume_size, blob_size_fraction=0.1, transparency_ratio=0.5, sigma=90.): + """Generates random blobs smoothed with Gaussian filter in a volume. + + Generates several blobs of random size in a volume using function from scikit-image, + which are subsequently smoothed with a Gaussian filter of a large sigma to imitate + 3D uneven illumination of the volume. + + Parameters + ---------- + volume_size : tuple + Indicates the size of the volume. + + blob_size_fraction : float + Indicates the fraction of volume occupied by blobs. + + transparency_ratio : float + Indicates the transparency of blobs in a range [0, 1]. + + sigma : float + Indicates the sigma of Gaussian filter. + + Returns + ------- + blobs_smeared : ndarray + The volume with smoothed blobs of a specified transparency. + """ + blobs = skidata.binary_blobs(length=max(volume_size), + blob_size_fraction=blob_size_fraction, + n_dim=len(volume_size), + seed=1) + blobs = blobs.astype(np.float32) + + blobs_smeared = ndi.gaussian_filter(blobs, sigma) * transparency_ratio + return blobs_smeared + + +def generate_noised_data(datasets_path, noise_levels=[0.0, 0.15, 0.3], + smooth_levels=[0.0, 1.0, 2.0], blobs=None, use_median=True, + median_rad=3, output_dir=None, n_processes=9): + """Contaminates datasets with a speficied additive Gaussian noise and smoothing level. + + Contaminates the datasets (generated with `generate_datasets` function) with the specified + level of additive Gaussian noise and smoothing, uneven illumination can be added if + `blobs` is provided. The contaminating process can be performed in a parallel `n_processes` + processes. + + Parameters + ---------- + datasets_path : str + Indicates the path to dataset. + + noise_levels : array + Indicates the array of standard deviations of noise. + + smooth_levels : array + Indicates the array of sigma values of Gaussian filter. + + blobs : ndarray + Indicates the volume of uneven illumination generated by `generate_blobs`. + + use_median : boolean + Specifies if the median filter should be applied after addition of noise. + + median_rad : integer + Indicates the size of median filter. + + output_dir : str + Indicates the path to the output folder where the data will be stored. + + n_processes : integer + Indicates the number of parallel processes. + + Returns + ------- + results : array of dicts + The array of dictionaries containing paths to contaminated datasets, and other + properties. + """ + datasets_names = np.load(datasets_path).item().keys() + n_datasets = len(datasets_names) + + dataset_filename = os.path.basename(datasets_path) + dataset_filename = os.path.splitext(dataset_filename)[0] + output_dir = os.path.join(output_dir, dataset_filename) + + data_items = [(dname, dpath, blb, odir) + for dname, dpath, blb, odir in zip(datasets_names, + [datasets_path]*n_datasets, + [blobs]*n_datasets, + [output_dir]*n_datasets)] + + params = [data_items, noise_levels, smooth_levels] + params = [p for p in itertools.product(*params)] + + proc_pool = Pool(processes=n_processes) + results = proc_pool.map(unpack_additive_noise, params) + proc_pool.close() + proc_pool.join() + proc_pool.terminate() + + np.save(os.path.join(output_dir, 'params.npy'), results) + + return results + + +def unpack_additive_noise(args): + """Unpack arguments and return result of `additive_noise` function. + """ + return additive_noise(*args) + + +def additive_noise(params, noise_lvl, smooth_lvl, use_median=True, median_rad=3): + """Contaminates datasets with a speficied additive Gaussian noise and smoothing level. + + Contaminates the datasets (generated with `generate_datasets` function) with the specified + level of additive Gaussian noise and smoothing, uneven illumination can be added by + extracting `blobs` from `params` tuple with some other arguments. + + Parameters + ---------- + params : tuple + Contains `name`, `dataset_path`, `blobs`, `output_dir` arguments passed from + `generate_noised_data`. + + noise_level : float + Indicates the standard deviations of noise. + + smooth_level : float + Indicates the sigma value of Gaussian filter. + + use_median : boolean + Specifies if the median filter should be applied after addition of noise. + + median_rad : integer + Indicates the size of median filter. + + Returns + ------- + datasets_props : dict + The dictionary containing the path to the reference dataset, the path to the + contaminated dataset, the generated name, the SNR level, the precision, the recall + and f1-score, and the level of noise and smoothing. + """ + name, dataset_path, blobs, output_dir = params + + datasets = np.load(dataset_path).item() + data = datasets[name]['data'] + data_skel = datasets[name]['skeleton'] + + def median_fltr(data, footprint): + out = np.zeros_like(data, dtype=np.uint8) + + for i in xrange(data.shape[0]): + out[i] = filters.rank.median(data[i], selem=footprint) + + return out + + def threshold_dataset(data): + data_seg = np.zeros_like(data, dtype=np.uint8) + data_8bit = exposure.rescale_intensity(data, in_range='image', + out_range=np.uint8).astype(np.uint8) + + for i in xrange(data_seg.shape[0]): + dslice = data_8bit[i] + th_val = filters.threshold_otsu(dslice) + data_seg[i] = (dslice > th_val).astype(np.uint8) + + return data_seg + + print('{}: Noise: {} | Smooth: {}'.format(name, noise_lvl, smooth_lvl)) + + data_ref = data.astype(np.float32) + data_ref_skel = exposure.rescale_intensity(data_skel, + out_range=(0, 1)).astype(np.uint8) + + data_noised = data_ref + noise_lvl * np.random.randn(*data_ref.shape) + + if (blobs is not None) and (noise_lvl != 0.) and (smooth_lvl != 0.): + data_noised += blobs + + if smooth_lvl != 0: + data_noised = ndi.gaussian_filter(data_noised, smooth_lvl) + + snr = np.mean(data_noised[data_ref != 0]) / np.std(data_noised[data_ref == 0]) + data_noised = exposure.rescale_intensity(data_noised, out_range=np.uint8).astype(np.uint8) + + if use_median and (noise_lvl != 0.) and (smooth_lvl != 0.): + data_noised = median_fltr(data_noised, morphology.disk(median_rad)) + + data_noised_seg = threshold_dataset(data_noised) + data_noised_skel = morphology.skeletonize_3d(data_noised_seg) + + precision, recall, fbeta_score, support = \ + metrics.precision_recall_fscore_support(data_ref_skel.flatten(), + data_noised_skel.flatten(), + beta=1.0, + pos_label=1, + average='binary') + + print('Precision: {}, Recall: {}, F1-score: {}'.format(precision, recall, fbeta_score)) + + data_out = {'data': data_ref, + 'data_noised': data_noised, + 'skeleton': data_ref_skel, + 'skeleton_noised': data_noised_skel, + 'seg_noised': data_noised_seg} + + dataset_outpath = os.path.join(output_dir, + 'dataset_noised_fibers_{}_nl{}_sl{}.npy'. + format(name, noise_lvl, smooth_lvl)) + + datasets_props = {'ref_dataset_path': dataset_path, + 'dataset_path': dataset_outpath, + 'name': name, + 'snr': snr, + 'precision': precision, + 'recall': recall, + 'f1_score': fbeta_score, + 'adgauss_std': noise_lvl, + 'smooth_sigma': smooth_lvl} + + if output_dir is not None: + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + np.save(dataset_outpath, data_out) + + return datasets_props diff --git a/quanfima/quanfima/test.txt b/quanfima/quanfima/test.txt new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/quanfima/quanfima/test.txt @@ -0,0 +1 @@ + diff --git a/quanfima/quanfima/test/test.txt b/quanfima/quanfima/test/test.txt new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/quanfima/quanfima/test/test.txt @@ -0,0 +1 @@ + diff --git a/quanfima/quanfima/utils.py b/quanfima/quanfima/utils.py new file mode 100644 index 0000000..1d21215 --- /dev/null +++ b/quanfima/quanfima/utils.py @@ -0,0 +1,150 @@ +from __future__ import print_function +import os +import operator +import numpy as np +from scipy import ndimage as ndi +from skimage import morphology +from matplotlib import colors +from statsmodels.stats.multicomp import MultiComparison +from statsmodels.stats.libqsturng import psturng + + +def prepare_data(data, dilate_iterations=1, sigma=0.5): + """Returns the given binary data, its skeleton and the thickened skeleton. + + The skeleton of a given 2D or 3D array is computed, then it is thickened + using morphological dilation with `dilate_iterations` and smoothed with + help of Gaussian filter of specified `sigma`. + + Parameters + ---------- + data : ndarray + 2D or 3D binary array which will be processed. + + dilate_iterations : integer + Indicates the number of iterations for thickenning the skeleton. + + sigma : float + Indicates the sigma of Gaussian filter used in smoothing of skeleton. + + Returns + ------- + arrays : tuple of 2D or 3D arrays + The original array, its skeleton and the thickened skeleton. + """ + data_8bit = data.astype(np.uint8) + data_8bit = ndi.binary_fill_holes(data_8bit).astype(np.uint8) + + if data.ndim == 3: + skeleton = morphology.skeletonize_3d(data_8bit) + elif data.ndim == 2: + skeleton = morphology.skeletonize(data_8bit) + else: + raise ValueError('Incorrect number of data dimensions, it supports from 2 to 3 dimensions.') + + skeleton_thick = ndi.binary_dilation(skeleton, iterations=dilate_iterations).astype(np.float32) + skeleton_thick = ndi.filters.gaussian_filter(skeleton_thick, sigma) + + return (data, skeleton, skeleton_thick) + + +def geo2rgb(lat, azth, azth_max=np.pi, lat_max=np.pi): + """Translates geo-coordinates to color in RGB color space. + + Parameters + ---------- + lat : float + Indicates latitude or elevation component of a given geo-coordinates. + + azth : float + Indicates azimuth component of a given geo-coordinates. + + azth_max : float + Indicates the normalization value for the azimuth component. + + lat_max : float + Indicates the normalization value for the latitude or elevation component. + + Returns + ------- + array : tuple of values + The tuple of RGB values [R, G, B]. + """ + return colors.hsv_to_rgb([azth/azth_max, lat/lat_max, 1.0]) + + +def calculate_tukey_posthoc(df, column, type_column='type', verbose=True, write=False, name=None, output_dir=None): + """Computes p-values using ANOVA with post-hoc Tukey HSD for a given DataFrame. + + Estimates p-values for a given DataFrame assuming that the sample + type is named as `type_column`. + + Parameters + ---------- + df : pandas DataFrame + Contains the table of values ans corresponding types or classes. + + column : str + Indicates the column of values. + + type_column : str + Indicates the column of sample kind. + + verbose : boolean + Specifies if the output should be printed into a terminal. + + write : boolean + Specifies if the output should be written into a text file. + + name : str + Indicates the name of the output file. + + output_dir : str + Indicates the output dir where the file will be written. + + Returns + ------- + dict : sample typles and p-values + The dict of sample types and cooresponding p-values. + """ + mc = MultiComparison(df[column], df[type_column]) + tt = mc.tukeyhsd() + st_range = np.abs(tt.meandiffs) / tt.std_pairs + + fout = None + if write and output_dir: + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + fout = open(os.path.join(output_dir, name + '.txt'), 'w') + print(os.path.join(output_dir, name + '.txt')) + + if write: + print('Tukey post-hoc ({0})'.format(column), end="", file=fout) + print(tt, end="", file=fout) + print(mc.groupsunique, end="", file=fout) + + if verbose: + print('Tukey post-hoc ({0})'%(column)) + print(tt) + print(mc.groupsunique) + + pvals = psturng(st_range, len(tt.groupsunique), tt.df_total) + + out = {} + groups = mc.groupsunique + g1idxs, g2idxs = mc.pairindices + + for g1i, g2i, p in zip(g1idxs, g2idxs, pvals): + gname = '{}-{}'.format(groups[g1i], groups[g2i]) + out[gname] = p + + min_item = min(out.iteritems(), key=operator.itemgetter(1)) + + for grp, p in out.items(): + if fout and write: + print >> fout, '{}: {}'.format(grp, p) + if verbose: + print(grp, ': ', p) + + return out, min_item diff --git a/quanfima/quanfima/utils2.py b/quanfima/quanfima/utils2.py new file mode 100644 index 0000000..a5bb792 --- /dev/null +++ b/quanfima/quanfima/utils2.py @@ -0,0 +1,156 @@ +from __future__ import print_function +import os +import operator +import numpy as np +from scipy import ndimage as ndi +from skimage import morphology +from matplotlib import colors +from statsmodels.stats.multicomp import MultiComparison +from statsmodels.stats.libqsturng import psturng +import matplotlib.pyplot as plt +from skimage import data +from skimage.filters import threshold_otsu +import cv2 + + +def prepare_data(data, dilate_iterations=1, sigma=0.5): + """Returns the given binary data, its skeleton and the thickened skeleton. + + The skeleton of a given 2D or 3D array is computed, then it is thickened + using morphological dilation with `dilate_iterations` and smoothed with + help of Gaussian filter of specified `sigma`. + + Parameters + ---------- + data : ndarray + 2D or 3D binary array which will be processed. + + dilate_iterations : integer + Indicates the number of iterations for thickenning the skeleton. + + sigma : float + Indicates the sigma of Gaussian filter used in smoothing of skeleton. + + Returns + ------- + arrays : tuple of 2D or 3D arrays + The original array, its skeleton and the thickened skeleton. + """ + data_8bit = data.astype(np.uint8) + data_8bit = ndi.binary_fill_holes(data_8bit).astype(np.uint8) + + if data.ndim == 3: + skeleton = data#morphology.skeletonize_3d(data_8bit) + plt.imshow(skeleton) + elif data.ndim == 2: + skeleton = data# morphology.skeletonize(data_8bit) + else: + raise ValueError('Incorrect number of data dimensions, it supports from 2 to 3 dimensions.') + + skeleton_thick = ndi.binary_dilation(skeleton, iterations=dilate_iterations).astype(np.float32) + skeleton_thick = ndi.filters.gaussian_filter(skeleton_thick, sigma) + plt.imshow(skeleton_thick) + + return (data, skeleton, skeleton_thick) + + +def geo2rgb(lat, azth, azth_max=np.pi, lat_max=np.pi): + """Translates geo-coordinates to color in RGB color space. + + Parameters + ---------- + lat : float + Indicates latitude or elevation component of a given geo-coordinates. + + azth : float + Indicates azimuth component of a given geo-coordinates. + + azth_max : float + Indicates the normalization value for the azimuth component. + + lat_max : float + Indicates the normalization value for the latitude or elevation component. + + Returns + ------- + array : tuple of values + The tuple of RGB values [R, G, B]. + """ + return colors.hsv_to_rgb([azth/azth_max, lat/lat_max, 1.0]) + + +def calculate_tukey_posthoc(df, column, type_column='type', verbose=True, write=False, name=None, output_dir=None): + """Computes p-values using ANOVA with post-hoc Tukey HSD for a given DataFrame. + + Estimates p-values for a given DataFrame assuming that the sample + type is named as `type_column`. + + Parameters + ---------- + df : pandas DataFrame + Contains the table of values ans corresponding types or classes. + + column : str + Indicates the column of values. + + type_column : str + Indicates the column of sample kind. + + verbose : boolean + Specifies if the output should be printed into a terminal. + + write : boolean + Specifies if the output should be written into a text file. + + name : str + Indicates the name of the output file. + + output_dir : str + Indicates the output dir where the file will be written. + + Returns + ------- + dict : sample typles and p-values + The dict of sample types and cooresponding p-values. + """ + mc = MultiComparison(df[column], df[type_column]) + tt = mc.tukeyhsd() + st_range = np.abs(tt.meandiffs) / tt.std_pairs + + fout = None + if write and output_dir: + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + fout = open(os.path.join(output_dir, name + '.txt'), 'w') + print(os.path.join(output_dir, name + '.txt')) + + if write: + print('Tukey post-hoc ({0})'.format(column), end="", file=fout) + print(tt, end="", file=fout) + print(mc.groupsunique, end="", file=fout) + + if verbose: + print('Tukey post-hoc ({0})'%(column)) + print(tt) + print(mc.groupsunique) + + pvals = psturng(st_range, len(tt.groupsunique), tt.df_total) + + out = {} + groups = mc.groupsunique + g1idxs, g2idxs = mc.pairindices + + for g1i, g2i, p in zip(g1idxs, g2idxs, pvals): + gname = '{}-{}'.format(groups[g1i], groups[g2i]) + out[gname] = p + + min_item = min(out.iteritems(), key=operator.itemgetter(1)) + + for grp, p in out.items(): + if fout and write: + print >> fout, '{}: {}'.format(grp, p) + if verbose: + print(grp, ': ', p) + + return out, min_item diff --git a/quanfima/quanfima/visualization.py b/quanfima/quanfima/visualization.py new file mode 100644 index 0000000..7d3c156 --- /dev/null +++ b/quanfima/quanfima/visualization.py @@ -0,0 +1,1082 @@ +from __future__ import print_function +import os +import math +import numpy as np +import matplotlib as mpl +from skimage import io, filters +import matplotlib.colors as colors +import matplotlib.pyplot as plt +import matplotlib.ticker as ticker +import mpl_toolkits.axisartist.floating_axes as floating_axes +from matplotlib import cm +from scipy import ndimage as ndi, interpolate +from skimage import morphology +from matplotlib.projections import PolarAxes +from matplotlib.transforms import Affine2D +from mpl_toolkits.axisartist.grid_finder import DictFormatter, FixedLocator +from mpl_toolkits.axes_grid1 import make_axes_locatable +from mpl_toolkits.axisartist.floating_axes import GridHelperCurveLinear +from utils import geo2rgb +from quanfima import visvis_available +if visvis_available: + import visvis as vv + + +def plot_orientation_map(orient_map, fiber_skel, radius_structure_elem=1, + figsize=(12, 12), cmap='hsv', dpi=200, min_label='0', + max_label='180', name=None, output_dir=None): + """Plots the orientation map with the color wheel. + + Plots the orientation map from the provided angles `orient_map` and skeleton + `fiber_skel` of size `figsize` using the colormap `cmap` and writes as a png + file with DPI of `dpi` to the folder specified by `output_dir`. + + Parameters + ---------- + orient_map : ndarray + 2D array of orientation at every point of the skeleton. + + fiber_skel : ndarray + The binary skeleton extracted from the binary data. + + radius_structure_elem : integer + Indicates the size of the structure element of the dilation process to + thicken the skeleton. + + figsize : tuple of integers + Indicates the size of the output figure. + + cmap : str + Indicates the name of a colormap used to map angles to colors. + + dpi : integer + Indicates the DPI of the output image. + + min_label : str + Indicates the label of minimum degree. + + max_label : str + Indicates the label of minimum degree. + + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + """ + disk = morphology.disk(radius_structure_elem) + orient_map = ndi.grey_dilation(orient_map, structure=disk).astype(np.float32) + fiber_skel = ndi.binary_dilation(fiber_skel, structure=disk).astype(np.float32) + + fig = plt.figure(figsize=figsize) + masked_orient_map = np.ma.masked_where(fiber_skel == 0, orient_map) + + quant_steps = 2056 + cmap_obj = cm.get_cmap(cmap, quant_steps) + cmap_obj.set_bad(color='black') + + omap_ax = fig.add_axes([0.0, 0.0, 1.0, 1.0]) + omap_ax.set_axis_off() + omap_ax.imshow(masked_orient_map, cmap=cmap_obj, interpolation=None) + + display_axes = fig.add_axes([0.780, -0.076, 0.2, 0.2], projection='polar') + display_axes._direction = np.pi + + norm = mpl.colors.Normalize(0.0, np.pi) + cb = mpl.colorbar.ColorbarBase(display_axes, cmap=cmap_obj, norm=norm, orientation='horizontal') + + display_axes.text(0.09, 0.56, min_label, color='white', fontsize=20, weight='bold', horizontalalignment='center', + verticalalignment='center', transform=display_axes.transAxes) + + display_axes.text(0.85, 0.56, max_label, color='white', fontsize=20, weight='bold', horizontalalignment='center', + verticalalignment='center', transform=display_axes.transAxes) + + cb.outline.set_visible(False) + display_axes.set_axis_off() + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + fig.savefig(os.path.join(output_dir, '{}_orientation_map.png'.format(name)), + transparent=False, + bbox_inches='tight', + pad_inches=0.1, + dpi=dpi) + + plt.show() + + +def plot_diameter_map(thickness_map, fiber_skel, radius_structure_elem=1, + figsize=(15, 15), cmap='hsv', tmin=None, tmax=None, dpi=200, + labelsize=20, label='Diameter, [pixels]', name=None, + output_dir=None): + """Plots the diameter map with the colorbar. + + Plots the diameter map from the provided diameters `thickness_map` and + skeleton `fiber_skel` of size `figsize` using the colormap `cmap` and the + limits of the colorbar specified by `tmin` and `tmax`, and writes as a png + file with DPI of `dpi` to the folder specified by `output_dir`. + + Parameters + ---------- + thickness_map : ndarray + 2D array of diameter at every point of the skeleton. + + fiber_skel : ndarray + The binary skeleton extracted from the binary data. + + radius_structure_elem : integer + Indicates the size of the structure element of the dilation process to + thicken the skeleton. + + figsize : tuple of integers + Indicates the size of the output figure. + + cmap : str + Indicates the name of a colormap used to map angles to colors. + + tmin : float + Indicates the minimum value of the colorbar. + + tmax : float + Indicates the maximum value of the colorbar. + + dpi : integer + Indicates the DPI of the output image. + + labelsize : integer + Indicates the fontsize of the label of the colorbar. + + label : str + Indicates the label of the colorbar. + + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + """ + disk = morphology.disk(radius_structure_elem) + thickness_map = ndi.grey_dilation(thickness_map, structure=disk).astype(np.float32) + fiber_skel = ndi.binary_dilation(fiber_skel, structure=disk).astype(np.float32) + + masked_thickness_map = np.ma.masked_where(fiber_skel == 0, thickness_map) + + cmap_obj = cm.get_cmap(cmap) + cmap_obj.set_bad(color='black') + + fig = plt.figure(figsize=figsize) + ax = plt.subplot(111) + ax.set_axis_off() + + im = ax.imshow(masked_thickness_map, cmap=cmap_obj, vmin=tmin, vmax=tmax) + divider = make_axes_locatable(ax) + cax = divider.append_axes('right', size="2.5%", pad=0.05) + cbar = plt.colorbar(im, cax=cax) + cbar.ax.tick_params(labelsize=labelsize) + cbar.set_label(label, fontsize=labelsize) + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + fig.savefig(os.path.join(output_dir, '{}_diameter_map.png'.format(name)), + transparent=False, + bbox_inches='tight', + pad_inches=0.1, + dpi=dpi) + plt.show() + + +def gather_polar_errors(datasets_path, estimation_path, + azth_rng=np.arange(-90, 91, step=3), + lat_rng=np.arange(0, 91, step=3)): + """Computes the absolute angular error in ranges between estimated datasets. + + Calculates the absolute angular error in ranges `azth_rng` and `lat_rng` + between estimated orientation datasets localed at paths `datasets_path` and + `estimation_path`. + + Parameters + ---------- + datasets_path : str + Indicates the path to the reference / estimated dataset. + + estimation_path : str + Indicates the path to the estimated / reference dataset. + + azth_rng : array + Indicates the ranges of azimuth angles where the error is accumulated. + + lat_rng : array + Indicates the ranges of latitude or elevation angles where the error + is accumulated. + + Returns + ------- + array : 2D array + The 2D array of accumulated errors within the combinations of the specified + angular ranges. + """ + reference_dataset = np.load(datasets_path).item() + estimated_dataset = np.load(estimation_path).item() + idxs = estimated_dataset['indices'] + + print(reference_dataset.keys()) + print(estimated_dataset.keys()) + + ref_azth, ref_lat = reference_dataset['azth'][idxs], reference_dataset['lat'][idxs] + res_azth, res_lat = estimated_dataset['azth'][idxs], estimated_dataset['lat'][idxs] + + abs_err_azth, abs_err_lat = np.abs(ref_azth - res_azth), np.abs(ref_lat - res_lat) + + out = _angle_err2mean_abs_err(abs_err_azth, abs_err_lat, + ref_azth, ref_lat, + azth_rng=azth_rng, lat_rng=lat_rng) + + return out + + +def _angle_err2mean_abs_err(azth_err, lat_err, azth_ref, lat_ref, + azth_rng=np.arange(-90, 91, step=5), + lat_rng=np.arange(0, 91, step=5)): + """Computes the absolute angular error in ranges between estimated datasets. + + Parameters + ---------- + azth_err : 3D array + Indicates the error of azimuth components. + + lat_err : 3D array + Indicates the error of latitude / elevation components. + + azth_ref : 3D array + Indicates the reference dataset of estimated azimuth component. + + lat_ref : 3D array + Indicates the reference dataset of estimated latitude or elevation component. + + azth_rng : array + Indicates the ranges of azimuth angles where the error is accumulated. + + lat_rng : array + Indicates the ranges of latitude and elevation angles where the error + is accumulated. + + Returns + ------- + array : 2D array + The 2D array of accumulated errors within the combinations of the specified + angular ranges. + """ + out = np.zeros((len(azth_rng) - 1, len(lat_rng) - 1), dtype=np.float32) + + for i in xrange(len(azth_rng) - 1): + for j in xrange(len(lat_rng) - 1): + rng_azth = (np.deg2rad(azth_rng[i]), np.deg2rad(azth_rng[i+1])) + rng_lat = (np.deg2rad(lat_rng[j]), np.deg2rad(lat_rng[j+1])) + + idxs_azth = np.where((azth_ref >= rng_azth[0]) & (azth_ref < rng_azth[1])) + idxs_lat = np.where((lat_ref >= rng_lat[0]) & (lat_ref < rng_lat[1])) + + idxs = np.intersect1d(idxs_azth, idxs_lat) + + if len(idxs): + mu_azth, mu_lat = np.rad2deg(np.mean(azth_err[idxs])), \ + np.rad2deg(np.mean(lat_err[idxs])) + + out[i, j] = mu_azth + mu_lat + else: + out[i, j] = 0.0 + + return out + + +def plot_polar_heatmap(data, name, interp_factor=5., color_limits=False, + hide_colorbar=False, vmin=None, vmax=None, log_scale=True, + dpi=200, output_dir=None): + """Plots the polar heatmap describing azimuth and latitude / elevation components. + + Plots the polar heatmap where each cell of the heatmap corresponds to + the specific element of the array provided by `gather_polar_errors` + function. + + Parameters + ---------- + data : 2D array + Indicates the array containing the sum of angular errors within the + specified angular ranges. It is usually provided by `gather_polar_errors` + function. + + name : str + Indicates the name of the output png file. + + interp_factor : float + Indicates the interpolation factor of the heatmap. + + color_limits : boolean + Specifies if the determined intensity limits should be returned. + + hide_colorbar : boolean + Specifies if the colorbar should be hidden. + + vmin : float + Indicates the minimum value of the colorbar. + + vmax : float + Indicates the maximum value of the colorbar. + + log_scale : float + Specifies if the heatmap sould be in the logarithmic scale. + + dpi : integer + Indicates the DPI of the output image. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + """ + th0, th1 = 0., 180. + r0, r1 = 0, 90 + thlabel, rlabel = 'Azimuth', 'Elevation' + + tr_scale = Affine2D().scale(np.pi/180., 1.) + tr = tr_scale + PolarAxes.PolarTransform() + + lat_ticks = [(.0*90., '0$^{\circ}$'), + (.33*90., '30$^{\circ}$'), + (.66*90., '60$^{\circ}$'), + (1.*90., '90$^{\circ}$')] + r_grid_locator = FixedLocator([v for v, s in lat_ticks]) + r_grid_formatter = DictFormatter(dict(lat_ticks)) + + angle_ticks = [(0*180., '90$^{\circ}$'), + (.25*180., '45$^{\circ}$'), + (.5*180., '0$^{\circ}$'), + (.75*180., '-45$^{\circ}$'), + (1.*180., '-90$^{\circ}$')] + theta_grid_locator = FixedLocator([v for v, s in angle_ticks]) + theta_tick_formatter = DictFormatter(dict(angle_ticks)) + + grid_helper = GridHelperCurveLinear(tr, + extremes=(th0, th1, r0, r1), + grid_locator1=theta_grid_locator, + grid_locator2=r_grid_locator, + tick_formatter1=theta_tick_formatter, + tick_formatter2=r_grid_formatter) + + fig = plt.figure() + ax = floating_axes.FloatingSubplot(fig, 111, grid_helper=grid_helper) + fig.add_subplot(ax) + ax.set_facecolor('white') + + ax.axis["bottom"].set_visible(False) + ax.axis["top"].toggle(ticklabels=True, label=True) + ax.axis["top"].set_axis_direction("bottom") + ax.axis["top"].major_ticklabels.set_axis_direction("top") + ax.axis["top"].label.set_axis_direction("top") + + ax.axis["left"].set_axis_direction("bottom") + ax.axis["right"].set_axis_direction("top") + + ax.axis["top"].label.set_text(thlabel) + ax.axis["left"].label.set_text(rlabel) + + aux_ax = ax.get_aux_axes(tr) + aux_ax.patch = ax.patch + ax.patch.zorder = 0.9 + + rad = np.linspace(0, 90, data.shape[1]) + azm = np.linspace(0, 180, data.shape[0]) + + f = interpolate.interp2d(rad, azm, data, kind='linear', bounds_error=True, fill_value=0) + + new_rad = np.linspace(0, 90, 180*interp_factor) + new_azm = np.linspace(0, 180, 360*interp_factor) + new_data_angle_dist = f(new_rad, new_azm) + new_r, new_th = np.meshgrid(new_rad, new_azm) + new_data_angle_dist += 1. + + if log_scale: + data_mesh = aux_ax.pcolormesh(new_th, new_r, new_data_angle_dist, cmap='jet', + norm=colors.LogNorm(vmin=1. if vmin is None else vmin, + vmax=new_data_angle_dist.max() if vmax is None else vmax)) + else: + data_mesh = aux_ax.pcolormesh(new_th, new_r, new_data_angle_dist, cmap='jet', vmin=vmin, vmax=vmax) + + cbar = plt.colorbar(data_mesh, orientation='vertical', shrink=.88, pad=.1, aspect=15) + cbar.ax.set_ylabel('Absolute error, [deg.]') + + if hide_colorbar: + cbar.remove() + + ax.grid(False) + + plt.show() + + if output_dir is not None: + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + fig.savefig(os.path.join(output_dir, '{}_chart.png'.format(name)), + transparent=False, + bbox_inches='tight', + pad_inches=0.1, + dpi=dpi) + + if color_limits: + return 1., new_data_angle_dist.max() + + +def plot_histogram_fig(data, num_bins, xticks, color, splot_index=111, output_dir=None, + xlim=(None, None), ylim=(None, None), name=None, in_percent=False, + bar_width=0.8, ticks_pad=7, xticks_fontsize=22, + yticks_fontsize=22, xlabel=None, ylabel=None, + labels_fontsize=20, grid_alpha=0.3, title_fontsize=22, + exp_fontsize=15, type=None, figsize=(12, 8), dpi=200): + """Plots the histogram from a given data. + + Parameters + ---------- + data : 1D array + Indicates the array containing the values. + + num_bins : integer + Indicates the number of histogram bins. + + xticks : array + Indicates the array of ticks of the X-axis. + + color : str + Indicates the color of the histogram. + + Returns + ------- + tuple : tuple of Figure and Axis objects. + """ + fig = plt.figure(figsize=figsize) + ax = fig.add_subplot(splot_index) + + weights = np.ones_like(data)/float(len(data)) + + n, bins, patches = ax.hist(data, num_bins, color=color, + rwidth=bar_width, weights=weights) + + ax.tick_params(axis='x', + labelsize=xticks_fontsize, + colors='#000000', + which='both', + direction='out', + length=6, + width=2, + pad=ticks_pad) + ax.tick_params(axis='x', + labelsize=xticks_fontsize, + colors='#000000', + which='minor', + direction='out', + length=4, + width=2, + pad=ticks_pad) + + ax.tick_params(axis='y', + labelsize=yticks_fontsize, + colors='#000000', + which='major', + direction='out', + length=6, + width=2, + pad=ticks_pad) + ax.tick_params(axis='y', + labelsize=xticks_fontsize, + colors='#000000', + which='minor', + direction='out', + length=4, + width=2, + pad=ticks_pad) + + ax.xaxis.set_minor_locator(ticker.MultipleLocator(5)) + ax.yaxis.set_minor_locator(ticker.MultipleLocator(5)) + + plt.xticks(xticks, xticks) + + ax.xaxis.offsetText.set_fontsize(exp_fontsize) + ax.yaxis.offsetText.set_fontsize(exp_fontsize) + + ax.ticklabel_format(axis='y', style='sci', scilimits=(-2, 2), + fontsize=yticks_fontsize) + ax.yaxis.set_major_locator(ticker.MaxNLocator(nbins=6, min_n_ticks=6)) + + ax.xaxis.grid(False) + ax.yaxis.grid(False) + + ax.set_ylabel(ylabel, labelpad=2, fontsize=labels_fontsize, color='black') + ax.set_xlabel(xlabel, labelpad=2, fontsize=labels_fontsize, color='black') + + ax.set_ylim(ylim) + ax.set_xlim(xlim) + + if in_percent: + vals = ax.get_yticks() + ax.set_yticklabels(['{:3.0f}'.format(x*100) for x in vals]) + + ax.set_axisbelow(True) + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + plt.tight_layout() + fig.savefig(os.path.join(output_dir, '{}_chart.png'.format(name)), + transparent=False, bbox_inches='tight', pad_inches=0.1, dpi=dpi) + + return fig, ax + + +def create_pie_chart(data, rngs, colors=['#244268', '#426084', '#67809F', '#95A9C1', '#C6D2E0'], + unit_scale=1.0, measure_quantity='m^3', figsize=(33, 15), + legend_loc=(0.383, -0.25), zebra_color=(False, 3), + legend_fontsize=50, chart_fontsize=60, dpi=72, name=None, + output_dir=None): + """Plots the piechart of from a given data. + + Parameters + ---------- + data : 1D array + Indicates the array containing the values. + + rngs : tuple of tuples + Indicates the ranges of the piechart. + + colors : array + Indicates the color for the region of the piechart corresponding to the + specific range. + + unit_scale : float + Indicates the scale factor of the data values. + + measure_quantity : str + Indicates the name of measure of the values. + + figsize : tuple of integers + Indicates the size of the output figure. + + legend_loc : tuple + Indicates the position of the legend of the figure. + + zebra_color : tuple + Allows to change the text color of the region to white from the first to + the speficied index of the region (True, reg_index). + + legend_fontsize : integer + Indicates the fontsize of the legend. + + chart_fontsize : integer + Indicates the fontsize of the figure. + + dpi : integer + Indicates the DPI of the output image. + + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + """ + def number(val): + if val < 1000: + return '%d' % val + + sv = str(val) + return '$\mathregular{10^{%d}}$' % (len(sv)-2) if val % 10 == 0 else '%0.0e' % val + + def get_title(v1, v2, measure_quantity): + ftm = '%s $\minus$ %s %s' + return ftm % (number(v1), number(v2), measure_quantity) + + data_ranges = [] + df = data * unit_scale + for rng in rngs: + rng_min, rng_max = rng[0], rng[1] + data_rng = df[(df > rng_min) & (df < rng_max)] + data_ranges.append(data_rng) + + num_elem = [len(p) for p in data_ranges] + se = sum(num_elem) + + print('Num of particles: {}'%(se)) + + proc_particles = [n/float(se) * 100.0 for n in num_elem] + + for size, rng in zip(num_elem, rngs): + print('{}-{}: {}'%(rng[0], rng[1], size)) + + titles = [get_title(minv, maxv, measure_quantity) for minv, maxv in rngs] + + textprops = {'fontsize': chart_fontsize, + 'weight': 'normal', + 'family': 'sans-serif'} + pie_width = 0.5 + fig, ax = plt.subplots(figsize=figsize) + ax.axis('equal') + patches, texts, autotexts = ax.pie(proc_particles, + textprops=textprops, + colors=colors, + autopct='%1.1f%%', + radius=1, + pctdistance=1-pie_width/2) + + if (zebra_color is not None) and (zebra_color[0]): + for tt in autotexts[:zebra_color[1]]: + tt.set_color('white') + + plt.setp(patches, + width=pie_width, + edgecolor='white') + + plt.legend(patches, titles, loc=legend_loc, fontsize=legend_fontsize) + + _d, _offset, _di = [1, -1], [0.45, 0.45], 0 + + for t, p in zip(autotexts, proc_particles): + if p < 2.0: + pos = list(t.get_position()) + pos[0] = pos[0] + _d[_di] * _offset[_di] + + t.set_position(pos) + _di += 1 + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + plt.tight_layout() + fig.savefig(os.path.join(output_dir, '{}_chart.png'.format(name)), + bbox_inches='tight', transparent=True, pad_inches=0.1, dpi=dpi) + + +def _bbox_3D(img): + """Crops the non-zero part of a volume + """ + r = np.any(img, axis=(1, 2)) + c = np.any(img, axis=(0, 2)) + z = np.any(img, axis=(0, 1)) + + rmin, rmax = np.where(r)[0][[0, -1]] + cmin, cmax = np.where(c)[0][[0, -1]] + zmin, zmax = np.where(z)[0][[0, -1]] + + return rmin, rmax, cmin, cmax, zmin, zmax + + +def plot_3d_orientation_map(name, lat_data, azth_data, radius_structure_elem=1, + output_dir=None, width=512, height=512, + camera_azth=44.5, camera_elev=35.8, camera_roll=0.0, + camera_fov=35.0, camera_zoom=0.0035, + camera_loc=(67.0, 81.6, 45.2), xlabel='', ylabel='', + zlabel='', axis_color='w', background_color='k'): + """Renders orientation data in 3D with RGB angular color-coding. + + Parameters + ---------- + name : str + Indicates the name of the output png file. + + lat_data : 3D array + Indicates the 3D array containing latitude / elevation angle at every point of + the skeleton in radians. + + azth_data : 3D array + Indicates the 3D array containing azimuth angle at every point of the skeleton + in radians. + + radius_structure_elem : integer + Indicates the size of the structure element of the dilation process to + thicken the skeleton. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + + width : int + Indicates the width of the visualization window. + + height : int + Indicates the width of the visualization window. + + camera_azth : float + Indicates the azimuth angle of the camera. + + camera_elev : float + Indicates the latitude / elevation angle of the camera. + + camera_roll : float + Indicates the roll angle of the camera. + + camera_fov : float + Indicates the field of view of the camera. + + camera_zoom : float + Indicates the zoom level of the camera. + + camera_loc : tuple + Indicates the camera location. + + xlabel : str + Indicates the label along the x-axis. + + ylabel : str + Indicates the label along the y-axis. + + zlabel : str + Indicates the label along the z-axis. + + axis_color : str + Indicates the color of axes. + + background_color : str + Indicates the background color of the figure. + """ + if not visvis_available: + print('The visvis package is not found. The visualization cannot be done.') + return + + rmin, rmax, cmin, cmax, zmin, zmax = _bbox_3D(azth_data) + + azth, lat = azth_data[rmin:rmax, cmin:cmax, zmin:zmax], \ + np.abs(lat_data[rmin:rmax, cmin:cmax, zmin:zmax]) + + skel = azth.copy().astype(np.float32) + skel[skel.nonzero()] = 1. + + azth = ndi.grey_dilation(azth, structure=morphology.ball(radius_structure_elem)) + lat = ndi.grey_dilation(lat, structure=morphology.ball(radius_structure_elem)) + skel = ndi.binary_dilation(skel, structure=morphology.ball(radius_structure_elem)) + + Z, Y, X = skel.nonzero() + vol_orient = np.zeros(skel.shape + (3,), dtype=np.float32) + + print(vol_orient.size, vol_orient[skel.nonzero()].size) + + for z, y, x in zip(Z, Y, X): + vol_orient[z, y, x] = geo2rgb(lat[z, y, x], azth[z, y, x]) + + app = vv.use() + + fig = vv.figure() + fig._currentAxes = None + fig.relativeFontSize = 2. + fig.position.w = width + fig.position.h = height + + t = vv.volshow(vol_orient[:, :, :], renderStyle='iso') + t.isoThreshold = 0.5 + + a = vv.gca() + a.camera.azimuth = camera_azth + a.camera.elevation = camera_elev + a.camera.roll = camera_roll + a.camera.fov = camera_fov + a.camera.zoom = camera_zoom + a.camera.loc = camera_loc + + a.bgcolor = background_color + a.axis.axisColor = axis_color + a.axis.xLabel = xlabel + a.axis.yLabel = ylabel + a.axis.zLabel = zlabel + + # def mouseUp(event): + # print 'mouseUp!!' + # a = vv.gca() + # print a.camera.GetViewParams() + # + # a.eventMouseUp.Bind(mouseUp) + # fig.eventMouseUp.Bind(mouseUp) + # + # a.Draw() + # fig.DrawNow() + + if output_dir is not None: + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + vv.screenshot(os.path.join(output_dir, '{}_3d_orientation.png'.format(name)), + sf=1, bg=background_color) + app.Run() + + +def plot_3d_diameter_map(name, data, unit_scale=1.0, measure_quantity='vox', + radius_structure_elem=1, output_dir=None, width=512, + height=512, camera_azth=44.5, camera_elev=35.8, + camera_roll=0.0, camera_fov=35.0, camera_zoom=0.0035, + camera_loc=(67.0, 81.6, 45.2), xlabel='', ylabel='', + zlabel='', axis_color='w', background_color='k', + cb_x_offset=10): + """Renders orientation data in 3D with RGB angular color-coding. + + Parameters + ---------- + name : str + Indicates the name of the output png file. + + data : 3D array + Indicates the 3D array containing diameter at every point of the skeleton. + + unit_scale : float + Indicates the scale factor of the data values. + + measure_quantity : str + Indicates the name of measure of the values. + + radius_structure_elem : integer + Indicates the size of the structure element of the dilation process to + thicken the skeleton. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + + camera_azth : float + Indicates the azimuth angle of the camera. + + width : int + Indicates the width of the visualization window. + + height : int + Indicates the width of the visualization window. + + camera_elev : float + Indicates the latitude / elevation angle of the camera. + + camera_roll : float + Indicates the roll angle of the camera. + + camera_fov : float + Indicates the field of view of the camera. + + camera_zoom : float + Indicates the zoom level of the camera. + + camera_loc : tuple + Indicates the camera location. + + xlabel : str + Indicates the label along the x-axis. + + ylabel : str + Indicates the label along the y-axis. + + zlabel : str + Indicates the label along the z-axis. + + axis_color : str + Indicates the color of axes. + + background_color : str + Indicates the background color of the figure. + + cb_x_offset : int + Indicates the offset of the colorbar from the right window side. + """ + if not visvis_available: + print('The visvis package is not found. The visualization cannot be done.') + return + + rmin, rmax, cmin, cmax, zmin, zmax = _bbox_3D(data) + dmtr = data[rmin:rmax, cmin:cmax, zmin:zmax] * unit_scale + skel = np.zeros_like(dmtr, dtype=np.uint8) + skel[dmtr.nonzero()] = 1 + + dmtr = ndi.grey_dilation(dmtr, structure=morphology.ball(radius_structure_elem)) + skel = ndi.binary_dilation(skel, structure=morphology.ball(radius_structure_elem)).astype(np.float32) + skel[skel.nonzero()] = 1. + + dmtr = dmtr * skel + + app = vv.use() + + fig = vv.figure() + fig._currentAxes = None + fig.relativeFontSize = 2. + fig.position.w = width + fig.position.h = height + + t = vv.volshow(dmtr[:, :, :], renderStyle='iso') + t.isoThreshold = 0.5 + t.colormap = vv.CM_JET + + a = vv.gca() + a.camera.azimuth = camera_azth + a.camera.elevation = camera_elev + a.camera.roll = camera_roll + a.camera.fov = camera_fov + a.camera.zoom = camera_zoom + a.camera.loc = camera_loc + + a.bgcolor = background_color + a.axis.axisColor = axis_color + a.axis.xLabel = xlabel + a.axis.yLabel = ylabel + a.axis.zLabel = zlabel + + cb = vv.colorbar() + cb.SetLabel('Diameter, [{}]'.format(measure_quantity)) + cb._label.position.x += cb_x_offset + + if output_dir is not None: + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + vv.screenshot(os.path.join(output_dir, '{}_3d_diameter.png'.format(name)), + sf=1, bg='w') + app.Run() + + +def plot_color_wheel(name, output_dir=None, dpi=500, xlabel='Elevation', + ylabel='Azimuth', fontsize=10, num_xticks=4, yticks=(-90, 90)): + """Plots the color wheel for visualizations of 3D orintation. + + Parameters + ---------- + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + + dpi : integer + Indicates the DPI of the output image. + + xlabel : str + Indicates the text along the x-axis. + + ylabel : str + Indicates the text along the y-axis. + + fontsize : int + Indicates the font size of labels along axes. + + num_xticks : int + Indicates the number of ticks along axes. + + yticks : tuple + Indicates the range of minimum and maximum values along the y-axis. + """ + + azth, lat = np.linspace(0., 1., num=180), np.linspace(0., 1., num=90) + rgb_arr = np.zeros((len(azth), len(lat), 3)) + for i in xrange(len(azth)): + for j in xrange(len(lat)): + rgb_arr[i, j, :] = colors.hsv_to_rgb([azth[i], lat[j], 1.0]) + + fig, ax = plt.subplots(figsize=(2, 2)) + ax.set_facecolor('red') + + ax.set_xlim([0, 90]) + ax.set_ylim([0, 181]) + + ax.set_yticks(np.linspace(0, 180, num=7).astype(np.int32)) + ax.set_yticklabels(np.linspace(yticks[0], yticks[1], num=7).astype(np.int32)) + + ax.set_xlabel(xlabel, fontsize=fontsize, labelpad=0, color='w') + ax.set_ylabel(ylabel, fontsize=fontsize, labelpad=0, color='w') + + xmajorlocator = ticker.LinearLocator(num_xticks) + ax.xaxis.set_major_locator(xmajorlocator) + ax.tick_params(direction='out', length=2, width=0, labelsize=fontsize, pad=0, colors='w') + + ax.imshow(rgb_arr) + + if output_dir is not None: + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + fig.savefig(os.path.join(output_dir, '{}_color_bar.png'.format(name)), + transparent=True, + bbox_inches='tight', + pad_inches=0.1, + dpi=dpi) + plt.show() + +def plot_fourier_orientation(data, orient_blocks, block_shape, figsize=(12,12), + cmap='gray', line_length=20, line_width=2.5, + line_color='red', line_style='-', name=None, + output_dir=None, dpi=200): + """Plots the orientation map in a block-wise manner. + + Plots the orientation vector over the image `data` at the center of each + block of the subdivided image. The orientation at every block 'orient_blocks' + and its size `block_shape` are specified by orientation estimation method. + The result can be and stored as a png file with DPI of `dpi` to the folder + specified by `output_dir`. + + Parameters + ---------- + data : ndarray + An image on top of which, the orientation will be plotted.. + + orient_blocks : ndarray + 2D array of orientation at every block of the subdivided image. + + block_shape : tuple + Indicates the block size within which the orientation is calculated. + + figsize : tuple of integers + Indicates the size of the output figure. + + cmap : str + Indicates the name of a colormap used for image. + + line_length : integer + Indicates the length of the orientation vector at each block. + + line_width : float + Indicates the line width of the orientation vector at each block. + + line_color : str + Indicates the line color of the orientation vector at each block. + + line_style : str + Indicates the line style of the orientation vector at each block. + + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + + dpi : integer + Indicates the DPI of the output image. + """ + fig, ax = plt.subplots(figsize=figsize) + + ax.xaxis.grid(False) + ax.yaxis.grid(False) + ax.set_axis_off() + + ax.imshow(data, cmap=cmap) + + for i in xrange(orient_blocks.shape[0]): + for j in xrange(orient_blocks.shape[1]): + y0, x0 = block_shape[0] * j + block_shape[0]/2, \ + block_shape[1] * i + block_shape[1]/2 + + orientation = orient_blocks[i ,j] + + x2 = x0 - math.sin(orientation) * line_length + y2 = y0 - math.cos(orientation) * line_length + + x3 = x0 + math.sin(orientation) * line_length + y3 = y0 + math.cos(orientation) * line_length + + ax.plot((x2, x3), (y2, y3), linestyle=line_style, + linewidth=line_width, color=line_color) + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + plt.tight_layout() + fig.savefig(os.path.join(output_dir, '{}_fourier_orientation_map.png'.format(name)), + bbox_inches='tight', transparent=True, pad_inches=0.05, dpi=dpi) + + plt.show() diff --git a/quanfima/quanfima/visualizationDE.py b/quanfima/quanfima/visualizationDE.py new file mode 100644 index 0000000..168efe6 --- /dev/null +++ b/quanfima/quanfima/visualizationDE.py @@ -0,0 +1,1081 @@ +from __future__ import print_function +import os +import math +import numpy as np +import matplotlib as mpl +from skimage import io, filters +import matplotlib.colors as colors +import matplotlib.pyplot as plt +import matplotlib.ticker as ticker +import mpl_toolkits.axisartist.floating_axes as floating_axes +from matplotlib import cm +from scipy import ndimage as ndi, interpolate +from skimage import morphology +from matplotlib.projections import PolarAxes +from matplotlib.transforms import Affine2D +from mpl_toolkits.axisartist.grid_finder import DictFormatter, FixedLocator +from mpl_toolkits.axes_grid1 import make_axes_locatable +from mpl_toolkits.axisartist.floating_axes import GridHelperCurveLinear +from quanfima.quanfima.utils2 import geo2rgb +from quanfima.quanfima import visvis_available +if visvis_available: + import visvis as vv + + +def plot_orientation_map(orient_map, fiber_skel, radius_structure_elem=1, + figsize=(20, 20), cmap='hsv', dpi=200, min_label=u'0°', + max_label=u'180°', name=None, output_dir=None): + """Plots the orientation map with the color wheel. + + Plots the orientation map from the provided angles `orient_map` and skeleton + `fiber_skel` of size `figsize` using the colormap `cmap` and writes as a png + file with DPI of `dpi` to the folder specified by `output_dir`. + + Parameters + ---------- + orient_map : ndarray + 2D array of orientation at every point of the skeleton. + + fiber_skel : ndarray + The binary skeleton extracted from the binary data. + + radius_structure_elem : integer + Indicates the size of the structure element of the dilation process to + thicken the skeleton. + + figsize : tuple of integers + Indicates the size of the output figure. + + cmap : str + Indicates the name of a colormap used to map angles to colors. + + dpi : integer + Indicates the DPI of the output image. + + min_label : str + Indicates the label of minimum degree. + + max_label : str + Indicates the label of minimum degree. + + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + """ + disk = morphology.disk(radius_structure_elem) + orient_map = ndi.grey_dilation(orient_map, structure=disk).astype(np.float32) + fiber_skel = ndi.binary_dilation(fiber_skel, structure=disk).astype(np.float32) + + fig = plt.figure(figsize=figsize) + masked_orient_map = np.ma.masked_where(fiber_skel == 0, orient_map) + + quant_steps = 2056 + cmap_obj = cm.get_cmap(cmap, quant_steps) + cmap_obj.set_bad(color='black') + + omap_ax = fig.add_axes([0.0, 0.0, 1.0, 1.0]) + omap_ax.set_axis_off() + omap_ax.imshow(masked_orient_map, cmap=cmap_obj, interpolation=None) + + #display_axes = fig.add_axes([0.780, -0.076, 0.2, 0.2], projection='polar') + #display_axes._direction = np.pi + + norm = mpl.colors.Normalize(0.0, np.pi) + #cb = mpl.colorbar.ColorbarBase(display_axes, cmap=cmap_obj, norm=norm, orientation='horizontal') + + #display_axes.text(0.09, 0.56, min_label, color='white', fontsize=20, weight='bold', horizontalalignment='center', + # verticalalignment='center', transform=display_axes.transAxes) + + #display_axes.text(0.85, 0.56, max_label, color='white', fontsize=20, weight='bold', horizontalalignment='center', + # verticalalignment='center', transform=display_axes.transAxes) + + #cb.outline.set_visible(False) + #display_axes.set_axis_off() + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + fig.savefig(os.path.join(output_dir, 'orientation_map.png'), + transparent=False, + bbox_inches='tight', + pad_inches=0.1, + dpi=dpi) + plt.show() + + +def plot_diameter_map(thickness_map, fiber_skel, radius_structure_elem=1, + figsize=(15, 15), cmap='hsv', tmin=None, tmax=None, dpi=200, + labelsize=20, label='Diameter, [pixels]', name=None, + output_dir=None): + """Plots the diameter map with the colorbar. + + Plots the diameter map from the provided diameters `thickness_map` and + skeleton `fiber_skel` of size `figsize` using the colormap `cmap` and the + limits of the colorbar specified by `tmin` and `tmax`, and writes as a png + file with DPI of `dpi` to the folder specified by `output_dir`. + + Parameters + ---------- + thickness_map : ndarray + 2D array of diameter at every point of the skeleton. + + fiber_skel : ndarray + The binary skeleton extracted from the binary data. + + radius_structure_elem : integer + Indicates the size of the structure element of the dilation process to + thicken the skeleton. + + figsize : tuple of integers + Indicates the size of the output figure. + + cmap : str + Indicates the name of a colormap used to map angles to colors. + + tmin : float + Indicates the minimum value of the colorbar. + + tmax : float + Indicates the maximum value of the colorbar. + + dpi : integer + Indicates the DPI of the output image. + + labelsize : integer + Indicates the fontsize of the label of the colorbar. + + label : str + Indicates the label of the colorbar. + + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + """ + disk = morphology.disk(radius_structure_elem) + thickness_map = ndi.grey_dilation(thickness_map, structure=disk).astype(np.float32) + fiber_skel = ndi.binary_dilation(fiber_skel, structure=disk).astype(np.float32) + + masked_thickness_map = np.ma.masked_where(fiber_skel == 0, thickness_map) + + cmap_obj = cm.get_cmap(cmap) + cmap_obj.set_bad(color='black') + + fig = plt.figure(figsize=figsize) + ax = plt.subplot(111) + ax.set_axis_off() + + im = ax.imshow(masked_thickness_map, cmap=cmap_obj, vmin=tmin, vmax=tmax) + divider = make_axes_locatable(ax) + cax = divider.append_axes('right', size="2.5%", pad=0.05) + cbar = plt.colorbar(im, cax=cax) + cbar.ax.tick_params(labelsize=labelsize) + cbar.set_label(label, fontsize=labelsize) + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + fig.savefig(os.path.join(output_dir, '{}_diameter_map.png'.format(name)), + transparent=False, + bbox_inches='tight', + pad_inches=0.1, + dpi=dpi) + plt.show() + + +def gather_polar_errors(datasets_path, estimation_path, + azth_rng=np.arange(-90, 91, step=3), + lat_rng=np.arange(0, 91, step=3)): + """Computes the absolute angular error in ranges between estimated datasets. + + Calculates the absolute angular error in ranges `azth_rng` and `lat_rng` + between estimated orientation datasets localed at paths `datasets_path` and + `estimation_path`. + + Parameters + ---------- + datasets_path : str + Indicates the path to the reference / estimated dataset. + + estimation_path : str + Indicates the path to the estimated / reference dataset. + + azth_rng : array + Indicates the ranges of azimuth angles where the error is accumulated. + + lat_rng : array + Indicates the ranges of latitude or elevation angles where the error + is accumulated. + + Returns + ------- + array : 2D array + The 2D array of accumulated errors within the combinations of the specified + angular ranges. + """ + reference_dataset = np.load(datasets_path).item() + estimated_dataset = np.load(estimation_path).item() + idxs = estimated_dataset['indices'] + + print(reference_dataset.keys()) + print(estimated_dataset.keys()) + + ref_azth, ref_lat = reference_dataset['azth'][idxs], reference_dataset['lat'][idxs] + res_azth, res_lat = estimated_dataset['azth'][idxs], estimated_dataset['lat'][idxs] + + abs_err_azth, abs_err_lat = np.abs(ref_azth - res_azth), np.abs(ref_lat - res_lat) + + out = _angle_err2mean_abs_err(abs_err_azth, abs_err_lat, + ref_azth, ref_lat, + azth_rng=azth_rng, lat_rng=lat_rng) + + return out + + +def _angle_err2mean_abs_err(azth_err, lat_err, azth_ref, lat_ref, + azth_rng=np.arange(-90, 91, step=5), + lat_rng=np.arange(0, 91, step=5)): + """Computes the absolute angular error in ranges between estimated datasets. + + Parameters + ---------- + azth_err : 3D array + Indicates the error of azimuth components. + + lat_err : 3D array + Indicates the error of latitude / elevation components. + + azth_ref : 3D array + Indicates the reference dataset of estimated azimuth component. + + lat_ref : 3D array + Indicates the reference dataset of estimated latitude or elevation component. + + azth_rng : array + Indicates the ranges of azimuth angles where the error is accumulated. + + lat_rng : array + Indicates the ranges of latitude and elevation angles where the error + is accumulated. + + Returns + ------- + array : 2D array + The 2D array of accumulated errors within the combinations of the specified + angular ranges. + """ + out = np.zeros((len(azth_rng) - 1, len(lat_rng) - 1), dtype=np.float32) + + for i in xrange(len(azth_rng) - 1): + for j in xrange(len(lat_rng) - 1): + rng_azth = (np.deg2rad(azth_rng[i]), np.deg2rad(azth_rng[i+1])) + rng_lat = (np.deg2rad(lat_rng[j]), np.deg2rad(lat_rng[j+1])) + + idxs_azth = np.where((azth_ref >= rng_azth[0]) & (azth_ref < rng_azth[1])) + idxs_lat = np.where((lat_ref >= rng_lat[0]) & (lat_ref < rng_lat[1])) + + idxs = np.intersect1d(idxs_azth, idxs_lat) + + if len(idxs): + mu_azth, mu_lat = np.rad2deg(np.mean(azth_err[idxs])), \ + np.rad2deg(np.mean(lat_err[idxs])) + + out[i, j] = mu_azth + mu_lat + else: + out[i, j] = 0.0 + + return out + + +def plot_polar_heatmap(data, name, interp_factor=5., color_limits=False, + hide_colorbar=False, vmin=None, vmax=None, log_scale=True, + dpi=200, output_dir=None): + """Plots the polar heatmap describing azimuth and latitude / elevation components. + + Plots the polar heatmap where each cell of the heatmap corresponds to + the specific element of the array provided by `gather_polar_errors` + function. + + Parameters + ---------- + data : 2D array + Indicates the array containing the sum of angular errors within the + specified angular ranges. It is usually provided by `gather_polar_errors` + function. + + name : str + Indicates the name of the output png file. + + interp_factor : float + Indicates the interpolation factor of the heatmap. + + color_limits : boolean + Specifies if the determined intensity limits should be returned. + + hide_colorbar : boolean + Specifies if the colorbar should be hidden. + + vmin : float + Indicates the minimum value of the colorbar. + + vmax : float + Indicates the maximum value of the colorbar. + + log_scale : float + Specifies if the heatmap sould be in the logarithmic scale. + + dpi : integer + Indicates the DPI of the output image. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + """ + th0, th1 = 0., 180. + r0, r1 = 0, 90 + thlabel, rlabel = 'Azimuth', 'Elevation' + + tr_scale = Affine2D().scale(np.pi/180., 1.) + tr = tr_scale + PolarAxes.PolarTransform() + + lat_ticks = [(.0*90., '0$^{\circ}$'), + (.33*90., '30$^{\circ}$'), + (.66*90., '60$^{\circ}$'), + (1.*90., '90$^{\circ}$')] + r_grid_locator = FixedLocator([v for v, s in lat_ticks]) + r_grid_formatter = DictFormatter(dict(lat_ticks)) + + angle_ticks = [(0*180., '90$^{\circ}$'), + (.25*180., '45$^{\circ}$'), + (.5*180., '0$^{\circ}$'), + (.75*180., '-45$^{\circ}$'), + (1.*180., '-90$^{\circ}$')] + theta_grid_locator = FixedLocator([v for v, s in angle_ticks]) + theta_tick_formatter = DictFormatter(dict(angle_ticks)) + + grid_helper = GridHelperCurveLinear(tr, + extremes=(th0, th1, r0, r1), + grid_locator1=theta_grid_locator, + grid_locator2=r_grid_locator, + tick_formatter1=theta_tick_formatter, + tick_formatter2=r_grid_formatter) + + fig = plt.figure() + ax = floating_axes.FloatingSubplot(fig, 111, grid_helper=grid_helper) + fig.add_subplot(ax) + ax.set_facecolor('white') + + ax.axis["bottom"].set_visible(False) + ax.axis["top"].toggle(ticklabels=True, label=True) + ax.axis["top"].set_axis_direction("bottom") + ax.axis["top"].major_ticklabels.set_axis_direction("top") + ax.axis["top"].label.set_axis_direction("top") + + ax.axis["left"].set_axis_direction("bottom") + ax.axis["right"].set_axis_direction("top") + + ax.axis["top"].label.set_text(thlabel) + ax.axis["left"].label.set_text(rlabel) + + aux_ax = ax.get_aux_axes(tr) + aux_ax.patch = ax.patch + ax.patch.zorder = 0.9 + + rad = np.linspace(0, 90, data.shape[1]) + azm = np.linspace(0, 180, data.shape[0]) + + f = interpolate.interp2d(rad, azm, data, kind='linear', bounds_error=True, fill_value=0) + + new_rad = np.linspace(0, 90, 180*interp_factor) + new_azm = np.linspace(0, 180, 360*interp_factor) + new_data_angle_dist = f(new_rad, new_azm) + new_r, new_th = np.meshgrid(new_rad, new_azm) + new_data_angle_dist += 1. + + if log_scale: + data_mesh = aux_ax.pcolormesh(new_th, new_r, new_data_angle_dist, cmap='jet', + norm=colors.LogNorm(vmin=1. if vmin is None else vmin, + vmax=new_data_angle_dist.max() if vmax is None else vmax)) + else: + data_mesh = aux_ax.pcolormesh(new_th, new_r, new_data_angle_dist, cmap='jet', vmin=vmin, vmax=vmax) + + cbar = plt.colorbar(data_mesh, orientation='vertical', shrink=.88, pad=.1, aspect=15) + cbar.ax.set_ylabel('Absolute error, [deg.]') + + if hide_colorbar: + cbar.remove() + + ax.grid(False) + + plt.show() + + if output_dir is not None: + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + fig.savefig(os.path.join(output_dir, '{}_chart.png'.format(name)), + transparent=False, + bbox_inches='tight', + pad_inches=0.1, + dpi=dpi) + + if color_limits: + return 1., new_data_angle_dist.max() + + +def plot_histogram_fig(data, num_bins, xticks, color, splot_index=111, output_dir=None, + xlim=(None, None), ylim=(None, None), name=None, in_percent=False, + bar_width=0.8, ticks_pad=7, xticks_fontsize=22, + yticks_fontsize=22, xlabel=None, ylabel=None, + labels_fontsize=20, grid_alpha=0.3, title_fontsize=22, + exp_fontsize=15, type=None, figsize=(12, 8), dpi=200): + """Plots the histogram from a given data. + + Parameters + ---------- + data : 1D array + Indicates the array containing the values. + + num_bins : integer + Indicates the number of histogram bins. + + xticks : array + Indicates the array of ticks of the X-axis. + + color : str + Indicates the color of the histogram. + + Returns + ------- + tuple : tuple of Figure and Axis objects. + """ + fig = plt.figure(figsize=figsize) + ax = fig.add_subplot(splot_index) + + weights = np.ones_like(data)/float(len(data)) + + n, bins, patches = ax.hist(data, num_bins, color=color, + rwidth=bar_width, weights=weights) + + ax.tick_params(axis='x', + labelsize=xticks_fontsize, + colors='#000000', + which='both', + direction='out', + length=6, + width=2, + pad=ticks_pad) + ax.tick_params(axis='x', + labelsize=xticks_fontsize, + colors='#000000', + which='minor', + direction='out', + length=4, + width=2, + pad=ticks_pad) + + ax.tick_params(axis='y', + labelsize=yticks_fontsize, + colors='#000000', + which='major', + direction='out', + length=6, + width=2, + pad=ticks_pad) + ax.tick_params(axis='y', + labelsize=xticks_fontsize, + colors='#000000', + which='minor', + direction='out', + length=4, + width=2, + pad=ticks_pad) + + ax.xaxis.set_minor_locator(ticker.MultipleLocator(5)) + ax.yaxis.set_minor_locator(ticker.MultipleLocator(5)) + + plt.xticks(xticks, xticks) + + ax.xaxis.offsetText.set_fontsize(exp_fontsize) + ax.yaxis.offsetText.set_fontsize(exp_fontsize) + + ax.ticklabel_format(axis='y', style='sci', scilimits=(-2, 2), + fontsize=yticks_fontsize) + ax.yaxis.set_major_locator(ticker.MaxNLocator(nbins=6, min_n_ticks=6)) + + ax.xaxis.grid(False) + ax.yaxis.grid(False) + + ax.set_ylabel(ylabel, labelpad=2, fontsize=labels_fontsize, color='black') + ax.set_xlabel(xlabel, labelpad=2, fontsize=labels_fontsize, color='black') + + ax.set_ylim(ylim) + ax.set_xlim(xlim) + + if in_percent: + vals = ax.get_yticks() + ax.set_yticklabels(['{:3.0f}'.format(x*100) for x in vals]) + + ax.set_axisbelow(True) + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + plt.tight_layout() + fig.savefig(os.path.join(output_dir, '{}_chart.png'.format(name)), + transparent=False, bbox_inches='tight', pad_inches=0.1, dpi=dpi) + + return fig, ax + + +def create_pie_chart(data, rngs, colors=['#244268', '#426084', '#67809F', '#95A9C1', '#C6D2E0'], + unit_scale=1.0, measure_quantity='m^3', figsize=(33, 15), + legend_loc=(0.383, -0.25), zebra_color=(False, 3), + legend_fontsize=50, chart_fontsize=60, dpi=72, name=None, + output_dir=None): + """Plots the piechart of from a given data. + + Parameters + ---------- + data : 1D array + Indicates the array containing the values. + + rngs : tuple of tuples + Indicates the ranges of the piechart. + + colors : array + Indicates the color for the region of the piechart corresponding to the + specific range. + + unit_scale : float + Indicates the scale factor of the data values. + + measure_quantity : str + Indicates the name of measure of the values. + + figsize : tuple of integers + Indicates the size of the output figure. + + legend_loc : tuple + Indicates the position of the legend of the figure. + + zebra_color : tuple + Allows to change the text color of the region to white from the first to + the speficied index of the region (True, reg_index). + + legend_fontsize : integer + Indicates the fontsize of the legend. + + chart_fontsize : integer + Indicates the fontsize of the figure. + + dpi : integer + Indicates the DPI of the output image. + + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + """ + def number(val): + if val < 1000: + return '%d' % val + + sv = str(val) + return '$\mathregular{10^{%d}}$' % (len(sv)-2) if val % 10 == 0 else '%0.0e' % val + + def get_title(v1, v2, measure_quantity): + ftm = '%s $\minus$ %s %s' + return ftm % (number(v1), number(v2), measure_quantity) + + data_ranges = [] + df = data * unit_scale + for rng in rngs: + rng_min, rng_max = rng[0], rng[1] + data_rng = df[(df > rng_min) & (df < rng_max)] + data_ranges.append(data_rng) + + num_elem = [len(p) for p in data_ranges] + se = sum(num_elem) + + print('Num of particles: {}'%(se)) + + proc_particles = [n/float(se) * 100.0 for n in num_elem] + + for size, rng in zip(num_elem, rngs): + print('{}-{}: {}'%(rng[0], rng[1], size)) + + titles = [get_title(minv, maxv, measure_quantity) for minv, maxv in rngs] + + textprops = {'fontsize': chart_fontsize, + 'weight': 'normal', + 'family': 'sans-serif'} + pie_width = 0.5 + fig, ax = plt.subplots(figsize=figsize) + ax.axis('equal') + patches, texts, autotexts = ax.pie(proc_particles, + textprops=textprops, + colors=colors, + autopct='%1.1f%%', + radius=1, + pctdistance=1-pie_width/2) + + if (zebra_color is not None) and (zebra_color[0]): + for tt in autotexts[:zebra_color[1]]: + tt.set_color('white') + + plt.setp(patches, + width=pie_width, + edgecolor='white') + + plt.legend(patches, titles, loc=legend_loc, fontsize=legend_fontsize) + + _d, _offset, _di = [1, -1], [0.45, 0.45], 0 + + for t, p in zip(autotexts, proc_particles): + if p < 2.0: + pos = list(t.get_position()) + pos[0] = pos[0] + _d[_di] * _offset[_di] + + t.set_position(pos) + _di += 1 + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + plt.tight_layout() + fig.savefig(os.path.join(output_dir, '{}_chart.png'.format(name)), + bbox_inches='tight', transparent=True, pad_inches=0.1, dpi=dpi) + + +def _bbox_3D(img): + """Crops the non-zero part of a volume + """ + r = np.any(img, axis=(1, 2)) + c = np.any(img, axis=(0, 2)) + z = np.any(img, axis=(0, 1)) + + rmin, rmax = np.where(r)[0][[0, -1]] + cmin, cmax = np.where(c)[0][[0, -1]] + zmin, zmax = np.where(z)[0][[0, -1]] + + return rmin, rmax, cmin, cmax, zmin, zmax + + +def plot_3d_orientation_map(name, lat_data, azth_data, radius_structure_elem=1, + output_dir=None, width=512, height=512, + camera_azth=44.5, camera_elev=35.8, camera_roll=0.0, + camera_fov=35.0, camera_zoom=0.0035, + camera_loc=(67.0, 81.6, 45.2), xlabel='', ylabel='', + zlabel='', axis_color='w', background_color='k'): + """Renders orientation data in 3D with RGB angular color-coding. + + Parameters + ---------- + name : str + Indicates the name of the output png file. + + lat_data : 3D array + Indicates the 3D array containing latitude / elevation angle at every point of + the skeleton in radians. + + azth_data : 3D array + Indicates the 3D array containing azimuth angle at every point of the skeleton + in radians. + + radius_structure_elem : integer + Indicates the size of the structure element of the dilation process to + thicken the skeleton. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + + width : int + Indicates the width of the visualization window. + + height : int + Indicates the width of the visualization window. + + camera_azth : float + Indicates the azimuth angle of the camera. + + camera_elev : float + Indicates the latitude / elevation angle of the camera. + + camera_roll : float + Indicates the roll angle of the camera. + + camera_fov : float + Indicates the field of view of the camera. + + camera_zoom : float + Indicates the zoom level of the camera. + + camera_loc : tuple + Indicates the camera location. + + xlabel : str + Indicates the label along the x-axis. + + ylabel : str + Indicates the label along the y-axis. + + zlabel : str + Indicates the label along the z-axis. + + axis_color : str + Indicates the color of axes. + + background_color : str + Indicates the background color of the figure. + """ + if not visvis_available: + print('The visvis package is not found. The visualization cannot be done.') + return + + rmin, rmax, cmin, cmax, zmin, zmax = _bbox_3D(azth_data) + + azth, lat = azth_data[rmin:rmax, cmin:cmax, zmin:zmax], \ + np.abs(lat_data[rmin:rmax, cmin:cmax, zmin:zmax]) + + skel = azth.copy().astype(np.float32) + skel[skel.nonzero()] = 1. + + azth = ndi.grey_dilation(azth, structure=morphology.ball(radius_structure_elem)) + lat = ndi.grey_dilation(lat, structure=morphology.ball(radius_structure_elem)) + skel = ndi.binary_dilation(skel, structure=morphology.ball(radius_structure_elem)) + + Z, Y, X = skel.nonzero() + vol_orient = np.zeros(skel.shape + (3,), dtype=np.float32) + + print(vol_orient.size, vol_orient[skel.nonzero()].size) + + for z, y, x in zip(Z, Y, X): + vol_orient[z, y, x] = geo2rgb(lat[z, y, x], azth[z, y, x]) + + app = vv.use() + + fig = vv.figure() + fig._currentAxes = None + fig.relativeFontSize = 2. + fig.position.w = width + fig.position.h = height + + t = vv.volshow(vol_orient[:, :, :], renderStyle='iso') + t.isoThreshold = 0.5 + + a = vv.gca() + a.camera.azimuth = camera_azth + a.camera.elevation = camera_elev + a.camera.roll = camera_roll + a.camera.fov = camera_fov + a.camera.zoom = camera_zoom + a.camera.loc = camera_loc + + a.bgcolor = background_color + a.axis.axisColor = axis_color + a.axis.xLabel = xlabel + a.axis.yLabel = ylabel + a.axis.zLabel = zlabel + + # def mouseUp(event): + # print 'mouseUp!!' + # a = vv.gca() + # print a.camera.GetViewParams() + # + # a.eventMouseUp.Bind(mouseUp) + # fig.eventMouseUp.Bind(mouseUp) + # + # a.Draw() + # fig.DrawNow() + + if output_dir is not None: + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + vv.screenshot(os.path.join(output_dir, '{}_3d_orientation.png'.format(name)), + sf=1, bg=background_color) + app.Run() + + +def plot_3d_diameter_map(name, data, unit_scale=1.0, measure_quantity='vox', + radius_structure_elem=1, output_dir=None, width=512, + height=512, camera_azth=44.5, camera_elev=35.8, + camera_roll=0.0, camera_fov=35.0, camera_zoom=0.0035, + camera_loc=(67.0, 81.6, 45.2), xlabel='', ylabel='', + zlabel='', axis_color='w', background_color='k', + cb_x_offset=10): + """Renders orientation data in 3D with RGB angular color-coding. + + Parameters + ---------- + name : str + Indicates the name of the output png file. + + data : 3D array + Indicates the 3D array containing diameter at every point of the skeleton. + + unit_scale : float + Indicates the scale factor of the data values. + + measure_quantity : str + Indicates the name of measure of the values. + + radius_structure_elem : integer + Indicates the size of the structure element of the dilation process to + thicken the skeleton. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + + camera_azth : float + Indicates the azimuth angle of the camera. + + width : int + Indicates the width of the visualization window. + + height : int + Indicates the width of the visualization window. + + camera_elev : float + Indicates the latitude / elevation angle of the camera. + + camera_roll : float + Indicates the roll angle of the camera. + + camera_fov : float + Indicates the field of view of the camera. + + camera_zoom : float + Indicates the zoom level of the camera. + + camera_loc : tuple + Indicates the camera location. + + xlabel : str + Indicates the label along the x-axis. + + ylabel : str + Indicates the label along the y-axis. + + zlabel : str + Indicates the label along the z-axis. + + axis_color : str + Indicates the color of axes. + + background_color : str + Indicates the background color of the figure. + + cb_x_offset : int + Indicates the offset of the colorbar from the right window side. + """ + if not visvis_available: + print('The visvis package is not found. The visualization cannot be done.') + return + + rmin, rmax, cmin, cmax, zmin, zmax = _bbox_3D(data) + dmtr = data[rmin:rmax, cmin:cmax, zmin:zmax] * unit_scale + skel = np.zeros_like(dmtr, dtype=np.uint8) + skel[dmtr.nonzero()] = 1 + + dmtr = ndi.grey_dilation(dmtr, structure=morphology.ball(radius_structure_elem)) + skel = ndi.binary_dilation(skel, structure=morphology.ball(radius_structure_elem)).astype(np.float32) + skel[skel.nonzero()] = 1. + + dmtr = dmtr * skel + + app = vv.use() + + fig = vv.figure() + fig._currentAxes = None + fig.relativeFontSize = 2. + fig.position.w = width + fig.position.h = height + + t = vv.volshow(dmtr[:, :, :], renderStyle='iso') + t.isoThreshold = 0.5 + t.colormap = vv.CM_JET + + a = vv.gca() + a.camera.azimuth = camera_azth + a.camera.elevation = camera_elev + a.camera.roll = camera_roll + a.camera.fov = camera_fov + a.camera.zoom = camera_zoom + a.camera.loc = camera_loc + + a.bgcolor = background_color + a.axis.axisColor = axis_color + a.axis.xLabel = xlabel + a.axis.yLabel = ylabel + a.axis.zLabel = zlabel + + cb = vv.colorbar() + cb.SetLabel('Diameter, [{}]'.format(measure_quantity)) + cb._label.position.x += cb_x_offset + + if output_dir is not None: + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + vv.screenshot(os.path.join(output_dir, '{}_3d_diameter.png'.format(name)), + sf=1, bg='w') + app.Run() + + +def plot_color_wheel(name, output_dir=None, dpi=500, xlabel='Elevation', + ylabel='Azimuth', fontsize=10, num_xticks=4, yticks=(-90, 90)): + """Plots the color wheel for visualizations of 3D orintation. + + Parameters + ---------- + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + + dpi : integer + Indicates the DPI of the output image. + + xlabel : str + Indicates the text along the x-axis. + + ylabel : str + Indicates the text along the y-axis. + + fontsize : int + Indicates the font size of labels along axes. + + num_xticks : int + Indicates the number of ticks along axes. + + yticks : tuple + Indicates the range of minimum and maximum values along the y-axis. + """ + + azth, lat = np.linspace(0., 1., num=180), np.linspace(0., 1., num=90) + rgb_arr = np.zeros((len(azth), len(lat), 3)) + for i in xrange(len(azth)): + for j in xrange(len(lat)): + rgb_arr[i, j, :] = colors.hsv_to_rgb([azth[i], lat[j], 1.0]) + + fig, ax = plt.subplots(figsize=(2, 2)) + ax.set_facecolor('red') + + ax.set_xlim([0, 90]) + ax.set_ylim([0, 181]) + + ax.set_yticks(np.linspace(0, 180, num=7).astype(np.int32)) + ax.set_yticklabels(np.linspace(yticks[0], yticks[1], num=7).astype(np.int32)) + + ax.set_xlabel(xlabel, fontsize=fontsize, labelpad=0, color='w') + ax.set_ylabel(ylabel, fontsize=fontsize, labelpad=0, color='w') + + xmajorlocator = ticker.LinearLocator(num_xticks) + ax.xaxis.set_major_locator(xmajorlocator) + ax.tick_params(direction='out', length=2, width=0, labelsize=fontsize, pad=0, colors='w') + + ax.imshow(rgb_arr) + + if output_dir is not None: + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + fig.savefig(os.path.join(output_dir, '{}_color_bar.png'.format(name)), + transparent=True, + bbox_inches='tight', + pad_inches=0.1, + dpi=dpi) + plt.show() + +def plot_fourier_orientation(data, orient_blocks, block_shape, figsize=(12,12), + cmap='gray', line_length=20, line_width=2.5, + line_color='red', line_style='-', name=None, + output_dir=None, dpi=200): + """Plots the orientation map in a block-wise manner. + + Plots the orientation vector over the image `data` at the center of each + block of the subdivided image. The orientation at every block 'orient_blocks' + and its size `block_shape` are specified by orientation estimation method. + The result can be and stored as a png file with DPI of `dpi` to the folder + specified by `output_dir`. + + Parameters + ---------- + data : ndarray + An image on top of which, the orientation will be plotted.. + + orient_blocks : ndarray + 2D array of orientation at every block of the subdivided image. + + block_shape : tuple + Indicates the block size within which the orientation is calculated. + + figsize : tuple of integers + Indicates the size of the output figure. + + cmap : str + Indicates the name of a colormap used for image. + + line_length : integer + Indicates the length of the orientation vector at each block. + + line_width : float + Indicates the line width of the orientation vector at each block. + + line_color : str + Indicates the line color of the orientation vector at each block. + + line_style : str + Indicates the line style of the orientation vector at each block. + + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + + dpi : integer + Indicates the DPI of the output image. + """ + fig, ax = plt.subplots(figsize=figsize) + + ax.xaxis.grid(False) + ax.yaxis.grid(False) + ax.set_axis_off() + + ax.imshow(data, cmap=cmap) + + for i in xrange(orient_blocks.shape[0]): + for j in xrange(orient_blocks.shape[1]): + y0, x0 = block_shape[0] * j + block_shape[0]/2, \ + block_shape[1] * i + block_shape[1]/2 + + orientation = orient_blocks[i ,j] + + x2 = x0 - math.sin(orientation) * line_length + y2 = y0 - math.cos(orientation) * line_length + + x3 = x0 + math.sin(orientation) * line_length + y3 = y0 + math.cos(orientation) * line_length + + ax.plot((x2, x3), (y2, y3), linestyle=line_style, + linewidth=line_width, color=line_color) + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + plt.tight_layout() + fig.savefig(os.path.join(output_dir, '{}_fourier_orientation_map.png'.format(name)), + bbox_inches='tight', transparent=True, pad_inches=0.05, dpi=dpi) + + plt.show() diff --git a/quanfima/quanfima/visualization_testO.py b/quanfima/quanfima/visualization_testO.py new file mode 100644 index 0000000..09af8a8 --- /dev/null +++ b/quanfima/quanfima/visualization_testO.py @@ -0,0 +1,1208 @@ +from __future__ import print_function +import os +import math +import numpy as np +import matplotlib as mpl +from skimage import io, filters +import matplotlib.colors as colors +import matplotlib.pyplot as plt +import matplotlib.ticker as ticker +import mpl_toolkits.axisartist.floating_axes as floating_axes +from matplotlib import cm +from scipy import ndimage as ndi, interpolate +from skimage import morphology +from matplotlib.projections import PolarAxes +from matplotlib.transforms import Affine2D +from mpl_toolkits.axisartist.grid_finder import DictFormatter, FixedLocator +from mpl_toolkits.axes_grid1 import make_axes_locatable +from mpl_toolkits.axisartist.floating_axes import GridHelperCurveLinear +from quanfima.quanfima.utils import geo2rgb +from quanfima.quanfima import visvis_available +if visvis_available: + import visvis as vv + +def plot_orientM_map(orient_map, fiber_skel, radius_structure_elem=1, + figsize=(15, 15), cmap='hsv', min_label=None, max_label=None, dpi=200, + labelsize=20, label='Orientation, [pixels]', name=None, + output_dir=None): + """Plots the diameter map with the colorbar. + + Plots the diameter map from the provided diameters `thickness_map` and + skeleton `fiber_skel` of size `figsize` using the colormap `cmap` and the + limits of the colorbar specified by `tmin` and `tmax`, and writes as a png + file with DPI of `dpi` to the folder specified by `output_dir`. + + Parameters + ---------- + thickness_map : ndarray + 2D array of diameter at every point of the skeleton. + + fiber_skel : ndarray + The binary skeleton extracted from the binary data. + + radius_structure_elem : integer + Indicates the size of the structure element of the dilation process to + thicken the skeleton. + + figsize : tuple of integers + Indicates the size of the output figure. + + cmap : str + Indicates the name of a colormap used to map angles to colors. + + tmin : float + Indicates the minimum value of the colorbar. + + tmax : float + Indicates the maximum value of the colorbar. + + dpi : integer + Indicates the DPI of the output image. + + labelsize : integer + Indicates the fontsize of the label of the colorbar. + + label : str + Indicates the label of the colorbar. + + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + """ + disk = morphology.disk(radius_structure_elem) + thickness_map = ndi.grey_dilation(orient_map, structure=disk).astype(np.float32) + fiber_skel = ndi.binary_dilation(fiber_skel, structure=disk).astype(np.float32) + + masked_thickness_map = np.ma.masked_where(fiber_skel == 0, thickness_map) + + cmap_obj = cm.get_cmap(cmap) + cmap_obj.set_bad(color='black') + + fig = plt.figure(figsize=figsize) + #ax = plt.subplot(111) + #ax.set_axis_off() + + #im = ax.imshow(masked_thickness_map, cmap=cmap_obj, vmin=tmin, vmax=tmax) + #divider = make_axes_locatable(ax) + #cax = divider.append_axes('right', size="2.5%", pad=0.05) + #cbar = plt.colorbar(im, cax=cax) + #cbar.ax.tick_params(labelsize=labelsize) + #cbar.set_label(label, fontsize=labelsize) + + + quant_steps = 2056 + cmap_obj = cm.get_cmap(cmap, quant_steps) + cmap_obj.set_bad(color='black') + + omap_ax = fig.add_axes([0.0, 0.0, 1.0, 1.0]) + omap_ax.set_axis_off() + omap_ax.imshow(masked_thickness_map, cmap=cmap_obj, interpolation=None) + + display_axes = fig.add_axes([0.780, -0.076, 0.2, 0.2], projection='polar') + display_axes_direction = np.pi + + norm = mpl.colors.Normalize(0.0, np.pi) + cb = mpl.colorbar.ColorbarBase(display_axes, cmap=cmap_obj, norm=norm, orientation='horizontal') + + display_axes.text(0.09, 0.56, min_label, color='white', fontsize=20, weight='bold', horizontalalignment='center', + verticalalignment='center', transform=display_axes.transAxes) + + display_axes.text(0.85, 0.56, max_label, color='white', fontsize=20, weight='bold', horizontalalignment='center', + verticalalignment='center', transform=display_axes.transAxes) + + cb.outline.set_visible(False) + display_axes.set_axis_off() + + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + fig.savefig(os.path.join(output_dir, '{}_diameter_map.png'.format(name)), + transparent=False, + bbox_inches='tight', + pad_inches=0.1, + dpi=dpi) + plt.show() + + + +def plot_orientation_test_map(orient_map, fiber_skel, radius_structure_elem=1, + figsize=(15, 15), cmap='hsv', dpi=200, min_label='0', + max_label='180', name=None, output_dir=None): + """Plots the orientation map with the color wheel. + + Plots the orientation map from the provided angles `orient_map` and skeleton + `fiber_skel` of size `figsize` using the colormap `cmap` and writes as a png + file with DPI of `dpi` to the folder specified by `output_dir`. + + Parameters + ---------- + orient_map : ndarray + 2D array of orientation at every point of the skeleton. + + fiber_skel : ndarray + The binary skeleton extracted from the binary data. + + radius_structure_elem : integer + Indicates the size of the structure element of the dilation process to + thicken the skeleton. + + figsize : tuple of integers + Indicates the size of the output figure. + + cmap : str + Indicates the name of a colormap used to map angles to colors. + + dpi : integer + Indicates the DPI of the output image. + + min_label : str + Indicates the label of minimum degree. + + max_label : str + Indicates the label of minimum degree. + + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + """ + disk = morphology.disk(radius_structure_elem) + orient_map = ndi.grey_dilation(orient_map, structure=disk).astype(np.float32) + fiber_skel = ndi.binary_dilation(fiber_skel, structure=disk).astype(np.float32) + + masked_orient_map = np.ma.masked_where(fiber_skel == 0, orient_map) + + fig = plt.figure(figsize=figsize) + ax = plt.subplot(111) + ax.set_axis_off() + + quant_steps = 2056 + + cmap_obj = cm.get_cmap(cmap, quant_steps) + cmap_obj.set_bad(color='black') + + im = ax.imshow(masked_thickness_map, cmap=cmap_obj, vmin=tmin, vmax=tmax) + divider = make_axes_locatable(ax) + cax = divider.append_axes('right', size="2.5%", pad=0.05) + cbar = plt.colorbar(im, cax=cax) + cbar.ax.tick_params(labelsize=labelsize) + cbar.set_label(label, fontsize=labelsize) + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + fig.savefig(os.path.join(output_dir, '{}_orient_map.png'.format(name)), + transparent=False, + bbox_inches='tight', + pad_inches=0.1, + dpi=dpi) + plt.show() + + #omap_ax = fig.add_axes([0.0, 0.0, 1.0, 1.0]) + #omap_ax.set_axis_off() + #omap_ax.imshow(masked_orient_map, cmap=cmap_obj, interpolation=None) + + #display_axes = fig.add_axes([0.780, -0.076, 0.2, 0.2], projection='polar') + #display_axes._direction = np.pi + + #norm = mpl.colors.Normalize(0.0, np.pi) + #cb = mpl.colorbar.ColorbarBase(display_axes, cmap=cmap_obj, norm=norm, orientation='horizontal') + + #display_axes.text(0.09, 0.56, min_label, color='white', fontsize=20, weight='bold', horizontalalignment='center', + # verticalalignment='center', transform=display_axes.transAxes) + + #display_axes.text(0.85, 0.56, max_label, color='white', fontsize=20, weight='bold', horizontalalignment='center', + # verticalalignment='center', transform=display_axes.transAxes) + + #cb.outline.set_visible(False) + #display_axes.set_axis_off() + + #if (output_dir is not None) and (name is not None): + # if not os.path.exists(output_dir): + # os.makedirs(output_dir) + # fig.savefig(os.path.join(output_dir, '{}_diameter_map.png'.format(name)), + # transparent=False, + # bbox_inches='tight', + # pad_inches=0.1, + # dpi=dpi) + + #plt.show() + + +def plot_diameter_map_test(thickness_map, fiber_skel, radius_structure_elem=1, + figsize=(15, 15), cmap='hsv', tmin=None, tmax=None, dpi=200, + labelsize=20, label='Diameter, [pixels]', name=None, + output_dir=None): + """Plots the diameter map with the colorbar. + + Plots the diameter map from the provided diameters `thickness_map` and + skeleton `fiber_skel` of size `figsize` using the colormap `cmap` and the + limits of the colorbar specified by `tmin` and `tmax`, and writes as a png + file with DPI of `dpi` to the folder specified by `output_dir`. + + Parameters + ---------- + thickness_map : ndarray + 2D array of diameter at every point of the skeleton. + + fiber_skel : ndarray + The binary skeleton extracted from the binary data. + + radius_structure_elem : integer + Indicates the size of the structure element of the dilation process to + thicken the skeleton. + + figsize : tuple of integers + Indicates the size of the output figure. + + cmap : str + Indicates the name of a colormap used to map angles to colors. + + tmin : float + Indicates the minimum value of the colorbar. + + tmax : float + Indicates the maximum value of the colorbar. + + dpi : integer + Indicates the DPI of the output image. + + labelsize : integer + Indicates the fontsize of the label of the colorbar. + + label : str + Indicates the label of the colorbar. + + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + """ + disk = morphology.disk(radius_structure_elem) + thickness_map = ndi.grey_dilation(thickness_map, structure=disk).astype(np.float32) + fiber_skel = ndi.binary_dilation(fiber_skel, structure=disk).astype(np.float32) + + masked_thickness_map = np.ma.masked_where(fiber_skel == 0, thickness_map) + + cmap_obj = cm.get_cmap(cmap) + cmap_obj.set_bad(color='black') + + fig = plt.figure(figsize=figsize) + ax = plt.subplot(111) + ax.set_axis_off() + + im = ax.imshow(masked_thickness_map, cmap=cmap_obj, vmin=tmin, vmax=tmax) + divider = make_axes_locatable(ax) + cax = divider.append_axes('right', size="2.5%", pad=0.05) + cbar = plt.colorbar(im, cax=cax) + cbar.ax.tick_params(labelsize=labelsize) + cbar.set_label(label, fontsize=labelsize) + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + fig.savefig(os.path.join(output_dir, '{}_diameter_map.png'.format(name)), + transparent=False, + bbox_inches='tight', + pad_inches=0.1, + dpi=dpi) + plt.show() + + +def gather_polar_errors(datasets_path, estimation_path, + azth_rng=np.arange(-90, 91, step=3), + lat_rng=np.arange(0, 91, step=3)): + """Computes the absolute angular error in ranges between estimated datasets. + + Calculates the absolute angular error in ranges `azth_rng` and `lat_rng` + between estimated orientation datasets localed at paths `datasets_path` and + `estimation_path`. + + Parameters + ---------- + datasets_path : str + Indicates the path to the reference / estimated dataset. + + estimation_path : str + Indicates the path to the estimated / reference dataset. + + azth_rng : array + Indicates the ranges of azimuth angles where the error is accumulated. + + lat_rng : array + Indicates the ranges of latitude or elevation angles where the error + is accumulated. + + Returns + ------- + array : 2D array + The 2D array of accumulated errors within the combinations of the specified + angular ranges. + """ + reference_dataset = np.load(datasets_path).item() + estimated_dataset = np.load(estimation_path).item() + idxs = estimated_dataset['indices'] + + print(reference_dataset.keys()) + print(estimated_dataset.keys()) + + ref_azth, ref_lat = reference_dataset['azth'][idxs], reference_dataset['lat'][idxs] + res_azth, res_lat = estimated_dataset['azth'][idxs], estimated_dataset['lat'][idxs] + + abs_err_azth, abs_err_lat = np.abs(ref_azth - res_azth), np.abs(ref_lat - res_lat) + + out = _angle_err2mean_abs_err(abs_err_azth, abs_err_lat, + ref_azth, ref_lat, + azth_rng=azth_rng, lat_rng=lat_rng) + + return out + + +def _angle_err2mean_abs_err(azth_err, lat_err, azth_ref, lat_ref, + azth_rng=np.arange(-90, 91, step=5), + lat_rng=np.arange(0, 91, step=5)): + """Computes the absolute angular error in ranges between estimated datasets. + + Parameters + ---------- + azth_err : 3D array + Indicates the error of azimuth components. + + lat_err : 3D array + Indicates the error of latitude / elevation components. + + azth_ref : 3D array + Indicates the reference dataset of estimated azimuth component. + + lat_ref : 3D array + Indicates the reference dataset of estimated latitude or elevation component. + + azth_rng : array + Indicates the ranges of azimuth angles where the error is accumulated. + + lat_rng : array + Indicates the ranges of latitude and elevation angles where the error + is accumulated. + + Returns + ------- + array : 2D array + The 2D array of accumulated errors within the combinations of the specified + angular ranges. + """ + out = np.zeros((len(azth_rng) - 1, len(lat_rng) - 1), dtype=np.float32) + + for i in xrange(len(azth_rng) - 1): + for j in xrange(len(lat_rng) - 1): + rng_azth = (np.deg2rad(azth_rng[i]), np.deg2rad(azth_rng[i+1])) + rng_lat = (np.deg2rad(lat_rng[j]), np.deg2rad(lat_rng[j+1])) + + idxs_azth = np.where((azth_ref >= rng_azth[0]) & (azth_ref < rng_azth[1])) + idxs_lat = np.where((lat_ref >= rng_lat[0]) & (lat_ref < rng_lat[1])) + + idxs = np.intersect1d(idxs_azth, idxs_lat) + + if len(idxs): + mu_azth, mu_lat = np.rad2deg(np.mean(azth_err[idxs])), \ + np.rad2deg(np.mean(lat_err[idxs])) + + out[i, j] = mu_azth + mu_lat + else: + out[i, j] = 0.0 + + return out + + +def plot_polar_heatmap(data, name, interp_factor=5., color_limits=False, + hide_colorbar=False, vmin=None, vmax=None, log_scale=True, + dpi=200, output_dir=None): + """Plots the polar heatmap describing azimuth and latitude / elevation components. + + Plots the polar heatmap where each cell of the heatmap corresponds to + the specific element of the array provided by `gather_polar_errors` + function. + + Parameters + ---------- + data : 2D array + Indicates the array containing the sum of angular errors within the + specified angular ranges. It is usually provided by `gather_polar_errors` + function. + + name : str + Indicates the name of the output png file. + + interp_factor : float + Indicates the interpolation factor of the heatmap. + + color_limits : boolean + Specifies if the determined intensity limits should be returned. + + hide_colorbar : boolean + Specifies if the colorbar should be hidden. + + vmin : float + Indicates the minimum value of the colorbar. + + vmax : float + Indicates the maximum value of the colorbar. + + log_scale : float + Specifies if the heatmap sould be in the logarithmic scale. + + dpi : integer + Indicates the DPI of the output image. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + """ + th0, th1 = 0., 180. + r0, r1 = 0, 90 + thlabel, rlabel = 'Azimuth', 'Elevation' + + tr_scale = Affine2D().scale(np.pi/180., 1.) + tr = tr_scale + PolarAxes.PolarTransform() + + lat_ticks = [(.0*90., '0$^{\circ}$'), + (.33*90., '30$^{\circ}$'), + (.66*90., '60$^{\circ}$'), + (1.*90., '90$^{\circ}$')] + r_grid_locator = FixedLocator([v for v, s in lat_ticks]) + r_grid_formatter = DictFormatter(dict(lat_ticks)) + + angle_ticks = [(0*180., '90$^{\circ}$'), + (.25*180., '45$^{\circ}$'), + (.5*180., '0$^{\circ}$'), + (.75*180., '-45$^{\circ}$'), + (1.*180., '-90$^{\circ}$')] + theta_grid_locator = FixedLocator([v for v, s in angle_ticks]) + theta_tick_formatter = DictFormatter(dict(angle_ticks)) + + grid_helper = GridHelperCurveLinear(tr, + extremes=(th0, th1, r0, r1), + grid_locator1=theta_grid_locator, + grid_locator2=r_grid_locator, + tick_formatter1=theta_tick_formatter, + tick_formatter2=r_grid_formatter) + + fig = plt.figure() + ax = floating_axes.FloatingSubplot(fig, 111, grid_helper=grid_helper) + fig.add_subplot(ax) + ax.set_facecolor('white') + + ax.axis["bottom"].set_visible(False) + ax.axis["top"].toggle(ticklabels=True, label=True) + ax.axis["top"].set_axis_direction("bottom") + ax.axis["top"].major_ticklabels.set_axis_direction("top") + ax.axis["top"].label.set_axis_direction("top") + + ax.axis["left"].set_axis_direction("bottom") + ax.axis["right"].set_axis_direction("top") + + ax.axis["top"].label.set_text(thlabel) + ax.axis["left"].label.set_text(rlabel) + + aux_ax = ax.get_aux_axes(tr) + aux_ax.patch = ax.patch + ax.patch.zorder = 0.9 + + rad = np.linspace(0, 90, data.shape[1]) + azm = np.linspace(0, 180, data.shape[0]) + + f = interpolate.interp2d(rad, azm, data, kind='linear', bounds_error=True, fill_value=0) + + new_rad = np.linspace(0, 90, 180*interp_factor) + new_azm = np.linspace(0, 180, 360*interp_factor) + new_data_angle_dist = f(new_rad, new_azm) + new_r, new_th = np.meshgrid(new_rad, new_azm) + new_data_angle_dist += 1. + + if log_scale: + data_mesh = aux_ax.pcolormesh(new_th, new_r, new_data_angle_dist, cmap='jet', + norm=colors.LogNorm(vmin=1. if vmin is None else vmin, + vmax=new_data_angle_dist.max() if vmax is None else vmax)) + else: + data_mesh = aux_ax.pcolormesh(new_th, new_r, new_data_angle_dist, cmap='jet', vmin=vmin, vmax=vmax) + + cbar = plt.colorbar(data_mesh, orientation='vertical', shrink=.88, pad=.1, aspect=15) + cbar.ax.set_ylabel('Absolute error, [deg.]') + + if hide_colorbar: + cbar.remove() + + ax.grid(False) + + plt.show() + + if output_dir is not None: + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + fig.savefig(os.path.join(output_dir, '{}_chart.png'.format(name)), + transparent=False, + bbox_inches='tight', + pad_inches=0.1, + dpi=dpi) + + if color_limits: + return 1., new_data_angle_dist.max() + + +def plot_histogram_fig(data, num_bins, xticks, color, splot_index=111, output_dir=None, + xlim=(None, None), ylim=(None, None), name=None, in_percent=False, + bar_width=0.8, ticks_pad=7, xticks_fontsize=22, + yticks_fontsize=22, xlabel=None, ylabel=None, + labels_fontsize=20, grid_alpha=0.3, title_fontsize=22, + exp_fontsize=15, type=None, figsize=(12, 8), dpi=200): + """Plots the histogram from a given data. + + Parameters + ---------- + data : 1D array + Indicates the array containing the values. + + num_bins : integer + Indicates the number of histogram bins. + + xticks : array + Indicates the array of ticks of the X-axis. + + color : str + Indicates the color of the histogram. + + Returns + ------- + tuple : tuple of Figure and Axis objects. + """ + fig = plt.figure(figsize=figsize) + ax = fig.add_subplot(splot_index) + + weights = np.ones_like(data)/float(len(data)) + + n, bins, patches = ax.hist(data, num_bins, color=color, + rwidth=bar_width, weights=weights) + + ax.tick_params(axis='x', + labelsize=xticks_fontsize, + colors='#000000', + which='both', + direction='out', + length=6, + width=2, + pad=ticks_pad) + ax.tick_params(axis='x', + labelsize=xticks_fontsize, + colors='#000000', + which='minor', + direction='out', + length=4, + width=2, + pad=ticks_pad) + + ax.tick_params(axis='y', + labelsize=yticks_fontsize, + colors='#000000', + which='major', + direction='out', + length=6, + width=2, + pad=ticks_pad) + ax.tick_params(axis='y', + labelsize=xticks_fontsize, + colors='#000000', + which='minor', + direction='out', + length=4, + width=2, + pad=ticks_pad) + + ax.xaxis.set_minor_locator(ticker.MultipleLocator(5)) + ax.yaxis.set_minor_locator(ticker.MultipleLocator(5)) + + plt.xticks(xticks, xticks) + + ax.xaxis.offsetText.set_fontsize(exp_fontsize) + ax.yaxis.offsetText.set_fontsize(exp_fontsize) + + ax.ticklabel_format(axis='y', style='sci', scilimits=(-2, 2), + fontsize=yticks_fontsize) + ax.yaxis.set_major_locator(ticker.MaxNLocator(nbins=6, min_n_ticks=6)) + + ax.xaxis.grid(False) + ax.yaxis.grid(False) + + ax.set_ylabel(ylabel, labelpad=2, fontsize=labels_fontsize, color='black') + ax.set_xlabel(xlabel, labelpad=2, fontsize=labels_fontsize, color='black') + + ax.set_ylim(ylim) + ax.set_xlim(xlim) + + if in_percent: + vals = ax.get_yticks() + ax.set_yticklabels(['{:3.0f}'.format(x*100) for x in vals]) + + ax.set_axisbelow(True) + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + plt.tight_layout() + fig.savefig(os.path.join(output_dir, '{}_chart.png'.format(name)), + transparent=False, bbox_inches='tight', pad_inches=0.1, dpi=dpi) + + return fig, ax + + +def create_pie_chart(data, rngs, colors=['#244268', '#426084', '#67809F', '#95A9C1', '#C6D2E0'], + unit_scale=1.0, measure_quantity='m^3', figsize=(33, 15), + legend_loc=(0.383, -0.25), zebra_color=(False, 3), + legend_fontsize=50, chart_fontsize=60, dpi=72, name=None, + output_dir=None): + """Plots the piechart of from a given data. + + Parameters + ---------- + data : 1D array + Indicates the array containing the values. + + rngs : tuple of tuples + Indicates the ranges of the piechart. + + colors : array + Indicates the color for the region of the piechart corresponding to the + specific range. + + unit_scale : float + Indicates the scale factor of the data values. + + measure_quantity : str + Indicates the name of measure of the values. + + figsize : tuple of integers + Indicates the size of the output figure. + + legend_loc : tuple + Indicates the position of the legend of the figure. + + zebra_color : tuple + Allows to change the text color of the region to white from the first to + the speficied index of the region (True, reg_index). + + legend_fontsize : integer + Indicates the fontsize of the legend. + + chart_fontsize : integer + Indicates the fontsize of the figure. + + dpi : integer + Indicates the DPI of the output image. + + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + """ + def number(val): + if val < 1000: + return '%d' % val + + sv = str(val) + return '$\mathregular{10^{%d}}$' % (len(sv)-2) if val % 10 == 0 else '%0.0e' % val + + def get_title(v1, v2, measure_quantity): + ftm = '%s $\minus$ %s %s' + return ftm % (number(v1), number(v2), measure_quantity) + + data_ranges = [] + df = data * unit_scale + for rng in rngs: + rng_min, rng_max = rng[0], rng[1] + data_rng = df[(df > rng_min) & (df < rng_max)] + data_ranges.append(data_rng) + + num_elem = [len(p) for p in data_ranges] + se = sum(num_elem) + + print('Num of particles: {}'%(se)) + + proc_particles = [n/float(se) * 100.0 for n in num_elem] + + for size, rng in zip(num_elem, rngs): + print('{}-{}: {}'%(rng[0], rng[1], size)) + + titles = [get_title(minv, maxv, measure_quantity) for minv, maxv in rngs] + + textprops = {'fontsize': chart_fontsize, + 'weight': 'normal', + 'family': 'sans-serif'} + pie_width = 0.5 + fig, ax = plt.subplots(figsize=figsize) + ax.axis('equal') + patches, texts, autotexts = ax.pie(proc_particles, + textprops=textprops, + colors=colors, + autopct='%1.1f%%', + radius=1, + pctdistance=1-pie_width/2) + + if (zebra_color is not None) and (zebra_color[0]): + for tt in autotexts[:zebra_color[1]]: + tt.set_color('white') + + plt.setp(patches, + width=pie_width, + edgecolor='white') + + plt.legend(patches, titles, loc=legend_loc, fontsize=legend_fontsize) + + _d, _offset, _di = [1, -1], [0.45, 0.45], 0 + + for t, p in zip(autotexts, proc_particles): + if p < 2.0: + pos = list(t.get_position()) + pos[0] = pos[0] + _d[_di] * _offset[_di] + + t.set_position(pos) + _di += 1 + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + plt.tight_layout() + fig.savefig(os.path.join(output_dir, '{}_chart.png'.format(name)), + bbox_inches='tight', transparent=True, pad_inches=0.1, dpi=dpi) + + +def _bbox_3D(img): + """Crops the non-zero part of a volume + """ + r = np.any(img, axis=(1, 2)) + c = np.any(img, axis=(0, 2)) + z = np.any(img, axis=(0, 1)) + + rmin, rmax = np.where(r)[0][[0, -1]] + cmin, cmax = np.where(c)[0][[0, -1]] + zmin, zmax = np.where(z)[0][[0, -1]] + + return rmin, rmax, cmin, cmax, zmin, zmax + + +def plot_3d_orientation_map(name, lat_data, azth_data, radius_structure_elem=1, + output_dir=None, width=512, height=512, + camera_azth=44.5, camera_elev=35.8, camera_roll=0.0, + camera_fov=35.0, camera_zoom=0.0035, + camera_loc=(67.0, 81.6, 45.2), xlabel='', ylabel='', + zlabel='', axis_color='w', background_color='k'): + """Renders orientation data in 3D with RGB angular color-coding. + + Parameters + ---------- + name : str + Indicates the name of the output png file. + + lat_data : 3D array + Indicates the 3D array containing latitude / elevation angle at every point of + the skeleton in radians. + + azth_data : 3D array + Indicates the 3D array containing azimuth angle at every point of the skeleton + in radians. + + radius_structure_elem : integer + Indicates the size of the structure element of the dilation process to + thicken the skeleton. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + + width : int + Indicates the width of the visualization window. + + height : int + Indicates the width of the visualization window. + + camera_azth : float + Indicates the azimuth angle of the camera. + + camera_elev : float + Indicates the latitude / elevation angle of the camera. + + camera_roll : float + Indicates the roll angle of the camera. + + camera_fov : float + Indicates the field of view of the camera. + + camera_zoom : float + Indicates the zoom level of the camera. + + camera_loc : tuple + Indicates the camera location. + + xlabel : str + Indicates the label along the x-axis. + + ylabel : str + Indicates the label along the y-axis. + + zlabel : str + Indicates the label along the z-axis. + + axis_color : str + Indicates the color of axes. + + background_color : str + Indicates the background color of the figure. + """ + if not visvis_available: + print('The visvis package is not found. The visualization cannot be done.') + return + + rmin, rmax, cmin, cmax, zmin, zmax = _bbox_3D(azth_data) + + azth, lat = azth_data[rmin:rmax, cmin:cmax, zmin:zmax], \ + np.abs(lat_data[rmin:rmax, cmin:cmax, zmin:zmax]) + + skel = azth.copy().astype(np.float32) + skel[skel.nonzero()] = 1. + + azth = ndi.grey_dilation(azth, structure=morphology.ball(radius_structure_elem)) + lat = ndi.grey_dilation(lat, structure=morphology.ball(radius_structure_elem)) + skel = ndi.binary_dilation(skel, structure=morphology.ball(radius_structure_elem)) + + Z, Y, X = skel.nonzero() + vol_orient = np.zeros(skel.shape + (3,), dtype=np.float32) + + print(vol_orient.size, vol_orient[skel.nonzero()].size) + + for z, y, x in zip(Z, Y, X): + vol_orient[z, y, x] = geo2rgb(lat[z, y, x], azth[z, y, x]) + + app = vv.use() + + fig = vv.figure() + fig._currentAxes = None + fig.relativeFontSize = 2. + fig.position.w = width + fig.position.h = height + + t = vv.volshow(vol_orient[:, :, :], renderStyle='iso') + t.isoThreshold = 0.5 + + a = vv.gca() + a.camera.azimuth = camera_azth + a.camera.elevation = camera_elev + a.camera.roll = camera_roll + a.camera.fov = camera_fov + a.camera.zoom = camera_zoom + a.camera.loc = camera_loc + + a.bgcolor = background_color + a.axis.axisColor = axis_color + a.axis.xLabel = xlabel + a.axis.yLabel = ylabel + a.axis.zLabel = zlabel + + # def mouseUp(event): + # print 'mouseUp!!' + # a = vv.gca() + # print a.camera.GetViewParams() + # + # a.eventMouseUp.Bind(mouseUp) + # fig.eventMouseUp.Bind(mouseUp) + # + # a.Draw() + # fig.DrawNow() + + if output_dir is not None: + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + vv.screenshot(os.path.join(output_dir, '{}_3d_orientation.png'.format(name)), + sf=1, bg=background_color) + app.Run() + + +def plot_3d_diameter_map(name, data, unit_scale=1.0, measure_quantity='vox', + radius_structure_elem=1, output_dir=None, width=512, + height=512, camera_azth=44.5, camera_elev=35.8, + camera_roll=0.0, camera_fov=35.0, camera_zoom=0.0035, + camera_loc=(67.0, 81.6, 45.2), xlabel='', ylabel='', + zlabel='', axis_color='w', background_color='k', + cb_x_offset=10): + """Renders orientation data in 3D with RGB angular color-coding. + + Parameters + ---------- + name : str + Indicates the name of the output png file. + + data : 3D array + Indicates the 3D array containing diameter at every point of the skeleton. + + unit_scale : float + Indicates the scale factor of the data values. + + measure_quantity : str + Indicates the name of measure of the values. + + radius_structure_elem : integer + Indicates the size of the structure element of the dilation process to + thicken the skeleton. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + + camera_azth : float + Indicates the azimuth angle of the camera. + + width : int + Indicates the width of the visualization window. + + height : int + Indicates the width of the visualization window. + + camera_elev : float + Indicates the latitude / elevation angle of the camera. + + camera_roll : float + Indicates the roll angle of the camera. + + camera_fov : float + Indicates the field of view of the camera. + + camera_zoom : float + Indicates the zoom level of the camera. + + camera_loc : tuple + Indicates the camera location. + + xlabel : str + Indicates the label along the x-axis. + + ylabel : str + Indicates the label along the y-axis. + + zlabel : str + Indicates the label along the z-axis. + + axis_color : str + Indicates the color of axes. + + background_color : str + Indicates the background color of the figure. + + cb_x_offset : int + Indicates the offset of the colorbar from the right window side. + """ + if not visvis_available: + print('The visvis package is not found. The visualization cannot be done.') + return + + rmin, rmax, cmin, cmax, zmin, zmax = _bbox_3D(data) + dmtr = data[rmin:rmax, cmin:cmax, zmin:zmax] * unit_scale + skel = np.zeros_like(dmtr, dtype=np.uint8) + skel[dmtr.nonzero()] = 1 + + dmtr = ndi.grey_dilation(dmtr, structure=morphology.ball(radius_structure_elem)) + skel = ndi.binary_dilation(skel, structure=morphology.ball(radius_structure_elem)).astype(np.float32) + skel[skel.nonzero()] = 1. + + dmtr = dmtr * skel + + app = vv.use() + + fig = vv.figure() + fig._currentAxes = None + fig.relativeFontSize = 2. + fig.position.w = width + fig.position.h = height + + t = vv.volshow(dmtr[:, :, :], renderStyle='iso') + t.isoThreshold = 0.5 + t.colormap = vv.CM_JET + + a = vv.gca() + a.camera.azimuth = camera_azth + a.camera.elevation = camera_elev + a.camera.roll = camera_roll + a.camera.fov = camera_fov + a.camera.zoom = camera_zoom + a.camera.loc = camera_loc + + a.bgcolor = background_color + a.axis.axisColor = axis_color + a.axis.xLabel = xlabel + a.axis.yLabel = ylabel + a.axis.zLabel = zlabel + + cb = vv.colorbar() + cb.SetLabel('Diameter, [{}]'.format(measure_quantity)) + cb._label.position.x += cb_x_offset + + if output_dir is not None: + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + vv.screenshot(os.path.join(output_dir, '{}_3d_diameter.png'.format(name)), + sf=1, bg='w') + app.Run() + + +def plot_color_wheel(name, output_dir=None, dpi=500, xlabel='Elevation', + ylabel='Azimuth', fontsize=10, num_xticks=4, yticks=(-90, 90)): + """Plots the color wheel for visualizations of 3D orintation. + + Parameters + ---------- + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + + dpi : integer + Indicates the DPI of the output image. + + xlabel : str + Indicates the text along the x-axis. + + ylabel : str + Indicates the text along the y-axis. + + fontsize : int + Indicates the font size of labels along axes. + + num_xticks : int + Indicates the number of ticks along axes. + + yticks : tuple + Indicates the range of minimum and maximum values along the y-axis. + """ + + azth, lat = np.linspace(0., 1., num=180), np.linspace(0., 1., num=90) + rgb_arr = np.zeros((len(azth), len(lat), 3)) + for i in xrange(len(azth)): + for j in xrange(len(lat)): + rgb_arr[i, j, :] = colors.hsv_to_rgb([azth[i], lat[j], 1.0]) + + fig, ax = plt.subplots(figsize=(2, 2)) + ax.set_facecolor('red') + + ax.set_xlim([0, 90]) + ax.set_ylim([0, 181]) + + ax.set_yticks(np.linspace(0, 180, num=7).astype(np.int32)) + ax.set_yticklabels(np.linspace(yticks[0], yticks[1], num=7).astype(np.int32)) + + ax.set_xlabel(xlabel, fontsize=fontsize, labelpad=0, color='w') + ax.set_ylabel(ylabel, fontsize=fontsize, labelpad=0, color='w') + + xmajorlocator = ticker.LinearLocator(num_xticks) + ax.xaxis.set_major_locator(xmajorlocator) + ax.tick_params(direction='out', length=2, width=0, labelsize=fontsize, pad=0, colors='w') + + ax.imshow(rgb_arr) + + if output_dir is not None: + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + fig.savefig(os.path.join(output_dir, '{}_color_bar.png'.format(name)), + transparent=True, + bbox_inches='tight', + pad_inches=0.1, + dpi=dpi) + plt.show() + +def plot_fourier_orientation(data, orient_blocks, block_shape, figsize=(12,12), + cmap='gray', line_length=20, line_width=2.5, + line_color='red', line_style='-', name=None, + output_dir=None, dpi=200): + """Plots the orientation map in a block-wise manner. + + Plots the orientation vector over the image `data` at the center of each + block of the subdivided image. The orientation at every block 'orient_blocks' + and its size `block_shape` are specified by orientation estimation method. + The result can be and stored as a png file with DPI of `dpi` to the folder + specified by `output_dir`. + + Parameters + ---------- + data : ndarray + An image on top of which, the orientation will be plotted.. + + orient_blocks : ndarray + 2D array of orientation at every block of the subdivided image. + + block_shape : tuple + Indicates the block size within which the orientation is calculated. + + figsize : tuple of integers + Indicates the size of the output figure. + + cmap : str + Indicates the name of a colormap used for image. + + line_length : integer + Indicates the length of the orientation vector at each block. + + line_width : float + Indicates the line width of the orientation vector at each block. + + line_color : str + Indicates the line color of the orientation vector at each block. + + line_style : str + Indicates the line style of the orientation vector at each block. + + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + + dpi : integer + Indicates the DPI of the output image. + """ + fig, ax = plt.subplots(figsize=figsize) + + ax.xaxis.grid(False) + ax.yaxis.grid(False) + ax.set_axis_off() + + ax.imshow(data, cmap=cmap) + + for i in xrange(orient_blocks.shape[0]): + for j in xrange(orient_blocks.shape[1]): + y0, x0 = block_shape[0] * j + block_shape[0]/2, \ + block_shape[1] * i + block_shape[1]/2 + + orientation = orient_blocks[i ,j] + + x2 = x0 - math.sin(orientation) * line_length + y2 = y0 - math.cos(orientation) * line_length + + x3 = x0 + math.sin(orientation) * line_length + y3 = y0 + math.cos(orientation) * line_length + + ax.plot((x2, x3), (y2, y3), linestyle=line_style, + linewidth=line_width, color=line_color) + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + plt.tight_layout() + fig.savefig(os.path.join(output_dir, '{}_fourier_orientation_map.png'.format(name)), + bbox_inches='tight', transparent=True, pad_inches=0.05, dpi=dpi) + + plt.show() diff --git a/quanfima/quanfima/visualization_testQ.py b/quanfima/quanfima/visualization_testQ.py new file mode 100644 index 0000000..9e44485 --- /dev/null +++ b/quanfima/quanfima/visualization_testQ.py @@ -0,0 +1,1208 @@ +from __future__ import print_function +import os +import math +import numpy as np +import matplotlib as mpl +from skimage import io, filters +import matplotlib.colors as colors +import matplotlib.pyplot as plt +import matplotlib.ticker as ticker +import mpl_toolkits.axisartist.floating_axes as floating_axes +from matplotlib import cm +from scipy import ndimage as ndi, interpolate +from skimage import morphology +from matplotlib.projections import PolarAxes +from matplotlib.transforms import Affine2D +from mpl_toolkits.axisartist.grid_finder import DictFormatter, FixedLocator +from mpl_toolkits.axes_grid1 import make_axes_locatable +from mpl_toolkits.axisartist.floating_axes import GridHelperCurveLinear +from quanfima.quanfima.utils2 import geo2rgb +from quanfima.quanfima import visvis_available +if visvis_available: + import visvis as vv + +def plot_orientM_map(orient_map, fiber_skel, radius_structure_elem=1, + figsize=(15, 15), cmap='hsv', min_label=None, max_label=None, dpi=200, + labelsize=20, label='Orientation, [pixels]', name=None, + output_dir=None): + """Plots the diameter map with the colorbar. + + Plots the diameter map from the provided diameters `thickness_map` and + skeleton `fiber_skel` of size `figsize` using the colormap `cmap` and the + limits of the colorbar specified by `tmin` and `tmax`, and writes as a png + file with DPI of `dpi` to the folder specified by `output_dir`. + + Parameters + ---------- + thickness_map : ndarray + 2D array of diameter at every point of the skeleton. + + fiber_skel : ndarray + The binary skeleton extracted from the binary data. + + radius_structure_elem : integer + Indicates the size of the structure element of the dilation process to + thicken the skeleton. + + figsize : tuple of integers + Indicates the size of the output figure. + + cmap : str + Indicates the name of a colormap used to map angles to colors. + + tmin : float + Indicates the minimum value of the colorbar. + + tmax : float + Indicates the maximum value of the colorbar. + + dpi : integer + Indicates the DPI of the output image. + + labelsize : integer + Indicates the fontsize of the label of the colorbar. + + label : str + Indicates the label of the colorbar. + + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + """ + disk = morphology.disk(radius_structure_elem) + thickness_map = ndi.grey_dilation(orient_map, structure=disk).astype(np.float32) + fiber_skel = ndi.binary_dilation(fiber_skel, structure=disk).astype(np.float32) + + masked_thickness_map = np.ma.masked_where(fiber_skel == 0, thickness_map) + + cmap_obj = cm.get_cmap(cmap) + cmap_obj.set_bad(color='black') + + fig = plt.figure(figsize=figsize) + #ax = plt.subplot(111) + #ax.set_axis_off() + + #im = ax.imshow(masked_thickness_map, cmap=cmap_obj, vmin=tmin, vmax=tmax) + #divider = make_axes_locatable(ax) + #cax = divider.append_axes('right', size="2.5%", pad=0.05) + #cbar = plt.colorbar(im, cax=cax) + #cbar.ax.tick_params(labelsize=labelsize) + #cbar.set_label(label, fontsize=labelsize) + + + quant_steps = 2056 + cmap_obj = cm.get_cmap(cmap, quant_steps) + cmap_obj.set_bad(color='black') + + omap_ax = fig.add_axes([0.0, 0.0, 1.0, 1.0]) + omap_ax.set_axis_off() + omap_ax.imshow(orient_map, cmap=cmap_obj, interpolation=None) + + display_axes = fig.add_axes([0.780, -0.076, 0.2, 0.2], projection='polar') + display_axes_direction = np.pi + + norm = mpl.colors.Normalize(0.0, np.pi) + cb = mpl.colorbar.ColorbarBase(display_axes, cmap=cmap_obj, norm=norm, orientation='horizontal') + + display_axes.text(0.09, 0.56, min_label, color='white', fontsize=20, weight='bold', horizontalalignment='center', + verticalalignment='center', transform=display_axes.transAxes) + + display_axes.text(0.85, 0.56, max_label, color='white', fontsize=20, weight='bold', horizontalalignment='center', + verticalalignment='center', transform=display_axes.transAxes) + + cb.outline.set_visible(False) + display_axes.set_axis_off() + + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + fig.savefig(os.path.join(output_dir, '{}_diameter_map.png'.format(name)), + transparent=False, + bbox_inches='tight', + pad_inches=0.1, + dpi=dpi) + plt.show() + + + +def plot_orientation_test_map(orient_map, fiber_skel, radius_structure_elem=1, + figsize=(15, 15), cmap='hsv', dpi=200, min_label='0', + max_label='180', name=None, output_dir=None): + """Plots the orientation map with the color wheel. + + Plots the orientation map from the provided angles `orient_map` and skeleton + `fiber_skel` of size `figsize` using the colormap `cmap` and writes as a png + file with DPI of `dpi` to the folder specified by `output_dir`. + + Parameters + ---------- + orient_map : ndarray + 2D array of orientation at every point of the skeleton. + + fiber_skel : ndarray + The binary skeleton extracted from the binary data. + + radius_structure_elem : integer + Indicates the size of the structure element of the dilation process to + thicken the skeleton. + + figsize : tuple of integers + Indicates the size of the output figure. + + cmap : str + Indicates the name of a colormap used to map angles to colors. + + dpi : integer + Indicates the DPI of the output image. + + min_label : str + Indicates the label of minimum degree. + + max_label : str + Indicates the label of minimum degree. + + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + """ + disk = morphology.disk(radius_structure_elem) + orient_map = ndi.grey_dilation(orient_map, structure=disk).astype(np.float32) + fiber_skel = ndi.binary_dilation(fiber_skel, structure=disk).astype(np.float32) + + masked_orient_map = np.ma.masked_where(fiber_skel == 0, orient_map) + + fig = plt.figure(figsize=figsize) + ax = plt.subplot(111) + ax.set_axis_off() + + quant_steps = 2056 + + cmap_obj = cm.get_cmap(cmap, quant_steps) + cmap_obj.set_bad(color='black') + + im = ax.imshow(orient_map, cmap=cmap_obj, vmin=tmin, vmax=tmax) + divider = make_axes_locatable(ax) + cax = divider.append_axes('right', size="2.5%", pad=0.05) + cbar = plt.colorbar(im, cax=cax) + cbar.ax.tick_params(labelsize=labelsize) + cbar.set_label(label, fontsize=labelsize) + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + fig.savefig(os.path.join(output_dir, '{}_orient_map.png'.format(name)), + transparent=False, + bbox_inches='tight', + pad_inches=0.1, + dpi=dpi) + plt.show() + + #omap_ax = fig.add_axes([0.0, 0.0, 1.0, 1.0]) + #omap_ax.set_axis_off() + #omap_ax.imshow(masked_orient_map, cmap=cmap_obj, interpolation=None) + + #display_axes = fig.add_axes([0.780, -0.076, 0.2, 0.2], projection='polar') + #display_axes._direction = np.pi + + #norm = mpl.colors.Normalize(0.0, np.pi) + #cb = mpl.colorbar.ColorbarBase(display_axes, cmap=cmap_obj, norm=norm, orientation='horizontal') + + #display_axes.text(0.09, 0.56, min_label, color='white', fontsize=20, weight='bold', horizontalalignment='center', + # verticalalignment='center', transform=display_axes.transAxes) + + #display_axes.text(0.85, 0.56, max_label, color='white', fontsize=20, weight='bold', horizontalalignment='center', + # verticalalignment='center', transform=display_axes.transAxes) + + #cb.outline.set_visible(False) + #display_axes.set_axis_off() + + #if (output_dir is not None) and (name is not None): + # if not os.path.exists(output_dir): + # os.makedirs(output_dir) + # fig.savefig(os.path.join(output_dir, '{}_diameter_map.png'.format(name)), + # transparent=False, + # bbox_inches='tight', + # pad_inches=0.1, + # dpi=dpi) + + #plt.show() + + +def plot_diameter_map_test(thickness_map, fiber_skel, radius_structure_elem=1, + figsize=(15, 15), cmap='hsv', tmin=None, tmax=None, dpi=200, + labelsize=20, label='Diameter, [pixels]', name=None, + output_dir=None): + """Plots the diameter map with the colorbar. + + Plots the diameter map from the provided diameters `thickness_map` and + skeleton `fiber_skel` of size `figsize` using the colormap `cmap` and the + limits of the colorbar specified by `tmin` and `tmax`, and writes as a png + file with DPI of `dpi` to the folder specified by `output_dir`. + + Parameters + ---------- + thickness_map : ndarray + 2D array of diameter at every point of the skeleton. + + fiber_skel : ndarray + The binary skeleton extracted from the binary data. + + radius_structure_elem : integer + Indicates the size of the structure element of the dilation process to + thicken the skeleton. + + figsize : tuple of integers + Indicates the size of the output figure. + + cmap : str + Indicates the name of a colormap used to map angles to colors. + + tmin : float + Indicates the minimum value of the colorbar. + + tmax : float + Indicates the maximum value of the colorbar. + + dpi : integer + Indicates the DPI of the output image. + + labelsize : integer + Indicates the fontsize of the label of the colorbar. + + label : str + Indicates the label of the colorbar. + + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + """ + disk = morphology.disk(radius_structure_elem) + thickness_map = ndi.grey_dilation(thickness_map, structure=disk).astype(np.float32) + fiber_skel = ndi.binary_dilation(fiber_skel, structure=disk).astype(np.float32) + + masked_thickness_map = np.ma.masked_where(fiber_skel == 0, thickness_map) + + cmap_obj = cm.get_cmap(cmap) + cmap_obj.set_bad(color='black') + + fig = plt.figure(figsize=figsize) + ax = plt.subplot(111) + ax.set_axis_off() + + im = ax.imshow(masked_thickness_map, cmap=cmap_obj, vmin=tmin, vmax=tmax) + divider = make_axes_locatable(ax) + cax = divider.append_axes('right', size="2.5%", pad=0.05) + cbar = plt.colorbar(im, cax=cax) + cbar.ax.tick_params(labelsize=labelsize) + cbar.set_label(label, fontsize=labelsize) + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + fig.savefig(os.path.join(output_dir, '{}_diameter_map.png'.format(name)), + transparent=False, + bbox_inches='tight', + pad_inches=0.1, + dpi=dpi) + plt.show() + + +def gather_polar_errors(datasets_path, estimation_path, + azth_rng=np.arange(-90, 91, step=3), + lat_rng=np.arange(0, 91, step=3)): + """Computes the absolute angular error in ranges between estimated datasets. + + Calculates the absolute angular error in ranges `azth_rng` and `lat_rng` + between estimated orientation datasets localed at paths `datasets_path` and + `estimation_path`. + + Parameters + ---------- + datasets_path : str + Indicates the path to the reference / estimated dataset. + + estimation_path : str + Indicates the path to the estimated / reference dataset. + + azth_rng : array + Indicates the ranges of azimuth angles where the error is accumulated. + + lat_rng : array + Indicates the ranges of latitude or elevation angles where the error + is accumulated. + + Returns + ------- + array : 2D array + The 2D array of accumulated errors within the combinations of the specified + angular ranges. + """ + reference_dataset = np.load(datasets_path).item() + estimated_dataset = np.load(estimation_path).item() + idxs = estimated_dataset['indices'] + + print(reference_dataset.keys()) + print(estimated_dataset.keys()) + + ref_azth, ref_lat = reference_dataset['azth'][idxs], reference_dataset['lat'][idxs] + res_azth, res_lat = estimated_dataset['azth'][idxs], estimated_dataset['lat'][idxs] + + abs_err_azth, abs_err_lat = np.abs(ref_azth - res_azth), np.abs(ref_lat - res_lat) + + out = _angle_err2mean_abs_err(abs_err_azth, abs_err_lat, + ref_azth, ref_lat, + azth_rng=azth_rng, lat_rng=lat_rng) + + return out + + +def _angle_err2mean_abs_err(azth_err, lat_err, azth_ref, lat_ref, + azth_rng=np.arange(-90, 91, step=5), + lat_rng=np.arange(0, 91, step=5)): + """Computes the absolute angular error in ranges between estimated datasets. + + Parameters + ---------- + azth_err : 3D array + Indicates the error of azimuth components. + + lat_err : 3D array + Indicates the error of latitude / elevation components. + + azth_ref : 3D array + Indicates the reference dataset of estimated azimuth component. + + lat_ref : 3D array + Indicates the reference dataset of estimated latitude or elevation component. + + azth_rng : array + Indicates the ranges of azimuth angles where the error is accumulated. + + lat_rng : array + Indicates the ranges of latitude and elevation angles where the error + is accumulated. + + Returns + ------- + array : 2D array + The 2D array of accumulated errors within the combinations of the specified + angular ranges. + """ + out = np.zeros((len(azth_rng) - 1, len(lat_rng) - 1), dtype=np.float32) + + for i in xrange(len(azth_rng) - 1): + for j in xrange(len(lat_rng) - 1): + rng_azth = (np.deg2rad(azth_rng[i]), np.deg2rad(azth_rng[i+1])) + rng_lat = (np.deg2rad(lat_rng[j]), np.deg2rad(lat_rng[j+1])) + + idxs_azth = np.where((azth_ref >= rng_azth[0]) & (azth_ref < rng_azth[1])) + idxs_lat = np.where((lat_ref >= rng_lat[0]) & (lat_ref < rng_lat[1])) + + idxs = np.intersect1d(idxs_azth, idxs_lat) + + if len(idxs): + mu_azth, mu_lat = np.rad2deg(np.mean(azth_err[idxs])), \ + np.rad2deg(np.mean(lat_err[idxs])) + + out[i, j] = mu_azth + mu_lat + else: + out[i, j] = 0.0 + + return out + + +def plot_polar_heatmap(data, name, interp_factor=5., color_limits=False, + hide_colorbar=False, vmin=None, vmax=None, log_scale=True, + dpi=200, output_dir=None): + """Plots the polar heatmap describing azimuth and latitude / elevation components. + + Plots the polar heatmap where each cell of the heatmap corresponds to + the specific element of the array provided by `gather_polar_errors` + function. + + Parameters + ---------- + data : 2D array + Indicates the array containing the sum of angular errors within the + specified angular ranges. It is usually provided by `gather_polar_errors` + function. + + name : str + Indicates the name of the output png file. + + interp_factor : float + Indicates the interpolation factor of the heatmap. + + color_limits : boolean + Specifies if the determined intensity limits should be returned. + + hide_colorbar : boolean + Specifies if the colorbar should be hidden. + + vmin : float + Indicates the minimum value of the colorbar. + + vmax : float + Indicates the maximum value of the colorbar. + + log_scale : float + Specifies if the heatmap sould be in the logarithmic scale. + + dpi : integer + Indicates the DPI of the output image. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + """ + th0, th1 = 0., 180. + r0, r1 = 0, 90 + thlabel, rlabel = 'Azimuth', 'Elevation' + + tr_scale = Affine2D().scale(np.pi/180., 1.) + tr = tr_scale + PolarAxes.PolarTransform() + + lat_ticks = [(.0*90., '0$^{\circ}$'), + (.33*90., '30$^{\circ}$'), + (.66*90., '60$^{\circ}$'), + (1.*90., '90$^{\circ}$')] + r_grid_locator = FixedLocator([v for v, s in lat_ticks]) + r_grid_formatter = DictFormatter(dict(lat_ticks)) + + angle_ticks = [(0*180., '90$^{\circ}$'), + (.25*180., '45$^{\circ}$'), + (.5*180., '0$^{\circ}$'), + (.75*180., '-45$^{\circ}$'), + (1.*180., '-90$^{\circ}$')] + theta_grid_locator = FixedLocator([v for v, s in angle_ticks]) + theta_tick_formatter = DictFormatter(dict(angle_ticks)) + + grid_helper = GridHelperCurveLinear(tr, + extremes=(th0, th1, r0, r1), + grid_locator1=theta_grid_locator, + grid_locator2=r_grid_locator, + tick_formatter1=theta_tick_formatter, + tick_formatter2=r_grid_formatter) + + fig = plt.figure() + ax = floating_axes.FloatingSubplot(fig, 111, grid_helper=grid_helper) + fig.add_subplot(ax) + ax.set_facecolor('white') + + ax.axis["bottom"].set_visible(False) + ax.axis["top"].toggle(ticklabels=True, label=True) + ax.axis["top"].set_axis_direction("bottom") + ax.axis["top"].major_ticklabels.set_axis_direction("top") + ax.axis["top"].label.set_axis_direction("top") + + ax.axis["left"].set_axis_direction("bottom") + ax.axis["right"].set_axis_direction("top") + + ax.axis["top"].label.set_text(thlabel) + ax.axis["left"].label.set_text(rlabel) + + aux_ax = ax.get_aux_axes(tr) + aux_ax.patch = ax.patch + ax.patch.zorder = 0.9 + + rad = np.linspace(0, 90, data.shape[1]) + azm = np.linspace(0, 180, data.shape[0]) + + f = interpolate.interp2d(rad, azm, data, kind='linear', bounds_error=True, fill_value=0) + + new_rad = np.linspace(0, 90, 180*interp_factor) + new_azm = np.linspace(0, 180, 360*interp_factor) + new_data_angle_dist = f(new_rad, new_azm) + new_r, new_th = np.meshgrid(new_rad, new_azm) + new_data_angle_dist += 1. + + if log_scale: + data_mesh = aux_ax.pcolormesh(new_th, new_r, new_data_angle_dist, cmap='jet', + norm=colors.LogNorm(vmin=1. if vmin is None else vmin, + vmax=new_data_angle_dist.max() if vmax is None else vmax)) + else: + data_mesh = aux_ax.pcolormesh(new_th, new_r, new_data_angle_dist, cmap='jet', vmin=vmin, vmax=vmax) + + cbar = plt.colorbar(data_mesh, orientation='vertical', shrink=.88, pad=.1, aspect=15) + cbar.ax.set_ylabel('Absolute error, [deg.]') + + if hide_colorbar: + cbar.remove() + + ax.grid(False) + + plt.show() + + if output_dir is not None: + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + fig.savefig(os.path.join(output_dir, '{}_chart.png'.format(name)), + transparent=False, + bbox_inches='tight', + pad_inches=0.1, + dpi=dpi) + + if color_limits: + return 1., new_data_angle_dist.max() + + +def plot_histogram_fig(data, num_bins, xticks, color, splot_index=111, output_dir=None, + xlim=(None, None), ylim=(None, None), name=None, in_percent=False, + bar_width=0.8, ticks_pad=7, xticks_fontsize=22, + yticks_fontsize=22, xlabel=None, ylabel=None, + labels_fontsize=20, grid_alpha=0.3, title_fontsize=22, + exp_fontsize=15, type=None, figsize=(12, 8), dpi=200): + """Plots the histogram from a given data. + + Parameters + ---------- + data : 1D array + Indicates the array containing the values. + + num_bins : integer + Indicates the number of histogram bins. + + xticks : array + Indicates the array of ticks of the X-axis. + + color : str + Indicates the color of the histogram. + + Returns + ------- + tuple : tuple of Figure and Axis objects. + """ + fig = plt.figure(figsize=figsize) + ax = fig.add_subplot(splot_index) + + weights = np.ones_like(data)/float(len(data)) + + n, bins, patches = ax.hist(data, num_bins, color=color, + rwidth=bar_width, weights=weights) + + ax.tick_params(axis='x', + labelsize=xticks_fontsize, + colors='#000000', + which='both', + direction='out', + length=6, + width=2, + pad=ticks_pad) + ax.tick_params(axis='x', + labelsize=xticks_fontsize, + colors='#000000', + which='minor', + direction='out', + length=4, + width=2, + pad=ticks_pad) + + ax.tick_params(axis='y', + labelsize=yticks_fontsize, + colors='#000000', + which='major', + direction='out', + length=6, + width=2, + pad=ticks_pad) + ax.tick_params(axis='y', + labelsize=xticks_fontsize, + colors='#000000', + which='minor', + direction='out', + length=4, + width=2, + pad=ticks_pad) + + ax.xaxis.set_minor_locator(ticker.MultipleLocator(5)) + ax.yaxis.set_minor_locator(ticker.MultipleLocator(5)) + + plt.xticks(xticks, xticks) + + ax.xaxis.offsetText.set_fontsize(exp_fontsize) + ax.yaxis.offsetText.set_fontsize(exp_fontsize) + + ax.ticklabel_format(axis='y', style='sci', scilimits=(-2, 2), + fontsize=yticks_fontsize) + ax.yaxis.set_major_locator(ticker.MaxNLocator(nbins=6, min_n_ticks=6)) + + ax.xaxis.grid(False) + ax.yaxis.grid(False) + + ax.set_ylabel(ylabel, labelpad=2, fontsize=labels_fontsize, color='black') + ax.set_xlabel(xlabel, labelpad=2, fontsize=labels_fontsize, color='black') + + ax.set_ylim(ylim) + ax.set_xlim(xlim) + + if in_percent: + vals = ax.get_yticks() + ax.set_yticklabels(['{:3.0f}'.format(x*100) for x in vals]) + + ax.set_axisbelow(True) + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + plt.tight_layout() + fig.savefig(os.path.join(output_dir, '{}_chart.png'.format(name)), + transparent=False, bbox_inches='tight', pad_inches=0.1, dpi=dpi) + + return fig, ax + + +def create_pie_chart(data, rngs, colors=['#244268', '#426084', '#67809F', '#95A9C1', '#C6D2E0'], + unit_scale=1.0, measure_quantity='m^3', figsize=(33, 15), + legend_loc=(0.383, -0.25), zebra_color=(False, 3), + legend_fontsize=50, chart_fontsize=60, dpi=72, name=None, + output_dir=None): + """Plots the piechart of from a given data. + + Parameters + ---------- + data : 1D array + Indicates the array containing the values. + + rngs : tuple of tuples + Indicates the ranges of the piechart. + + colors : array + Indicates the color for the region of the piechart corresponding to the + specific range. + + unit_scale : float + Indicates the scale factor of the data values. + + measure_quantity : str + Indicates the name of measure of the values. + + figsize : tuple of integers + Indicates the size of the output figure. + + legend_loc : tuple + Indicates the position of the legend of the figure. + + zebra_color : tuple + Allows to change the text color of the region to white from the first to + the speficied index of the region (True, reg_index). + + legend_fontsize : integer + Indicates the fontsize of the legend. + + chart_fontsize : integer + Indicates the fontsize of the figure. + + dpi : integer + Indicates the DPI of the output image. + + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + """ + def number(val): + if val < 1000: + return '%d' % val + + sv = str(val) + return '$\mathregular{10^{%d}}$' % (len(sv)-2) if val % 10 == 0 else '%0.0e' % val + + def get_title(v1, v2, measure_quantity): + ftm = '%s $\minus$ %s %s' + return ftm % (number(v1), number(v2), measure_quantity) + + data_ranges = [] + df = data * unit_scale + for rng in rngs: + rng_min, rng_max = rng[0], rng[1] + data_rng = df[(df > rng_min) & (df < rng_max)] + data_ranges.append(data_rng) + + num_elem = [len(p) for p in data_ranges] + se = sum(num_elem) + + print('Num of particles: {}'%(se)) + + proc_particles = [n/float(se) * 100.0 for n in num_elem] + + for size, rng in zip(num_elem, rngs): + print('{}-{}: {}'%(rng[0], rng[1], size)) + + titles = [get_title(minv, maxv, measure_quantity) for minv, maxv in rngs] + + textprops = {'fontsize': chart_fontsize, + 'weight': 'normal', + 'family': 'sans-serif'} + pie_width = 0.5 + fig, ax = plt.subplots(figsize=figsize) + ax.axis('equal') + patches, texts, autotexts = ax.pie(proc_particles, + textprops=textprops, + colors=colors, + autopct='%1.1f%%', + radius=1, + pctdistance=1-pie_width/2) + + if (zebra_color is not None) and (zebra_color[0]): + for tt in autotexts[:zebra_color[1]]: + tt.set_color('white') + + plt.setp(patches, + width=pie_width, + edgecolor='white') + + plt.legend(patches, titles, loc=legend_loc, fontsize=legend_fontsize) + + _d, _offset, _di = [1, -1], [0.45, 0.45], 0 + + for t, p in zip(autotexts, proc_particles): + if p < 2.0: + pos = list(t.get_position()) + pos[0] = pos[0] + _d[_di] * _offset[_di] + + t.set_position(pos) + _di += 1 + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + plt.tight_layout() + fig.savefig(os.path.join(output_dir, '{}_chart.png'.format(name)), + bbox_inches='tight', transparent=True, pad_inches=0.1, dpi=dpi) + + +def _bbox_3D(img): + """Crops the non-zero part of a volume + """ + r = np.any(img, axis=(1, 2)) + c = np.any(img, axis=(0, 2)) + z = np.any(img, axis=(0, 1)) + + rmin, rmax = np.where(r)[0][[0, -1]] + cmin, cmax = np.where(c)[0][[0, -1]] + zmin, zmax = np.where(z)[0][[0, -1]] + + return rmin, rmax, cmin, cmax, zmin, zmax + + +def plot_3d_orientation_map(name, lat_data, azth_data, radius_structure_elem=1, + output_dir=None, width=512, height=512, + camera_azth=44.5, camera_elev=35.8, camera_roll=0.0, + camera_fov=35.0, camera_zoom=0.0035, + camera_loc=(67.0, 81.6, 45.2), xlabel='', ylabel='', + zlabel='', axis_color='w', background_color='k'): + """Renders orientation data in 3D with RGB angular color-coding. + + Parameters + ---------- + name : str + Indicates the name of the output png file. + + lat_data : 3D array + Indicates the 3D array containing latitude / elevation angle at every point of + the skeleton in radians. + + azth_data : 3D array + Indicates the 3D array containing azimuth angle at every point of the skeleton + in radians. + + radius_structure_elem : integer + Indicates the size of the structure element of the dilation process to + thicken the skeleton. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + + width : int + Indicates the width of the visualization window. + + height : int + Indicates the width of the visualization window. + + camera_azth : float + Indicates the azimuth angle of the camera. + + camera_elev : float + Indicates the latitude / elevation angle of the camera. + + camera_roll : float + Indicates the roll angle of the camera. + + camera_fov : float + Indicates the field of view of the camera. + + camera_zoom : float + Indicates the zoom level of the camera. + + camera_loc : tuple + Indicates the camera location. + + xlabel : str + Indicates the label along the x-axis. + + ylabel : str + Indicates the label along the y-axis. + + zlabel : str + Indicates the label along the z-axis. + + axis_color : str + Indicates the color of axes. + + background_color : str + Indicates the background color of the figure. + """ + if not visvis_available: + print('The visvis package is not found. The visualization cannot be done.') + return + + rmin, rmax, cmin, cmax, zmin, zmax = _bbox_3D(azth_data) + + azth, lat = azth_data[rmin:rmax, cmin:cmax, zmin:zmax], \ + np.abs(lat_data[rmin:rmax, cmin:cmax, zmin:zmax]) + + skel = azth.copy().astype(np.float32) + skel[skel.nonzero()] = 1. + + azth = ndi.grey_dilation(azth, structure=morphology.ball(radius_structure_elem)) + lat = ndi.grey_dilation(lat, structure=morphology.ball(radius_structure_elem)) + skel = ndi.binary_dilation(skel, structure=morphology.ball(radius_structure_elem)) + + Z, Y, X = skel.nonzero() + vol_orient = np.zeros(skel.shape + (3,), dtype=np.float32) + + print(vol_orient.size, vol_orient[skel.nonzero()].size) + + for z, y, x in zip(Z, Y, X): + vol_orient[z, y, x] = geo2rgb(lat[z, y, x], azth[z, y, x]) + + app = vv.use() + + fig = vv.figure() + fig._currentAxes = None + fig.relativeFontSize = 2. + fig.position.w = width + fig.position.h = height + + t = vv.volshow(vol_orient[:, :, :], renderStyle='iso') + t.isoThreshold = 0.5 + + a = vv.gca() + a.camera.azimuth = camera_azth + a.camera.elevation = camera_elev + a.camera.roll = camera_roll + a.camera.fov = camera_fov + a.camera.zoom = camera_zoom + a.camera.loc = camera_loc + + a.bgcolor = background_color + a.axis.axisColor = axis_color + a.axis.xLabel = xlabel + a.axis.yLabel = ylabel + a.axis.zLabel = zlabel + + # def mouseUp(event): + # print 'mouseUp!!' + # a = vv.gca() + # print a.camera.GetViewParams() + # + # a.eventMouseUp.Bind(mouseUp) + # fig.eventMouseUp.Bind(mouseUp) + # + # a.Draw() + # fig.DrawNow() + + if output_dir is not None: + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + vv.screenshot(os.path.join(output_dir, '{}_3d_orientation.png'.format(name)), + sf=1, bg=background_color) + app.Run() + + +def plot_3d_diameter_map(name, data, unit_scale=1.0, measure_quantity='vox', + radius_structure_elem=1, output_dir=None, width=512, + height=512, camera_azth=44.5, camera_elev=35.8, + camera_roll=0.0, camera_fov=35.0, camera_zoom=0.0035, + camera_loc=(67.0, 81.6, 45.2), xlabel='', ylabel='', + zlabel='', axis_color='w', background_color='k', + cb_x_offset=10): + """Renders orientation data in 3D with RGB angular color-coding. + + Parameters + ---------- + name : str + Indicates the name of the output png file. + + data : 3D array + Indicates the 3D array containing diameter at every point of the skeleton. + + unit_scale : float + Indicates the scale factor of the data values. + + measure_quantity : str + Indicates the name of measure of the values. + + radius_structure_elem : integer + Indicates the size of the structure element of the dilation process to + thicken the skeleton. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + + camera_azth : float + Indicates the azimuth angle of the camera. + + width : int + Indicates the width of the visualization window. + + height : int + Indicates the width of the visualization window. + + camera_elev : float + Indicates the latitude / elevation angle of the camera. + + camera_roll : float + Indicates the roll angle of the camera. + + camera_fov : float + Indicates the field of view of the camera. + + camera_zoom : float + Indicates the zoom level of the camera. + + camera_loc : tuple + Indicates the camera location. + + xlabel : str + Indicates the label along the x-axis. + + ylabel : str + Indicates the label along the y-axis. + + zlabel : str + Indicates the label along the z-axis. + + axis_color : str + Indicates the color of axes. + + background_color : str + Indicates the background color of the figure. + + cb_x_offset : int + Indicates the offset of the colorbar from the right window side. + """ + if not visvis_available: + print('The visvis package is not found. The visualization cannot be done.') + return + + rmin, rmax, cmin, cmax, zmin, zmax = _bbox_3D(data) + dmtr = data[rmin:rmax, cmin:cmax, zmin:zmax] * unit_scale + skel = np.zeros_like(dmtr, dtype=np.uint8) + skel[dmtr.nonzero()] = 1 + + dmtr = ndi.grey_dilation(dmtr, structure=morphology.ball(radius_structure_elem)) + skel = ndi.binary_dilation(skel, structure=morphology.ball(radius_structure_elem)).astype(np.float32) + skel[skel.nonzero()] = 1. + + dmtr = dmtr * skel + + app = vv.use() + + fig = vv.figure() + fig._currentAxes = None + fig.relativeFontSize = 2. + fig.position.w = width + fig.position.h = height + + t = vv.volshow(dmtr[:, :, :], renderStyle='iso') + t.isoThreshold = 0.5 + t.colormap = vv.CM_JET + + a = vv.gca() + a.camera.azimuth = camera_azth + a.camera.elevation = camera_elev + a.camera.roll = camera_roll + a.camera.fov = camera_fov + a.camera.zoom = camera_zoom + a.camera.loc = camera_loc + + a.bgcolor = background_color + a.axis.axisColor = axis_color + a.axis.xLabel = xlabel + a.axis.yLabel = ylabel + a.axis.zLabel = zlabel + + cb = vv.colorbar() + cb.SetLabel('Diameter, [{}]'.format(measure_quantity)) + cb._label.position.x += cb_x_offset + + if output_dir is not None: + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + vv.screenshot(os.path.join(output_dir, '{}_3d_diameter.png'.format(name)), + sf=1, bg='w') + app.Run() + + +def plot_color_wheel(name, output_dir=None, dpi=500, xlabel='Elevation', + ylabel='Azimuth', fontsize=10, num_xticks=4, yticks=(-90, 90)): + """Plots the color wheel for visualizations of 3D orintation. + + Parameters + ---------- + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + + dpi : integer + Indicates the DPI of the output image. + + xlabel : str + Indicates the text along the x-axis. + + ylabel : str + Indicates the text along the y-axis. + + fontsize : int + Indicates the font size of labels along axes. + + num_xticks : int + Indicates the number of ticks along axes. + + yticks : tuple + Indicates the range of minimum and maximum values along the y-axis. + """ + + azth, lat = np.linspace(0., 1., num=180), np.linspace(0., 1., num=90) + rgb_arr = np.zeros((len(azth), len(lat), 3)) + for i in xrange(len(azth)): + for j in xrange(len(lat)): + rgb_arr[i, j, :] = colors.hsv_to_rgb([azth[i], lat[j], 1.0]) + + fig, ax = plt.subplots(figsize=(2, 2)) + ax.set_facecolor('red') + + ax.set_xlim([0, 90]) + ax.set_ylim([0, 181]) + + ax.set_yticks(np.linspace(0, 180, num=7).astype(np.int32)) + ax.set_yticklabels(np.linspace(yticks[0], yticks[1], num=7).astype(np.int32)) + + ax.set_xlabel(xlabel, fontsize=fontsize, labelpad=0, color='w') + ax.set_ylabel(ylabel, fontsize=fontsize, labelpad=0, color='w') + + xmajorlocator = ticker.LinearLocator(num_xticks) + ax.xaxis.set_major_locator(xmajorlocator) + ax.tick_params(direction='out', length=2, width=0, labelsize=fontsize, pad=0, colors='w') + + ax.imshow(rgb_arr) + + if output_dir is not None: + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + fig.savefig(os.path.join(output_dir, '{}_color_bar.png'.format(name)), + transparent=True, + bbox_inches='tight', + pad_inches=0.1, + dpi=dpi) + plt.show() + +def plot_fourier_orientation(data, orient_blocks, block_shape, figsize=(12,12), + cmap='gray', line_length=20, line_width=2.5, + line_color='red', line_style='-', name=None, + output_dir=None, dpi=200): + """Plots the orientation map in a block-wise manner. + + Plots the orientation vector over the image `data` at the center of each + block of the subdivided image. The orientation at every block 'orient_blocks' + and its size `block_shape` are specified by orientation estimation method. + The result can be and stored as a png file with DPI of `dpi` to the folder + specified by `output_dir`. + + Parameters + ---------- + data : ndarray + An image on top of which, the orientation will be plotted.. + + orient_blocks : ndarray + 2D array of orientation at every block of the subdivided image. + + block_shape : tuple + Indicates the block size within which the orientation is calculated. + + figsize : tuple of integers + Indicates the size of the output figure. + + cmap : str + Indicates the name of a colormap used for image. + + line_length : integer + Indicates the length of the orientation vector at each block. + + line_width : float + Indicates the line width of the orientation vector at each block. + + line_color : str + Indicates the line color of the orientation vector at each block. + + line_style : str + Indicates the line style of the orientation vector at each block. + + name : str + Indicates the name of the output png file. + + output_dir : str + Indicates the path to the output folder where the image will be stored. + + dpi : integer + Indicates the DPI of the output image. + """ + fig, ax = plt.subplots(figsize=figsize) + + ax.xaxis.grid(False) + ax.yaxis.grid(False) + ax.set_axis_off() + + ax.imshow(data, cmap=cmap) + + for i in xrange(orient_blocks.shape[0]): + for j in xrange(orient_blocks.shape[1]): + y0, x0 = block_shape[0] * j + block_shape[0]/2, \ + block_shape[1] * i + block_shape[1]/2 + + orientation = orient_blocks[i ,j] + + x2 = x0 - math.sin(orientation) * line_length + y2 = y0 - math.cos(orientation) * line_length + + x3 = x0 + math.sin(orientation) * line_length + y3 = y0 + math.cos(orientation) * line_length + + ax.plot((x2, x3), (y2, y3), linestyle=line_style, + linewidth=line_width, color=line_color) + + if (output_dir is not None) and (name is not None): + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + plt.tight_layout() + fig.savefig(os.path.join(output_dir, '{}_fourier_orientation_map.png'.format(name)), + bbox_inches='tight', transparent=True, pad_inches=0.05, dpi=dpi) + + plt.show() diff --git a/quanfima/requirements.txt b/quanfima/requirements.txt new file mode 100644 index 0000000..40dca1a --- /dev/null +++ b/quanfima/requirements.txt @@ -0,0 +1,7 @@ +matplotlib==2.0.2 +nose==1.3.7 +numpy==1.13.3 +pandas==0.19.2 +scikit-image==0.12.3 +scikit-learn==0.18.1 +scipy==0.19.0 diff --git a/quanfima/setup.py b/quanfima/setup.py new file mode 100644 index 0000000..e40cd36 --- /dev/null +++ b/quanfima/setup.py @@ -0,0 +1,36 @@ +from setuptools import setup + +def readme(): + with open('README.rst') as f: + return f.read() + +setup(name='quanfima', + version='0.1a3', + description='The package for morphological analysis and visualization of fibrous materials.', + long_description=readme(), + classifiers=[ + 'Development Status :: 3 - Alpha', + 'License :: OSI Approved :: MIT License', + 'Programming Language :: Python :: 2.7', + 'Topic :: Scientific/Engineering :: Information Analysis', + ], + keywords='biomaterials fiber analysis material science visualization', + url='http://github.com/rshkarin/quanfima', + author='Roman Shkarin, Andrei Shkarin', + author_email='roman.v.shkarin@gmail.com, andrei.shkarin@gmail.com', + license='MIT', + packages=['quanfima', 'docs'], + python_requires='>=2.6,<3', + install_requires=[ + 'matplotlib==2.0.2', + 'numpy>=1.13.3', + 'pandas>=0.19.2', + 'scikit-image>=0.12.3', + 'scikit-learn>=0.18.1', + 'scipy>=0.19.0', + 'imageio>=2.3.0' + ], + test_suite='nose.collector', + tests_require=['nose', 'nose-cover3'], + include_package_data=True, + zip_safe=False) diff --git a/quanfima/test.txt b/quanfima/test.txt new file mode 100644 index 0000000..9daeafb --- /dev/null +++ b/quanfima/test.txt @@ -0,0 +1 @@ +test diff --git a/quanfima/tox.ini b/quanfima/tox.ini new file mode 100644 index 0000000..74b39f6 --- /dev/null +++ b/quanfima/tox.ini @@ -0,0 +1,3 @@ +[flake8] +ignore = W391 +max-line-length = 100 diff --git a/quanfima_torun.py b/quanfima_torun.py new file mode 100644 index 0000000..549b02f --- /dev/null +++ b/quanfima_torun.py @@ -0,0 +1,40 @@ +import numpy as np +from skimage import io, filters +import cv2 +from quanfima.quanfima import morphology as mrph +from quanfima.quanfima import visualizationDE as vis +from quanfima.quanfima import visualization_testQ as visb +import matplotlib.pyplot as plt +#img_seg = cv2.imread('Data/MAX_5p_10_2.tif',0) +from os import listdir +from os.path import isfile, join +import numpy +import cv2 + + +def fiber_visualization(): + plt.ioff() + mypath='Data/' + onlyfiles = [ f for f in listdir(mypath) if isfile(join(mypath,f)) ] + images = numpy.empty(len(onlyfiles), dtype=object) + names = numpy.empty(len(onlyfiles), dtype=object) + for n in range(0, len(onlyfiles)): + images[n] = cv2.imread( join(mypath,onlyfiles[n]),0 ) + names[n]= (onlyfiles[n]) + img = numpy.empty(len(images), dtype=object) + for n in range(0,len(images)): + img_seg= images[n] + data = img_seg + skeleton = img_seg + skeleton_thick = img_seg + #skeleton, skeleton_thick = utils.prepare_data(img_seg) + cskel, fskel, omap, dmap, ovals, dvals = \ + mrph.estimate_fiber_properties(data, skeleton) + # plot results + vis.plot_diameter_map(dmap, cskel, figsize=(10,10), cmap='gist_rainbow', + name='DIAMETER_estimate', + output_dir='output/') + vis.plot_orientation_map(omap, fskel, min_label=None, max_label=None, + figsize=(10,10), + name='ORIENTATION_estimate', + output_dir='output/') \ No newline at end of file diff --git a/quanfima_torun11.py b/quanfima_torun11.py new file mode 100644 index 0000000..30a199e --- /dev/null +++ b/quanfima_torun11.py @@ -0,0 +1,100 @@ +from __future__ import print_function +import numpy as np +from skimage import io, filters +from quanfima.quanfima import morphology as mrph +from quanfima.quanfima import visualization as vis +from quanfima.quanfima.visualization_testQ import plot_orientation_test_map +from quanfima.quanfima import utils2 +#from quanfima.quanfima import utils +from skimage.morphology import watershed +from skimage.feature import peak_local_max +from skimage import measure +from skimage.segmentation import random_walker +import matplotlib.pyplot as plt +from scipy import ndimage +import matplotlib.pyplot as plt +from skimage import data +from skimage.filters import threshold_otsu +import cv2 +import math +import time +import os +import itertools +import numpy as np +from skimage import feature, measure, filters +from skimage.util.shape import view_as_blocks +from scipy import ndimage as ndi +from scipy.spatial import distance +import vigra +import pandas as pd +from multiprocessing import Pool +from quanfima.quanfima import cuda_available + +if cuda_available: + import pycuda.autoinit + import pycuda.driver as cuda + from pycuda.compiler import SourceModule + from pycuda import gpuarray + + +def colormaps(orient_type,paddding=20,diameter_window_radius=12,window_radius=12): + img = cv2.imread('Data/segmented_square.png',0) + #img = np.invert(img) + #img_seg = img + #plt.imshow(img) + binary = img #>125 + #plt.imshow(binary) + img_seg=binary + #binary.ndim + skeleton = img_seg + data = skeleton + skeleton_thick=skeleton + #plt.imshow(skeleton) + fiber_mask= data + fiber_skel = skeleton + #paddding=50 + #window_radius=12, + #orient_type='pca' + #diameter_window_radius=12 + orientation_vals, diameter_vals = [], [] + + padded_fiber_skel = np.pad(fiber_skel, pad_width=(paddding,), mode='constant', constant_values=0) + padded_fiber_mask = np.pad(fiber_mask, pad_width=(paddding,), mode='constant', constant_values=0) + #plt.imshow(padded_fiber_skel) + pcy, pcx = (10*2+1)/2., (10*2+1)/2. + + ycords, xcords = np.nonzero(padded_fiber_skel) + + clear_padded_fiber_skel = padded_fiber_skel.copy() + method = feature.corner_harris(clear_padded_fiber_skel, sigma=1)#sigma=1.5 + #plt.imshow(method) + corner_points = feature.corner_peaks(method, min_distance=1)#min_distance =3 + image = img + #fig = plt.figure() + #plt.imshow(image) + # Convert coordinates to x and y lists + #y_corner,x_corner = zip(*corner_points) + #plt.plot(x_corner,y_corner,'o') # Plot corners + #if title: + #plt.title('title') + #plt.xlim(0,data.shape[1]) + #plt.ylim(data.shape[0],0) # Images use weird axes + #fig.set_size_inches(np.array(fig.get_size_inches()) * 1.5) + #plt.show() + #print "Number of corners:",len(corners) + for i, (yy, xx) in enumerate(corner_points): + clear_padded_fiber_skel[yy-2:yy+4, xx-2:xx+4] = np.zeros((6, 6), + dtype=clear_padded_fiber_skel.dtype) + from quanfima.quanfima import morphology7 as mrph + cskel, fskel, omap, dmap, ovals, dvals = \ + mrph.estimate_fiber_properties(data, skeleton) + #plt.imshow(clear_padded_fiber_skel) + from quanfima.quanfima.visualization_testQ import plot_orientM_map, plot_orientation_test_map + # plot results + plot_orientM_map(omap, cskel, min_label=u'0°', max_label=u'180°', + figsize=(10,10), + name='2d_polymer', + output_dir='output/') + vis.plot_diameter_map(dmap, cskel, figsize=(10,10), cmap='gist_rainbow', + name='2d_polymer', + output_dir='/path/to/output/dir') diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..14a39ad --- /dev/null +++ b/setup.py @@ -0,0 +1,2 @@ +jupyter nbextension enable --py --sys-prefix widgetsnbextension +conda install -c ukoethe vigra=1.11.1 anaconda --yes; diff --git a/src/angles_csv.csv b/src/angles_csv.csv new file mode 100644 index 0000000..592fccb --- /dev/null +++ b/src/angles_csv.csv @@ -0,0 +1,180 @@ +-89,, +-88,, +-87,, +-86,, +-85,, +-84,, +-83,, +-82,, +-81,, +-80,, +-79,, +-78,, +-77,, +-76,, +-75,, +-74,, +-73,, +-72,, +-71,, +-70,, +-69,, +-68,, +-67,, +-66,, +-65,, +-64,, +-63,, +-62,, +-61,, +-60,, +-59,, +-58,, +-57,, +-57,, +-56,, +-55,, +-54,, +-53,, +-52,, +-51,, +-50,, +-49,, +-48,, +-47,, +-46,, +-45,, +-44,, +-43,, +-42,, +-41,, +-40,, +-39,, +-38,, +-37,, +-36,, +-35,, +-34,, +-33,, +-32,, +-31,, +-30,, +-29,, +-28,, +-27,, +-26,, +-25,, +-24,, +-23,, +-22,, +-21,, +-20,, +-19,, +-18,, +-17,, +-16,, +-15,, +-14,, +-13,, +-12,, +-11,, +-10,, +-9,, +-8,, +-7,, +-6,, +-5,, +-4,, +-3,, +-2,, +-1,, +0,, +1,, +2,, +3,, +4,, +5,, +6,, +7,, +8,, +9,, +10,, +11,, +12,, +13,, +14,, +15,, +16,, +17,, +18,, +19,, +20,, +21,, +22,, +23,, +24,, +25,, +26,, +27,, +28,, +29,, +30,, +31,, +32,, +33,, +34,, +35,, +36,, +37,, +38,, +39,, +40,, +41,, +42,, +43,, +44,, +45,, +46,, +47,, +48,, +49,, +50,, +51,, +52,, +53,, +54,, +55,, +56,, +57,, +58,, +59,, +60,, +61,, +62,, +63,, +64,, +65,, +66,, +67,, +68,, +69,, +70,, +71,, +72,, +73,, +74,, +75,, +76,, +77,, +78,, +79,, +80,, +81,, +82,, +83,, +84,, +85,, +86,, +87,, +88,, +89,, diff --git a/src/fibers_P.py b/src/fibers_P.py new file mode 100644 index 0000000..b387e41 --- /dev/null +++ b/src/fibers_P.py @@ -0,0 +1,509 @@ +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals +from __future__ import division +from __future__ import absolute_import + +import argparse, re, os, glob, sys, pprint, itertools +import inspect + +import matplotlib +matplotlib.rcParams['backend'] = "Qt5Agg" +from matplotlib.figure import Figure +import matplotlib.pyplot as plt + + +import numpy as np +from src.fiberfit_model.helpers import debug +import numpy as np +from numpy.linalg import eig, inv +import scipy.ndimage +from pylab import * + +import ipywidgets as widgets +from IPython.display import display +#btn = widgets.Button(description='Press Me') +#display(btn) +#display(QPixmap(saved_images_dir_name + "/angDist_" + ".png").scaled(300, 400, Qt.KeepAspectRatio, Qt.SmoothTransformation)) +import skimage.io as io +from skimage import data +from skimage.filters import sobel +from skimage.morphology import watershed +from skimage import feature +from skimage.color import rgb2gray +from skimage.feature import corner_harris,corner_peaks +from skimage.exposure import equalize_hist +from skimage import measure +from skimage.draw import ellipse +from skimage.measure import label, regionprops +from skimage.transform import rotate +from skimage.feature import corner_orientations +from skimage.filters import threshold_otsu +from skimage.segmentation import clear_border +from skimage.measure import label +from skimage.morphology import closing, square +from skimage.measure import regionprops +from skimage.color import label2rgb +import seaborn +import pandas +try: + from runtime import * +except ImportError: + from smop.runtime import * + +from src.fiberfit_model.helpers import debug +import numpy as np +from numpy.linalg import eig, inv +import scipy.ndimage +from pylab import * +from os import listdir +from os.path import isfile, join +import numpy +import cv2 + + + +import sys +import matplotlib +import numpy as np +import scipy +import scipy.ndimage +import scipy.interpolate +import scipy.optimize +import scipy.integrate +import scipy.stats +import math +import time +import glob +from pylab import * +from pandas import DataFrame +#from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas +import matplotlib.pyplot as plt +from src.foruploading5 import prepimages +from src.fiberfit_model.EllipseDirectFit import* +from src.fiberfit_model import helpers +from ipywidgets import IntProgress +from IPython.display import display +import time + +def fiber_analysisD(mypath,savepath, CO_lower,CO_upper,angleInc,radStep): + plt.ioff() + ## Upload Files + ## Upload Files + onlyfiles = [ f for f in listdir(mypath) if isfile(join(mypath,f)) ] + images = numpy.empty(len(onlyfiles), dtype=object) + names = numpy.empty(len(onlyfiles), dtype=object) + for n in range(0, len(onlyfiles)): + images[n] = cv2.imread( join(mypath,onlyfiles[n]) ) + names[n]= (onlyfiles[n]) + img = numpy.empty(len(images), dtype=object) + for n in range(0,len(images)): + img[n]= np.invert(images[n]) + ## Set up empty Matrices (can clean this up later ) + ell_data_test = numpy.empty(len(img),dtype=object) + Mirtheta1RadFinal_test =numpy.empty(len(img),dtype=object) + MirnormPower_test = numpy.empty(len(img),dtype=object) + normPower_test = numpy.empty(len(img),dtype=object) + theta1RadFinal_test = numpy.empty(len(img),dtype=object) + A_test = numpy.empty(len(img),dtype=object) + centroid_test = numpy.empty(len(img),dtype=object) + t_new_test = numpy.empty(len(img),dtype=object) + im = numpy.empty(len(img), dtype=object) + PF = numpy.empty(len(img),dtype=object) + PS = numpy.empty(len(img),dtype=object) + theta1Rad_test = numpy.empty(len(img),dtype=object) + PowerY_test = numpy.empty(len(img),dtype=object) + Mirtheta1RadFinal1_plot= [] + MirnormPower_plot=[] + th_plot=[] + r_line_plot=[] + normPower1_plot=[] + num=0 + index=num + + plt.rcParams.update({'figure.max_open_warning': 0}) ## avoid warning output + mpl.rcParams['figure.max_open_warning'] = 0 + maxiter = 100 + + ## begin analysis + print('Number of Images:',len(img)) + max_count=100 + tot_sum_progress = len(img) + PROGRESSBAR = IntProgress(min=0, max=max_count,description='Progress') + PROGRESSBAR.bar_style='info' + display(PROGRESSBAR) # display the bar + for n in range(0,len(img)): + PROGRESSBAR.value = int(n/tot_sum_progress*100) # signal to increment the progress bar + im=img[n][:,:,1] # select first image + M= im.shape[0] #get image shape + N1 = im.shape[1] + n1 = np.round(N1/ 2) - 1 + freq = np.arange(-n1, n1 + 1, 1) + x, y = freq, freq + + ## Fourier Transform + fft_result = np.fft.fft2(im) + Fshift = np.fft.fftshift(fft_result) + + ## adjust to power spectrum to remove imaginary components + magnitude_spectrum= 20*np.log(np.abs(Fshift)) + Pabs = np.abs(Fshift) ** 2 + + M= im.shape[0] + N1=im.shape[1] + + # shift in terms of image because power spectrum is the mirroR of lines so + # mirroring back in terms of image would give right allignment + PabsFlip1 = np.rot90(Pabs) + PabsFlip = np.flipud(PabsFlip1) + while (PabsFlip.shape[0]>x.shape[0]): + PabsFlip = np.delete(PabsFlip, (0), axis=0) + while (PabsFlip.shape[1]>y.shape[0]): + PabsFlip = np.delete(PabsFlip, (0), axis=1) + while (PabsFlip.shape[0] 0)[0]] + A = np.vstack((A1, dot(T, A1))) + A3 = A[3] - 2 * A[0] * centroid[0] - A[1] * centroid[1] + A4 = A[4] - 2 * A[2] * centroid[1] - A[1] * centroid[0] + A5 = A[5] + A[0] * centroid[0] ** 2 + A[2] * centroid[1] ** 2 + A[1] * centroid[0] * centroid[1] - A[3] * centroid[ + 0] - A[4] * centroid[1] + A[3] = A3 + A[4] = A4 + A[5] = A5 + A = A / linalg.norm(A) + if (abs(A[1]) < (1 * 10 ^ (-15))): + if ((A[0]) <= (A[2])): + # Ellipse is horizontal + angle = 0; + major = sqrt(1 / A[0]) + minor = sqrt(1 / A[2]) + else: + angle = np.pi/2; + major = sqrt(1 / A[2]) + minor = sqrt(1 / A[0]) + else: + R = ((A[2] - A[0]) / A[1]) + tg = R - sqrt((R * R) + 1) + angle = math.atan(tg) + P = (2 * tg) / (1 + (tg * tg)) + + if ((A[0] > 0 and A[1] > 0 and A[2] > 0)): + if (angle < (-pi / 4)): + angle = angle + np.pi + else: + angle = angle + elif ((A[1] / P <= (-A[1] / P))): + if (angle < 0): + angle = angle + np.pi / 2 + else: + angle = angle - np.pi / 2 + elif (A[0] < 0 and A[1] < 0 and A[2] < 0): + if (angle < 0): + angle = angle + np.pi + else: + angle = angle - np.pi + else: + # Switch + if (angle < 0): + angle = angle + np.pi / 2 + else: + angle = angle - np.pi / 2 + + t_New = (angle) #angle + t_new_test[n] = abs((angle* 180/np.pi)-90)-(2*angle* 180/np.pi) + Mirtheta1RadFinal1= Mirtheta1RadFinal_test[n] + #print('round2',Mirtheta1RadFinal1) + MirnormPower = MirnormPower_test[n] + #print('mirnormpower',MirnormPower) + t = t_new_test[n] + r_line = np.arange(0, max(MirnormPower) + .5, .5) + th = np.zeros(len(r_line)) + for i in range(0, len(r_line)): + th[i] = t + th = np.concatenate([th, (th + 180)]) + r_line = np.concatenate([r_line, r_line]) + ## Polar Graph + graphs = ["#71C4BE","#BBDBA0","#CE7E56","#BC5C81","#84355A"]; + palette = seaborn.set_palette(graphs); + seaborn.set_context("poster",font_scale=1.5); + seaborn.set_style('ticks') + csfont = {'fontname':'Times New Roman', + 'size':'14', + } + ticksfont = {'fontname':'Times New Roman'} + plt.figure(frameon=False, figsize=(15, 10)); + plt.polar(Mirtheta1RadFinal1, abs(MirnormPower), 'o',ms=5); + plt.suptitle(abs(t)-180, fontsize=20); + plt.polar(th * pi/180, r_line,color="#CE7E56",linewidth=4); + plt.savefig(savepath+'polarpositive'+names[n]); + if (max(MirnormPower)<2): + inc = 0.5 + elif (max(MirnormPower)<5): + inc = 1 + elif max(MirnormPower)<20: + inc = 5 + else: + inc = 10 + plt.yticks(np.arange(inc, max(MirnormPower), inc));#, **ticksfont) + + #plt.show() + t_final_rad = t*pi/180 + def fitted_func(thetas, c): + int_value, int_err = scipy.integrate.quadrature(func=lambda x: exp(c * cos(x)), a=0.0, b=np.pi) + return ((np.pi * (1.0 / np.pi * (int_value))) ** - 1) * \ + np.exp(c * np.cos(2 * (thetas - t_final_rad))) + c0 = 15 + xdata = theta1RadFinal_test[n] + theta1RadFinal = theta1RadFinal_test[n] + ydata= normPower_test[n] + normPower = normPower_test[n] + kappa, kappa_pcov = scipy.optimize.curve_fit(f=fitted_func, p0=(c0,), xdata=theta1RadFinal_test[n], ydata=normPower_test[n]) + + # Shift data for plotting purposes + + diff = abs(theta1RadFinal - (t * pi / 180)) + #print(diff,'diff') + #centerLoc=diff.find(diff==min(diff)) + centerLoc = diff[np.argmin(diff)]#find(diff== min(diff)) + #print(centerLoc,'centerLoc') + + num_angles = len(theta1RadFinal) + shift = ((round(num_angles / 2) - (num_angles - centerLoc))) + #print(shift,'shift') + + normPower1 = np.roll(normPower, -int(shift)) + normPower1 + theta1RadFinal1 = np.roll(theta1RadFinal, -int(shift)) + + if (shift > 0): + s = num_angles - shift + s= np.asscalar(s) + for k in range(s, num_angles): + theta1RadFinal1[k] = np.pi + theta1RadFinal1[k] + elif (shift < 0): + for k in range(0,90): + theta1RadFinal1[k] = -np.pi + theta1RadFinal1[k] + + ## Normalized Distribution Graph + cartDist = plt.figure(figsize=(10, 10)); # Creates a figure containing cartesian distribution. + + h2 = plt.bar((theta1RadFinal1 * 180 / pi), normPower1,color="#dedede"); + plt.xticks(np.arange(-360, 360, 45,), **ticksfont); + plt.xlim([-95, 95]); + p_act = fitted_func(theta1RadFinal1, kappa); + #print(p_act) + a = np.vstack([[p_act],[(theta1RadFinal1 * 180 / pi)]]); + a= a.T + b = a[a[:,1].argsort()]#np.sort(a,axis=1) + h3, = plt.plot(b[:,1],b[:,0], linewidth=4,color="#CE7E56"); + plt.xlabel('Angle (°)', **csfont); + plt.ylabel('Normalized Intensity', **csfont); + #print(max(normPower)) + if (max(normPower)<2): + inc = 0.5 + elif (max(normPower)<5): + inc = 1 + elif (max(normPower)<20): + inc = 5 + else: + inc = 10 + plt.yticks(np.arange(0, max(normPower1) + .3, inc), **ticksfont); + plt.ylim([0, max(normPower1) + .3]); + normPower1_plot.append(normPower1); + plt.savefig(savepath+'Normdist '+names[n]); + #plt.show() + #cartDist.savefig(dir + 'cartDist_' + number.__str__(), bbox_inches='tight') + plt.close(); + t_final=t + slope, intercept, rValue, pValue, stderr = scipy.stats.linregress(p_act, normPower1) + #print(rValue**2) + # Rounding results for Title of Figure + krnd = math.ceil(k * 1000) / 1000 + thrnd = abs(math.ceil(t_final * 1000) / 1000)-180 + krnd = math.trunc(krnd * 100) / 100 + thrnd = math.trunc(thrnd * 100) / 100 + a = 32.02 + b= -12.43 + c = 47.06 + d = -0.9185 + e = 19.43 + f = -0.07693 + x = kappa[0] + sig = math.exp(b*x) + c*math.exp(d*x) + e*exp(f*x) + #print("name",names[n],"kappa",kappa,"sigma",sig,"Rsq",rValue**2,"mu",t_final) + sum_NP = sum(normPower1) + #normPower2 = normPower1/sum(normPower1) + #sum_NP2 = sum(normPower2) + #print(sum_NP2) + ## Combined plot figure + palette = seaborn.set_palette(graphs); + seaborn.set_context("notebook",font_scale=1.5); + seaborn.set_style('ticks') + csfont = {'fontname':'Times New Roman', + 'size':'18', + } + ticksfont = {'fontname':'Times New Roman'} + fig = plt.figure(frameon=False,figsize=(20,12)); + fig.subplots_adjust(top=0.95); + fig.suptitle('%s '%(names[n]), fontsize=34, fontweight='bold'); + grid = plt.GridSpec(2, 5,hspace=0,wspace=0,width_ratios=[1, 1, 0.3, 1,1]); + cart_dist = plt.subplot(grid[1, 0:2]); + cart_dist.set_title('Von Mises Distribution (k)',fontsize=20); + h2 = cart_dist.bar((theta1RadFinal1 * 180 / pi), normPower1,color="#dedede"); + cart_dist.set_xlim([-92, 92]); + cart_dist.set_xlabel('Angle (°)',fontsize=18); + cart_dist.set_ylim([0, max(normPower1) + .3]); + cart_dist.set_ylabel('Normalized Intensity',fontsize=18); + + p_act = fitted_func(theta1RadFinal1, kappa) + a = np.vstack([[p_act],[(theta1RadFinal1 * 180 / pi)]]) + a= a.T + b = a[a[:,1].argsort()]#np.sort(a,axis=1) + h3, = cart_dist.plot(b[:,1],b[:,0], linewidth=4,color="#CE7E56"); + + if (max(normPower)<2): + inc = 0.5 + elif (max(normPower)<5): + inc = 1 + elif (max(normPower)<20): + inc = 5 + else: + inc = 10 + + polar_dist = plt.subplot(grid[0:2,3:5],projection='polar'); + polar_dist.plot(Mirtheta1RadFinal1, abs(MirnormPower), 'o'); + polar_dist.plot(th * pi/180, r_line,color="#CE7E56",linewidth=4); + polar_dist.set_title('Ellipse Fit (u)',fontsize=20).set_position([.5, 1.08]); + fig.text(0.75, 0.05, ' \n Kappa: %s '%(krnd)+' mu: %s \n'%(round((abs(t_final)-180),2))+' R^2: %s '%(round((rValue**2),2))+' sigma: %s \n'%(round(sig,2)),fontsize=20, fontweight='bold', bbox={'facecolor': '#BC5C81', 'alpha': 0.5, 'pad': 0.4}); + #polar_dist.text(np.radians(-110),1.95, ' \n Kappa: %s '%(krnd)+' mu: %s \n'%(round((abs(t_final)-180),2))+' R^2: %s '%(round((rValue**2),2))+' sigma: %s \n'%(round(sig,2)),fontsize=20, fontweight='bold', bbox={'facecolor': '#BC5C81', 'alpha': 0.5, 'pad': 0.4}) + orig_fig = plt.subplot(grid[0, 0]); + orig_fig.set_title('Original Image',fontsize=20).set_position([.5, 1.0]); + orig_fig.axis('off'); + fft_fig = plt.subplot(grid[0, 1]); + fft_fig.set_title('Fast Fourier Transform',fontsize=20).set_position([.5, 1.0]); + fft_fig.set_axis_off(); + fft_fig.axis('off'); + #fft_fig.imshow(log(PabsFlip), cmap='gray') + #orig_fig.imshow(images[n],cmap=plt.cm.viridis ) + #plt.tight_layout() + plt.savefig(savepath+'Grouped'+names[n]); + #print(n,'n') + #plt.show() + polar_dist_2 = plt.figure(figsize=(10, 10)); + #print('linebeforeindex') + #print(index,'currentindex') + #Mirtheta1RadFinal1_plot[index]=Mirtheta1RadFinal1; + #MirnormPower_plot[index]=MirnormPower; + #th_plot[index] = th * pi/180 + #print(th_plot[index],'thplotindex') + #r_line_plot[index]=r_line + Mirtheta1RadFinal1_plot.append(Mirtheta1RadFinal1); + MirnormPower_plot.append(MirnormPower); + th_plot.append(th * pi/180) + r_line_plot.append(r_line) + polar_dist.plot(Mirtheta1RadFinal1, MirnormPower); + polar_dist.plot(th * pi/180, r_line,linewidth=4); + polar_dist.set_title('Ellipse Fit (u)',fontsize=20).set_position([.5, 1.08]); + plt.close('all') + #index=index+1 + #print(index,'Index') +## export of cartesian coordinate data + arr = numpy.array(MirnormPower_plot) + arr2 = numpy.array(Mirtheta1RadFinal1_plot) + arr3 = numpy.array(normPower1_plot) + arr4 = numpy.array(theta1RadFinal1*180/pi) + a = numpy.transpose(arr3) + b = numpy.transpose(arr4) + b[b> 100] = 0 + b[b<-100] = 0 + export_data=np.hstack(((b[:,None]),a)) + numpy.savetxt(savepath+"data_export.csv", export_data, delimiter=",") + return export_data + diff --git a/src/fibers_QQQQQQ_invert.py b/src/fibers_QQQQQQ_invert.py new file mode 100644 index 0000000..1bd36eb --- /dev/null +++ b/src/fibers_QQQQQQ_invert.py @@ -0,0 +1,527 @@ +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals +from __future__ import division +from __future__ import absolute_import + +import argparse, re, os, glob, sys, pprint, itertools +import inspect + +import matplotlib +matplotlib.rcParams['backend'] = "Qt5Agg" +from matplotlib.figure import Figure +import matplotlib.pyplot as plt + +import pandas as pd +import numpy as np +from src.fiberfit_model.helpers import debug +import numpy as np +from numpy.linalg import eig, inv +import scipy.ndimage +from pylab import * + +import ipywidgets as widgets +from IPython.display import display +#btn = widgets.Button(description='Press Me') +#display(btn) +#display(QPixmap(saved_images_dir_name + "/angDist_" + ".png").scaled(300, 400, Qt.KeepAspectRatio, Qt.SmoothTransformation)) +import skimage.io as io +from skimage import data +from skimage.filters import sobel +from skimage.segmentation import watershed +#from skimage.morphology import watershed +from skimage import feature +from skimage.color import rgb2gray +from skimage.feature import corner_harris,corner_peaks +from skimage.exposure import equalize_hist +from skimage import measure +from skimage.draw import ellipse +from skimage.measure import label, regionprops +from skimage.transform import rotate +from skimage.feature import corner_orientations +from skimage.filters import threshold_otsu +from skimage.segmentation import clear_border +from skimage.measure import label +from skimage.morphology import closing, square +from skimage.measure import regionprops +from skimage.color import label2rgb +import seaborn +import pandas +try: + from runtime import * +except ImportError: + from smop.runtime import * + +from src.fiberfit_model.helpers import debug +import numpy as np +from numpy.linalg import eig, inv +import scipy.ndimage +from pylab import * +from os import listdir +from os.path import isfile, join +import numpy +import cv2 + + + +import sys +import matplotlib +import numpy as np +import scipy +import scipy.ndimage +import scipy.interpolate +import scipy.optimize +import scipy.integrate +import scipy.stats +import math +import time +import glob +from pylab import * +from pandas import DataFrame +#from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas +import matplotlib.pyplot as plt +from src.foruploading5 import prepimages +from src.fiberfit_model.EllipseDirectFit import* +from src.fiberfit_model import helpers +from ipywidgets import IntProgress +from IPython.display import display +import time + +def fiber_analysis(mypath,savepath, CO_lower,CO_upper): + angleInc = 1 + radStep = 0.5 + ##%matplotlib inline + ##plt.ioff() + ## Upload Files + ## Upload Files + onlyfiles = [ f for f in listdir(mypath) if isfile(join(mypath,f)) ] + images = numpy.empty(len(onlyfiles), dtype=object) + names = numpy.empty(len(onlyfiles), dtype=object) + for n in range(0, len(onlyfiles)): + images[n] = cv2.imread( join(mypath,onlyfiles[n]) ) + names[n]= (onlyfiles[n]) + img = numpy.empty(len(images), dtype=object) + for n in range(0,len(images)): + img[n]= np.invert(images[n]) + ## Set up empty Matrices (can clean this up later ) + ell_data_test = numpy.empty(len(img),dtype=object) + Mirtheta1RadFinal_test =numpy.empty(len(img),dtype=object) + MirnormPower_test = numpy.empty(len(img),dtype=object) + normPower_test = numpy.empty(len(img),dtype=object) + theta1RadFinal_test = numpy.empty(len(img),dtype=object) + A_test = numpy.empty(len(img),dtype=object) + centroid_test = numpy.empty(len(img),dtype=object) + t_new_test = numpy.empty(len(img),dtype=object) + im = numpy.empty(len(img), dtype=object) + PF = numpy.empty(len(img),dtype=object) + PS = numpy.empty(len(img),dtype=object) + theta1Rad_test = numpy.empty(len(img),dtype=object) + PowerY_test = numpy.empty(len(img),dtype=object) + Mirtheta1RadFinal1_plot= [] + MirnormPower_plot=[] + th_plot=[] + r_line_plot=[] + normPower1_plot=[] + num=0 + index=num + + plt.rcParams.update({'figure.max_open_warning': 0}) ## avoid warning output + mpl.rcParams['figure.max_open_warning'] = 0 + maxiter = 100 + + ## begin analysis + print('Number of Images:',len(img)) + max_count=100 + tot_sum_progress = len(img) + PROGRESSBAR = IntProgress(min=0, max=max_count,description='Progress') + PROGRESSBAR.bar_style='info' + display(PROGRESSBAR) # display the bar + for n in range(0,len(img)): + PROGRESSBAR.value = int(n/tot_sum_progress*100) # signal to increment the progress bar + im=img[n][:,:,1] # select first image + M= im.shape[0] #get image shape + N1 = im.shape[1] + n1 = np.round(N1/ 2) - 1 + freq = np.arange(-n1, n1 + 1, 1) + x, y = freq, freq + + ## Fourier Transform + fft_result = np.fft.fft2(im) + Fshift = np.fft.fftshift(fft_result) + + ## adjust to power spectrum to remove imaginary components + magnitude_spectrum= 20*np.log(np.abs(Fshift)) + Pabs = np.abs(Fshift) ** 2 + + M= im.shape[0] + N1=im.shape[1] + + # shift in terms of image because power spectrum is the mirroR of lines so + # mirroring back in terms of image would give right allignment + PabsFlip1 = np.rot90(Pabs) + PabsFlip = np.flipud(PabsFlip1) + while (PabsFlip.shape[0]>x.shape[0]): + PabsFlip = np.delete(PabsFlip, (0), axis=0) + while (PabsFlip.shape[1]>y.shape[0]): + PabsFlip = np.delete(PabsFlip, (0), axis=1) + while (PabsFlip.shape[0] 0)[0]] + A = np.vstack((A1, dot(T, A1))) + A3 = A[3] - 2 * A[0] * centroid[0] - A[1] * centroid[1] + A4 = A[4] - 2 * A[2] * centroid[1] - A[1] * centroid[0] + A5 = A[5] + A[0] * centroid[0] ** 2 + A[2] * centroid[1] ** 2 + A[1] * centroid[0] * centroid[1] - A[3] * centroid[ + 0] - A[4] * centroid[1] + A[3] = A3 + A[4] = A4 + A[5] = A5 + A = A / linalg.norm(A) + if (abs(A[1]) < (1 * 10 ^ (-15))): + if ((A[0]) <= (A[2])): + # Ellipse is horizontal + angle = 0; + major = sqrt(1 / A[0]) + minor = sqrt(1 / A[2]) + else: + angle = np.pi/2; + major = sqrt(1 / A[2]) + minor = sqrt(1 / A[0]) + else: + R = ((A[2] - A[0]) / A[1]) + tg = R - sqrt((R * R) + 1) + angle = math.atan(tg) + P = (2 * tg) / (1 + (tg * tg)) + + if ((A[0] > 0 and A[1] > 0 and A[2] > 0)): + if (angle < (-pi / 4)): + angle = angle + np.pi + else: + angle = angle + elif ((A[1] / P <= (-A[1] / P))): + if (angle < 0): + angle = angle + np.pi / 2 + else: + angle = angle - np.pi / 2 + elif (A[0] < 0 and A[1] < 0 and A[2] < 0): + if (angle < 0): + angle = angle + np.pi + else: + angle = angle - np.pi + else: + # Switch + if (angle < 0): + angle = angle + np.pi / 2 + else: + angle = angle - np.pi / 2 + + t_New = (angle) #angle + t_new_test[n] = abs((angle* 180/np.pi)-90)-(2*angle* 180/np.pi) + Mirtheta1RadFinal1= Mirtheta1RadFinal_test[n] + #print('round2',Mirtheta1RadFinal1) + MirnormPower = MirnormPower_test[n] + #print('mirnormpower',MirnormPower) + t = t_new_test[n] + r_line = np.arange(0, max(MirnormPower) + .5, .5) + th = np.zeros(len(r_line)) + for i in range(0, len(r_line)): + th[i] = t + th = np.concatenate([th, (th + 180)]) + r_line = np.concatenate([r_line, r_line]) + ## Polar Graph + graphs = ["#71C4BE","#BBDBA0","#CE7E56","#BC5C81","#84355A"]; + palette = seaborn.set_palette(graphs); + seaborn.set_context("poster",font_scale=1.5); + seaborn.set_style('ticks') + csfont = {'fontname':'Times New Roman', + 'size':'14', + } + ticksfont = {'fontname':'Times New Roman'} + plt.figure(frameon=False, figsize=(15, 10)); + plt.polar(Mirtheta1RadFinal1, abs(MirnormPower), 'o',ms=5); + plt.suptitle(abs(t)-180, fontsize=20); + plt.polar(th * pi/180, r_line,color="#CE7E56",linewidth=4); + plt.savefig(savepath+'polarpositive'+names[n]); + if (max(MirnormPower)<2): + inc = 0.5 + elif (max(MirnormPower)<5): + inc = 1 + elif max(MirnormPower)<20: + inc = 5 + else: + inc = 10 + plt.yticks(np.arange(inc, max(MirnormPower), inc));#, **ticksfont) + + #plt.show() + t_final_rad = t*pi/180 + def fitted_func(thetas, c): + int_value, int_err = scipy.integrate.quadrature(func=lambda x: exp(c * cos(x)), a=0.0, b=np.pi) + return ((np.pi * (1.0 / np.pi * (int_value))) ** - 1) * \ + np.exp(c * np.cos(2 * (thetas - t_final_rad))) + c0 = 15 + xdata = theta1RadFinal_test[n] + theta1RadFinal = theta1RadFinal_test[n] + ydata= normPower_test[n] + normPower = normPower_test[n] + kappa, kappa_pcov = scipy.optimize.curve_fit(f=fitted_func, p0=(c0,), xdata=theta1RadFinal_test[n], ydata=normPower_test[n]) + + # Shift data for plotting purposes + + diff = abs(theta1RadFinal - (t * pi / 180)) + #print(diff,'diff') + #centerLoc=diff.find(diff==min(diff)) + centerLoc = diff[np.argmin(diff)]#find(diff== min(diff)) + #print(centerLoc,'centerLoc') + + num_angles = len(theta1RadFinal) + shift = ((round(num_angles / 2) - (num_angles - centerLoc))) + #print(shift,'shift') + + normPower1 = np.roll(normPower, -int(shift)) + normPower1 + theta1RadFinal1 = np.roll(theta1RadFinal, -int(shift)) + + if (shift > 0): + s = num_angles - shift + s= np.asscalar(s) + for k in range(s, num_angles): + theta1RadFinal1[k] = np.pi + theta1RadFinal1[k] + elif (shift < 0): + for k in range(0,90): + theta1RadFinal1[k] = -np.pi + theta1RadFinal1[k] + + ## Normalized Distribution Graph + cartDist = plt.figure(figsize=(10, 10)); # Creates a figure containing cartesian distribution. + + h2 = plt.bar((theta1RadFinal1 * 180 / pi), normPower1,color="#dedede"); + plt.xticks(np.arange(-360, 360, 45,), **ticksfont); + plt.xlim([-95, 95]); + p_act = fitted_func(theta1RadFinal1, kappa); + #print(p_act) + a = np.vstack([[p_act],[(theta1RadFinal1 * 180 / pi)]]); + a= a.T + b = a[a[:,1].argsort()]#np.sort(a,axis=1) + h3, = plt.plot(b[:,1],b[:,0], linewidth=4,color="#CE7E56"); + plt.xlabel('Angle (°)', **csfont); + plt.ylabel('Normalized Intensity', **csfont); + #print(max(normPower)) + if (max(normPower)<2): + inc = 0.5 + elif (max(normPower)<5): + inc = 1 + elif (max(normPower)<20): + inc = 5 + else: + inc = 10 + plt.yticks(np.arange(0, max(normPower1) + .3, inc), **ticksfont); + plt.ylim([0, max(normPower1) + .3]); + normPower1_plot.append(normPower1); + plt.savefig(savepath+'Normdist '+names[n]); + #plt.show() + #cartDist.savefig(dir + 'cartDist_' + number.__str__(), bbox_inches='tight') + plt.close(); + t_final=t + slope, intercept, rValue, pValue, stderr = scipy.stats.linregress(p_act, normPower1) + #print(rValue**2) + # Rounding results for Title of Figure + krnd = math.ceil(k * 1000) / 1000 + thrnd = abs(math.ceil(t_final * 1000) / 1000)-180 + krnd = math.trunc(krnd * 100) / 100 + thrnd = math.trunc(thrnd * 100) / 100 + a = 32.02 + b= -12.43 + c = 47.06 + d = -0.9185 + e = 19.43 + f = -0.07693 + x = kappa[0] + sig = math.exp(b*x) + c*math.exp(d*x) + e*exp(f*x) + #print("name",names[n],"kappa",kappa,"sigma",sig,"Rsq",rValue**2,"mu",t_final) + sum_NP = sum(normPower1) + #normPower2 = normPower1/sum(normPower1) + #sum_NP2 = sum(normPower2) + #print(sum_NP2) + ## Combined plot figure + palette = seaborn.set_palette(graphs); + seaborn.set_context("notebook",font_scale=1.5); + seaborn.set_style('ticks') + csfont = {'fontname':'Times New Roman', + 'size':'18', + } + ticksfont = {'fontname':'Times New Roman'} + fig = plt.figure(frameon=False,figsize=(20,12)); + fig.subplots_adjust(top=0.95); + fig.suptitle('%s '%(names[n]), fontsize=34, fontweight='bold'); + grid = plt.GridSpec(2, 5,hspace=0,wspace=0,width_ratios=[1, 1, 0.3, 1,1]); + cart_dist = plt.subplot(grid[1, 0:2]); + cart_dist.set_title('Von Mises Distribution (k)',fontsize=20); + h2 = cart_dist.bar((theta1RadFinal1 * 180 / pi), normPower1,color="#dedede"); + cart_dist.set_xlim([-92, 92]); + cart_dist.set_xlabel('Angle (°)',fontsize=18); + cart_dist.set_ylim([0, max(normPower1) + .3]); + cart_dist.set_ylabel('Normalized Intensity',fontsize=18); + + p_act = fitted_func(theta1RadFinal1, kappa) + a = np.vstack([[p_act],[(theta1RadFinal1 * 180 / pi)]]) + a= a.T + b = a[a[:,1].argsort()]#np.sort(a,axis=1) + h3, = cart_dist.plot(b[:,1],b[:,0], linewidth=4,color="#CE7E56"); + + if (max(normPower)<2): + inc = 0.5 + elif (max(normPower)<5): + inc = 1 + elif (max(normPower)<20): + inc = 5 + else: + inc = 10 + + polar_dist = plt.subplot(grid[0:2,3:5],projection='polar'); + polar_dist.plot(Mirtheta1RadFinal1, abs(MirnormPower), 'o'); + polar_dist.plot(th * pi/180, r_line,color="#CE7E56",linewidth=4); + polar_dist.set_title('Ellipse Fit (u)',fontsize=20).set_position([.5, 1.08]); + fig.text(0.75, 0.05, ' \n Kappa: %s '%(krnd)+' mu: %s \n'%(round((abs(t_final)-180),2))+' R^2: %s '%(round((rValue**2),2))+' sigma: %s \n'%(round(sig,2)),fontsize=20, fontweight='bold', bbox={'facecolor': '#BC5C81', 'alpha': 0.5, 'pad': 0.4}); + #polar_dist.text(np.radians(-110),1.95, ' \n Kappa: %s '%(krnd)+' mu: %s \n'%(round((abs(t_final)-180),2))+' R^2: %s '%(round((rValue**2),2))+' sigma: %s \n'%(round(sig,2)),fontsize=20, fontweight='bold', bbox={'facecolor': '#BC5C81', 'alpha': 0.5, 'pad': 0.4}) + orig_fig = plt.subplot(grid[0, 0]); + orig_fig.set_title('Original Image',fontsize=20).set_position([.5, 1.0]); + orig_fig.axis('off'); + fft_fig = plt.subplot(grid[0, 1]); + fft_fig.set_title('Fast Fourier Transform',fontsize=20).set_position([.5, 1.0]); + fft_fig.set_axis_off(); + fft_fig.axis('off'); + fft_fig.imshow(log(PabsFlip), cmap='gray') + orig_fig.imshow(images[n],cmap=plt.cm.viridis ) + #plt.tight_layout() + plt.savefig(savepath+'Grouped'+names[n]); + #print(n,'n') + #plt.show() + polar_dist_2 = plt.figure(figsize=(10, 10)); + #print('linebeforeindex') + #print(index,'currentindex') + #Mirtheta1RadFinal1_plot[index]=Mirtheta1RadFinal1; + #MirnormPower_plot[index]=MirnormPower; + #th_plot[index] = th * pi/180 + #print(th_plot[index],'thplotindex') + #r_line_plot[index]=r_line + Mirtheta1RadFinal1_plot.append(Mirtheta1RadFinal1); + MirnormPower_plot.append(MirnormPower); + th_plot.append(th * pi/180) + r_line_plot.append(r_line) + polar_dist.plot(Mirtheta1RadFinal1, MirnormPower); + polar_dist.plot(th * pi/180, r_line,linewidth=4); + polar_dist.set_title('Ellipse Fit (u)',fontsize=20).set_position([.5, 1.08]); + plt.close('all') + #index=index+1 + #print(index,'Index') + if n == (tot_sum_progress-1): + PROGRESSBAR.value = 100 # signal to increment the progress bar +## export of cartesian coordinate data + arr = numpy.array(MirnormPower_plot) + arr2 = numpy.array(Mirtheta1RadFinal1_plot) + arr3 = numpy.array(normPower1_plot) + cv2.imread('src/angles_csv.csv') + angles = np.genfromtxt('src/angles_csv.csv',delimiter=',') + arr4 = angles[:,0] + #arr4 = numpy.array(theta1RadFinal1*180/pi) + a = numpy.transpose(arr3) + b = numpy.transpose(arr4) + b[b> 100] = 0 + b[b<-100] = 0 + export_data=np.hstack(((b[:,None]),a)) + #numpy.savetxt(savepath+"data_export.csv", export_data, delimiter=",") + my_input=['Angle'] + for n in range(0,len(img)): + my_input.append(names[n]) + column_values = my_input + df = pd.DataFrame(data = export_data, + columns = column_values) + df.to_csv(savepath+"data_export.csv") + return export_data + diff --git a/src/foruploading5.py b/src/foruploading5.py new file mode 100644 index 0000000..96adb8b --- /dev/null +++ b/src/foruploading5.py @@ -0,0 +1,32 @@ +import io +import ipywidgets as widgets +import numpy +import cv2 +import os +from os import listdir +from os.path import isfile, join +import numpy as np + +def ff(b): + for elem in b.new.values(): + #print(elem['metadata']['name']) + with open(os.path.join('Data/',elem['metadata']['name']), 'wb') as file: + file.write( elem['content']) + +def upload(s): + s.style.button_color = "lightblue" + s.observe(ff, names='value') + +def prepimages(mypath): + ##Path to folder where you would like images to be saved + savepath = 'output2/' + ## Upload Files + onlyfiles = [ f for f in listdir(mypath) if isfile(join(mypath,f)) ]; + images = numpy.empty(len(onlyfiles), dtype=object); + names = numpy.empty(len(onlyfiles), dtype=object); + for n in range(0, len(onlyfiles)): + images[n] = cv2.imread( join(mypath,onlyfiles[n]) ); + names[n]= (onlyfiles[n]); + img_tot = numpy.empty(len(images), dtype=object); + for n in range(0,len(images)): + img_tot[n]= np.invert(images[n]) diff --git a/test_images/not_square.jpg b/test_images/not_square.jpg deleted file mode 100644 index 50f7d79..0000000 Binary files a/test_images/not_square.jpg and /dev/null differ diff --git a/widgetsmp.py b/widgetsmp.py new file mode 100644 index 0000000..061d9fa --- /dev/null +++ b/widgetsmp.py @@ -0,0 +1,130 @@ +import pandas as pd +import matplotlib.pyplot as plt +#%matplotlib notebook +from ipywidgets import * +from IPython.display import display +from IPython.html import widgets +import base64 +res = 'computed results' +plt.style.use('ggplot') +from src.fibers_QQQQQQ_invert import fiber_analysis +from quanfima_torun import fiber_visualization +import os +import glob +#button = widgets.Button(description='My Button') +#out = widgets.Output() +#import ipyvuetify as v +from IPython.display import display, FileLink +# linking button and function together using a button's method +#button.on_click(on_button_clicked) +# displaying button and its output together +#widgets.VBox([button,out]) + +from src.foruploading5 import ff,upload,prepimages +import ipywidgets as widgets +from IPython.display import display, Markdown, clear_output +from ipywidgets import interact, interact_manual +def s(_): + s = widgets.FileUpload(accept='', multiple=True) + return s.data + +def run(_): + + button = widgets.Button( + description='Run Fiber Analysis', + layout={'width': '300px'}) + out= widgets.Output() + def on_button_clicked(_): + # "linking function with output" + with out: + # what happens when we press the button + clear_output() + print('Program executing. Output in output folder on home page') + fiber_analysis( CO_lower=int(text1.value), + CO_upper=int(text2.value), + mypath='Data/', + savepath='output/') + text1 = widgets.Text( + placeholder='20', + description="Max fiber Thickness (pixels)", + layout={'width':'300px'}, + style={'description_width':'175px'}) + text2 = widgets.Text( + placeholder='2', + description="Min fiber Thickness (pixels)", + layout={'width':'300px'}, + style={'description_width':'175px'}) + button.on_click(on_button_clicked) + box=widgets.VBox([text1,text2,button]) + display(box,out) +def run2(_): + + button = widgets.Button( + description='Run Fiber Orientation Visualization', + layout={'width': '300px'}) + out= widgets.Output() + def on_button_clicked(_): + # "linking function with output" + with out: + # what happens when we press the button + clear_output() + print('Program executing. Output in output folder on home page') + fiber_visualization() + + button.on_click(on_button_clicked) + box=widgets.VBox([button]) + display(box,out) + +from ipywidgets import HTML +HTML() + +def save(_): + filename = 'notebook.tar.gz' + b64 = base64.b64encode(res.encode()) + payload = b64.decode() + #BUTTONS + html_buttons = ''' + + + + + + + + + + ''' + html_button = html_buttons.format(payload=payload,filename=filename) + display(HTML(html_button)) + out= widgets.Output() + + +def clear(_): + button = widgets.Button( + description='Clear Data', + layout={'width': '300px'}) + out= widgets.Output() + def on_button_clicked(_): + # "linking function with output" + with out: + # what happens when we press the button + clear_output() + print('Program executing. Files in Data and output folder will be cleared') + files = glob.glob('Data/*',recursive=True) + for f in files: + try: + os.remove(f) + except OSError as e: + print("Error: %s : %s" % (f, e.strerror)) + + files = glob.glob('output/*',recursive=True) + + for f in files: + try: + os.remove(f) + except OSError as e: + print("Error: %s : %s" % (f, e.strerror)) + button.on_click(on_button_clicked) + box=widgets.VBox([button]) + display(box,out) + out= widgets.Output()