{ "cells": [ { "attachments": {}, "cell_type": "markdown", "id": "64cdaa1c", "metadata": {}, "source": [ "# Process results" ] }, { "cell_type": "code", "execution_count": 1, "id": "12040cb6-7143-404a-8e0c-5da728d0a8a4", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Requirement already satisfied: vtk in i:\\iisc-project\\venv\\lib\\site-packages (9.2.6)\n", "Requirement already satisfied: matplotlib>=2.0.0 in i:\\iisc-project\\venv\\lib\\site-packages (from vtk) (3.8.0)\n", "Requirement already satisfied: contourpy>=1.0.1 in i:\\iisc-project\\venv\\lib\\site-packages (from matplotlib>=2.0.0->vtk) (1.0.7)\n", "Requirement already satisfied: cycler>=0.10 in i:\\iisc-project\\venv\\lib\\site-packages (from matplotlib>=2.0.0->vtk) (0.11.0)\n", "Requirement already satisfied: fonttools>=4.22.0 in i:\\iisc-project\\venv\\lib\\site-packages (from matplotlib>=2.0.0->vtk) (4.38.0)\n", "Requirement already satisfied: kiwisolver>=1.0.1 in i:\\iisc-project\\venv\\lib\\site-packages (from matplotlib>=2.0.0->vtk) (1.4.4)\n", "Requirement already satisfied: numpy<2,>=1.21 in i:\\iisc-project\\venv\\lib\\site-packages (from matplotlib>=2.0.0->vtk) (1.24.2)\n", "Requirement already satisfied: packaging>=20.0 in i:\\iisc-project\\venv\\lib\\site-packages (from matplotlib>=2.0.0->vtk) (23.0)\n", "Requirement already satisfied: pillow>=6.2.0 in i:\\iisc-project\\venv\\lib\\site-packages (from matplotlib>=2.0.0->vtk) (9.4.0)\n", "Requirement already satisfied: pyparsing>=2.3.1 in i:\\iisc-project\\venv\\lib\\site-packages (from matplotlib>=2.0.0->vtk) (3.0.9)\n", "Requirement already satisfied: python-dateutil>=2.7 in i:\\iisc-project\\venv\\lib\\site-packages (from matplotlib>=2.0.0->vtk) (2.8.2)\n", "Requirement already satisfied: six>=1.5 in i:\\iisc-project\\venv\\lib\\site-packages (from python-dateutil>=2.7->matplotlib>=2.0.0->vtk) (1.16.0)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "[notice] A new release of pip is available: 23.0 -> 24.0\n", "[notice] To update, run: python.exe -m pip install --upgrade pip\n" ] } ], "source": [ "!pip install vtk" ] }, { "cell_type": "code", "execution_count": 2, "id": "d78f6cc0", "metadata": { "vscode": { "languageId": "plaintext" } }, "outputs": [ { "ename": "ModuleNotFoundError", "evalue": "No module named 'cv2'", "output_type": "error", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[1;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)", "Cell \u001b[1;32mIn[2], line 18\u001b[0m\n\u001b[0;32m 16\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mmatplotlib\u001b[39;00m\n\u001b[0;32m 17\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mmath\u001b[39;00m\n\u001b[1;32m---> 18\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mcv2\u001b[39;00m\n\u001b[0;32m 20\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mtyping\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m Dict, Any, List\n\u001b[0;32m 21\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mpathlib\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m Path\n", "\u001b[1;31mModuleNotFoundError\u001b[0m: No module named 'cv2'" ] } ], "source": [ "\"\"\"\n", "------------------------------------------------------------------------------------------------------------------------\n", "Creates figures for manuscript\n", "------------------------------------------------------------------------------------------------------------------------\n", "Creates geometry and morphological characteristic plot\n", "Creates plots for comparing comsol vs simgraph\n", "- pressure plot -> xlabel: segment position from the origin ylabel: pressure value\n", "- volummetric flow plot -> xlabel: segment position from the origin ylabel: velocity value\n", "\"\"\"\n", "import os\n", "import string\n", "import vtk\n", "import networkx as nx\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib\n", "import math\n", "import cv2\n", "\n", "from typing import Dict, Any, List\n", "from pathlib import Path\n", "from collections import OrderedDict\n", "from matplotlib import pyplot as plt\n", "from matplotlib.transforms import Bbox\n", "from matplotlib.lines import Line2D\n", "from matplotlib.patches import Rectangle\n", "from matplotlib.offsetbox import OffsetImage, AnnotationBbox\n", "\n", "from image_processing.centerline.read_mesh import mesh_network\n", "from matpancreas.utils import io_utils_py, graph_utils_py\n", "from types import SimpleNamespace\n", "from vedo import * #Lines, Points, show, Plotter, interactive, settings, Picture, screenshot, buildLUT, Axes, Text\n", "from vedo.pyplot import plot\n", "from matpancreas.pancreas.analysis.compare_lines import get_static_data, get_dynamic_data, get_bounds_static\n", "from matpancreas.pancreas.analysis.compare_lines_pscan import get_static_data as get_pscan_static_data,\\\n", " get_dynamic_data as get_pscan_dynamic_data\n", "from matpancreas.settings_model import RESULTS_SIMGRAPH_FILE, \\\n", " RESULTS_COMSOL_FILE, \\\n", " comsol, \\\n", " RESULTS_PLOTS_DIR, \\\n", " RESULTS_PSCAN_FILE, \\\n", " test_case as test, \\\n", " ESPECIES, \\\n", " single, \\\n", " pscan, \\\n", " pressure_scan, \\\n", " RESULTS_SIMGRAPH_DIR, \\\n", " RESULTS_COMSOL_DIR, \\\n", " test_cases, pb, nspecies, SPECIES, RESULTS_GSCAN_FILE\n", "\n", "from matpancreas.utils.graph_utils_py import get_graph, \\\n", " create_graph, \\\n", " get_eucledian_lengths_wrt_origin, \\\n", " draw_graph3d_vedo, \\\n", " plot_histogram, convert_graph_to_df, get_vals_segment\n", "\n", "# vedo settings\n", "from matpancreas.utils.io_utils_py import write_output\n", "\n", "settings.screenshotTransparentBackground = True\n", "settings.screeshotScale = 1\n", "settings.screeshotLargeImage = True\n", "settings.showRendererFrame = False\n", "\n", "# global settings for plots\n", "# plt.rcParams.update({\n", "# 'axes.labelsize': 'large',\n", "# 'axes.labelweight': 'bold'\n", "# })\n", "\n", "# -------------- test label settings ---------------------------------------------------------------------------\n", "left, width = .35, .5\n", "bottom, height = .40, .5\n", "right = left + width\n", "top = bottom + height\n", "fig_labels = dict(zip(test_cases, [\"c\", \"a\", \"b\", \"d\"]))\n", "# --------------------------------------------------------------------------------------------------------------\n", "\n", "\n", "def computeTicks (x, step=5):\n", " \"\"\"\n", " Computes domain with given step encompassing series x\n", " @ params\n", " x - Required - A list-like object of integers or floats\n", " step - Optional - Tick frequency\n", " Reference:\n", " https://stackoverflow.com/questions/12608788/changing-the-tick-frequency-on-x-or-y-axis-in-matplotlib\n", " \"\"\"\n", " xMax, xMin = math.ceil(max(x)), math.floor(min(x))\n", " dMax, dMin = xMax + abs((xMax % step) - step) + (step if (xMax % step != 0) else 0), xMin - abs((xMin % step))\n", " return range(dMin, dMax, step)\n", "\n", "\n", "def plot_histogram(ax, data: List, binwidth: int, density=False, bin=False):\n", " if bin:\n", " ax.hist(\n", " data,\n", " density=density,\n", " color='k',\n", " bins=np.arange(min(data), max(data) + binwidth, binwidth)\n", " )\n", " else:\n", " ax.hist(\n", " data,\n", " density=density,\n", " color='k'\n", " )\n", "\n", "\n", "def draw_weighted_connectivity_matrix(weights: List[List], test_cases):\n", " \"\"\"\n", " This function creates a heatmap\n", " - incidence matrix, M (V x E), is created\n", " - w is weights of edges\n", " - W = diag(w)\n", " - M x W -> column scaled incidence matirx\n", " : param G , graph\n", " :return:\n", " ref: https://github.com/mne-tools/mne-python/issues/5693 (to assign white color to zero)\n", " https://stackoverflow.com/questions/55306803/matplotlib-colorbar-need-to-force-scientific-notation-with-exponent-at-top-of-b\n", " \"\"\"\n", "\n", " fig, axes = plt.subplots(nrows=1, ncols=len(test_cases), figsize=(16, 5))\n", " fig.subplots_adjust(wspace=0.3, hspace=0.3)\n", "\n", " for test_case, weight in zip(test_cases, weights):\n", " i = test_cases.index(test_case)\n", " sheet = test_case.split(\"_pb\")[0]\n", "\n", " output = get_graph(sheet=sheet)\n", " ed_ls = output['ed_ls']\n", " G = create_graph(output=output)\n", " print(G.nodes())\n", " M = nx.incidence_matrix(G, nodelist=sorted(G.nodes), edgelist=ed_ls, oriented=True)\n", " W = np.diag(weight)\n", " MW = M @ W\n", " vmin, vmax = np.min(MW), np.max(MW)\n", "\n", " cbformat = matplotlib.ticker.ScalarFormatter() # create the formatter\n", " cbformat.set_powerlimits((-2, 2)) # set the limits for sci. not.\n", " im = axes[i].matshow(\n", " MW, cmap=plt.cm.bwr, #bwr_r\n", " vmin=-abs(max(vmin, vmax)),\n", " vmax=abs(max(vmin, vmax))\n", " )\n", " axes[i].set_ylabel('Nodes')\n", " axes[i].set_xlabel('Edges')\n", " axes[i].text(\n", " 0.9, 0.9,\n", " '(' + fig_labels[test_case] + ')',\n", " # string.ascii_uppercase[idx],\n", " transform=axes[i].transAxes,\n", " # weight='bold',\n", " size=14\n", " )\n", " plt.colorbar(\n", " im,\n", " ax=axes[i],\n", " fraction=0.057,\n", " pad=0.06,\n", " format=cbformat, #'%.0e'\n", " orientation=\"horizontal\"\n", " ).set_label('qdot ($\\mu$m$^3$/s)', rotation=0)\n", "\n", " return fig\n", "\n", "\n", "def get_plot_data(fpath_s: Path, fpath_c: Path, measure: str, sort_by_distance=False) -> Dict:\n", " \"\"\"\n", " :param sort_by_distance:\n", " :param fpath_s:\n", " :param fpath_c:\n", " :param measure:\n", " :return:\n", " \"\"\"\n", " labels = {'pressure': {'xlabel': 'Node', 'ylabel': \"Pressure (Pa)\"},\n", " 'velocity': {'xlabel': 'Edge', 'ylabel': \"Velocity ($\\mu$m/s)\"},\n", " 'pe_glc_ext': {'xlabel': 'Edge', 'ylabel': \"Axial peclet number\"},\n", " 'pe_lac_ext': {'xlabel': 'Edge', 'ylabel': \"Axial peclet number\"}\n", " }\n", " props = {'node': ['pressure'], 'edge': [\"velocity\"]}\n", "\n", " # xaxis data:\n", " G = get_eucledian_lengths_wrt_origin(fs=RESULTS_SIMGRAPH_FILE)\n", "\n", " if measure == \"pressure\":\n", " d = nx.get_node_attributes(G, 'ed_from_source') # eucledian distance from source node\n", " elif measure == \"velocity\":\n", " d = nx.get_edge_attributes(G, 'ed_from_source')\n", "\n", " # yaxis data\n", " if single:\n", " scan = get_static_data(fs=fpath_s, fc=fpath_c, interpolate=False)\n", " elif pscan:\n", " scan = get_pscan_static_data(f=fpath_s, interpolate=False)\n", "\n", " y = {}\n", " for key in scan.keys():\n", " data = scan[key]\n", " node_attr, edge_attr = data['node_attr'], data['edge_attr']\n", "\n", " if measure in props['node']:\n", " m = OrderedDict(zip(sorted(G.nodes), node_attr[measure]))\n", " elif measure in props['edge']:\n", " # m = OrderedDict(zip(sorted(G.edges), edge_attr[measure]))\n", " m = OrderedDict(zip(range(1, len(G.edges())), edge_attr[measure]))\n", " if sort_by_distance:\n", " temp = dict((v, k) for k, v in d.items()) # swaps and xdata will be position\n", " temp = OrderedDict(sorted(temp.items(), key=lambda kv: kv[0])) # sorted (pos, node) / (pos, edge)\n", " sorted_keys = list(temp.values())\n", " m = OrderedDict(sorted(m.items(), key=lambda pair: sorted_keys.index(pair[0]))) # sorted(nodes, pressure)\n", "\n", " x = m.keys() # xdata\n", " y[key] = m.values() # ydata\n", "\n", " return {'x': x, 'y': y, 'xlabel': labels[measure]['xlabel'], 'ylabel': labels[measure]['ylabel']}\n", "\n", "\n", "def plot_static_data(scan: List):\n", " \"\"\"\"\n", " plots comsol vs simgraph data\n", " :param data: x and multiple y values\n", " \"\"\"\n", " if single:\n", " marker = {'simgraph': 's', 'comsol': '^'}\n", " color = {'simgraph': 'k', 'comsol': 'grey'}\n", " linewidth = {'simgraph': 0, 'comsol': 0.5}\n", " elif pscan:\n", " symbols = ['p', '^', 's', 'o', 'v', 'd', '+', 'x', \"_\", \"1\"]\n", " colors = ['grey', 'grey', 'grey', 'grey', 'grey', 'grey', 'grey', 'grey', 'grey', 'grey']\n", " widths = [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5]\n", " marker = OrderedDict(zip(pressure_scan, symbols))\n", " color = OrderedDict(zip(pressure_scan, colors))\n", " linewidth = OrderedDict(zip(pressure_scan, widths))\n", "\n", " # fig = plt.figure()\n", " fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(16, 6))\n", " fig.subplots_adjust(wspace=0.3, hspace=0.3)\n", " axes = (ax1, ax2)\n", " for ax, data in zip(axes, scan):\n", " x = [str(i) for i in data['x']] # data['x']\n", " y = data['y']\n", "\n", " for key in data['y'].keys():\n", " if single:\n", " label = key\n", " elif pscan:\n", " label = f'$\\Delta$P {key} (Pa)'\n", " ax.plot(\n", " x, y[key],\n", " color=color[key],\n", " marker=marker[key],\n", " label=label,\n", " linestyle='dashed',\n", " linewidth=linewidth[key]\n", " )\n", "\n", " ax.set_xlabel(data['xlabel'], fontsize=12)\n", " ax.set_ylabel(data['ylabel'], fontsize=12)\n", " ax.text(\n", " # 0.95, 0.01, # (right, bottom)\n", " 0.9, 0.85,\n", " '(' + chr(scan.index(data)+97) + ')',\n", " # string.ascii_uppercase[idx],\n", " transform=ax.transAxes,\n", " # weight='bold',\n", " size=15\n", " )\n", " xdata = [int(i) for i in x]\n", " # ax.set_xticks(np.arange(0, len(x)+1, 4)) #rotation=90\n", " # ax.set_yticks()\n", " leg = ax.legend(\n", " loc='best',\n", " ncol=1,\n", " fontsize=12\n", " )\n", " leg.get_frame().set_alpha(0.5)\n", " plt.setp(ax.get_xticklabels(), rotation=90)\n", " plt.show()\n", "\n", " return fig\n", "\n", "\n", "def plot_peclet_distribution(fpath: Path, pdata: List = pressure_scan, greyscale=False):\n", " \"\"\"\n", " plots peclet distributions\n", " reference:\n", " ---------\n", " https://stackoverflow.com/questions/43872450/matplotlib-histogram-with-multiple-legend-entries\n", " https://stackoverflow.com/questions/20118258/matplotlib-coloring-line-plots-by-iteration-dependent-gray-scale\n", " \"\"\"\n", " scan = get_pscan_static_data(f=fpath, interpolate=False)\n", "\n", " fig, axes = plt.subplots(nrows=1, ncols=1) #, figsize=(16, 5))\n", " fig.subplots_adjust(wspace=0.3, hspace=0.3)\n", " bins = 10\n", " if greyscale:\n", " colors = [(t/2.0, t/2.0, t/2.0, 0.5) for t in np.arange(0., 2., 0.4)]\n", " fc = dict(zip(pdata, colors))\n", " linestyle = dict(zip(pdata, 'solid'))\n", "\n", " else:\n", " # fc = dict(zip(pdata, [(0, 0, 1, 0.5), (1, 0, 0, 0.5), (0, 0, 0, 0.5)])) #, (1, 0, 0.5, 0.25), (0.5, 0, 0.5, 0.25)]))\n", " fc = dict(zip(pdata, [(0, 0, 0, 0.5), (0, 0, 0, 0.5), (0, 0, 0, 0.5)]))\n", " linestyle = dict(zip(pdata, ['dashed', 'solid', 'dotted']))\n", "\n", " for i, key in enumerate(pdata):\n", " data = scan[key]\n", " node_attr, edge_attr = data['node_attr'], data['edge_attr']\n", " x = edge_attr['pe_glc_ext']\n", " plt.hist(\n", " x,\n", " bins=bins,\n", " stacked=True,\n", " density=True,\n", " fc=fc[key],\n", " # fill=False,\n", " histtype='step',\n", " linestyle=linestyle[key],\n", " color='k', # fc[key],\n", " linewidth=2\n", " )\n", "\n", " plt.ylabel(\"Density\", fontsize=10)\n", " plt.xlabel(\"Axial peclet number\", fontsize=10)\n", " plt.xticks(fontsize=10)\n", " plt.yticks(fontsize=10)\n", " # plt.xscale('log')\n", "\n", " # create legend\n", " # handles = [Rectangle((0, 0), 1, 1, color=c, ec=\"k\") for c in fc.values()]\n", " handles = [Line2D([0], [0], color='k', linewidth=3, linestyle=ls) for ls in ['dashed', 'solid']]\n", " labels = [f\"$\\Delta$P = {p} (Pa)\" for p in pdata]\n", " plt.legend(handles, labels)\n", " plt.show()\n", " return fig\n", "\n", "\n", "def plot_conc_pscan(fpath: Path):\n", " \"\"\" plots conc vs nodes \"\"\"\n", "\n", " # symbols = ['+', '^', 's', 'p', 'v', 'd', 'o', 'x', \"_\", \"1\"]\n", " # widths = [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5]\n", " # marker = OrderedDict(zip(pressure_scan, symbols))\n", " # linewidth = OrderedDict(zip(pressure_scan, widths))\n", "\n", " fig = plt.figure()\n", " time, scan = get_pscan_dynamic_data(f=fpath, resample=False)\n", " print(scan.keys())\n", " # nodes to plot (inlet, center, outlet) (34, 4 , 39)\n", " nodes = [27, 11, 38]\n", " color = OrderedDict(zip(nodes, ['r', 'g', 'b']))\n", " linestyle = OrderedDict(zip(scan.keys(), ['dashed', 'solid']))\n", "\n", " data = {}\n", " for idx, s in enumerate(scan):\n", " species = scan[s][\"glc_ext\"]\n", " data[s] = pd.DataFrame(species)\n", " for key in data.keys():\n", " print(key)\n", " for n in nodes:\n", " plt.plot(\n", " time, data[key][n],\n", " color=color[n],\n", " # marker=marker[key],\n", " markersize=5,\n", " label=key,\n", " linestyle=linestyle[key],\n", " # linewidth=linewidth[key]\n", " )\n", " plt.xlabel('Time (s)', fontsize=10)\n", " plt.ylabel('Concentration (mM)', fontsize=10)\n", " plt.show()\n", " return fig\n", "\n", "\n", "def load_G_H():\n", " G = create_graph(attributes=True)\n", " input = get_graph(sheet=\"test11\")\n", " H = create_graph(output=input, derived_attr=True)\n", " return G, input, H\n", "\n", "\n", "class ManuscriptFigures:\n", " def figures(self):\n", " # self.draw_graph3d_vedo()\n", " self.plot_3dhistogram()\n", " # self.plot_morphological_params()\n", " # self.plot_qdot()\n", " # self.Fig1() # static properties single run comparison of comsol vs. simgraph\n", " # self.Fig2() # static properties pscan run\n", " self.Fig3() # pe distribution\n", " # self.Fig4() # plot of pressure scan at diferent pe\n", " # self.Fig5()\n", " # self.Fig6() # conc plot of simgraph vs comsol\n", " # self.Fig7()\n", " # self.Fig8()\n", " # self.Fig9() #line plot of pressure and velocity\n", " # self.Fig10() # methods figure 2\n", " # self.Fig11() # volume interpolation plot for glc, lac\n", " # self.Fig12() # probe line plot for glc and lac exchange\n", " # self.Fig13() # plots steady-state/rise times of simgraph vs. comsol\n", " # self.Fig14() # plots steady-state/rise times of tumor design 1 vs. design 2\n", " # self.Fig15() # plots bv cell glucose and lactate time-varying concentration in 3D vasculature\n", " # self.Fig16() # glucose uptake and latate release flux\n", " # self.Fig17()\n", " # self.Fig18()\n", " # self.Fig19()\n", "\n", " def plot_3dhistogram(self):\n", " \"\"\"\n", " :param data: Attribute data to plot\n", " :param label: xaxis label\n", " :param binwidth:\n", " :param density:\n", " :return:\n", " reference:\n", " https://stackoverflow.com/questions/6871201/plot-two-histograms-on-single-chart-with-matplotlib\n", " example:\n", " import random\n", " import numpy\n", " from matplotlib import pyplot\n", " x = [random.gauss(3, 1) for _ in range(400)]\n", " y = [random.gauss(4, 2) for _ in range(400)]\n", " bins = numpy.linspace(-10, 10, 100)\n", " pyplot.hist(x, bins, alpha=0.5, label='x')\n", " pyplot.hist(y, bins, alpha=0.5, label='y')\n", " pyplot.legend(loc='upper right')\n", " pyplot.show()\n", " \"\"\"\n", " fc = dict(zip(test_cases, ['b', 'r', 'k', 'g']))\n", " linestyle = dict(zip(test_cases, ['solid', (0, (5, 10)), 'dotted', 'solid']))\n", "\n", " fig, axes = plt.subplots(nrows=2, ncols=1) # , figsize=(16, 5))\n", " # fig.subplots_adjust(wspace=0.3, hspace=0.3)\n", "\n", " data = {}\n", " for sheet in test_cases:\n", " data[sheet] = get_graph(sheet=sheet)\n", "\n", " param = {'d': {'bw': 5, 'label': 'Diameter', 'figlabel': 'e'}, 'l': {'bw': 10, 'label': 'Length', 'figlabel': 'f'}}\n", " for i, p in enumerate(param):\n", " for idx, key in enumerate(data):\n", " x = data[key][p]\n", " ax = axes[i]\n", " ax.hist(\n", " x,\n", " # bins=10,\n", " stacked=False,\n", " density=True,\n", " histtype='step',\n", " linestyle=linestyle[key],\n", " color=fc[key],\n", " linewidth=2\n", " )\n", " # create legend\n", " handles = [Line2D([0], [0], color=c, linewidth=3, linestyle=ls) for c, ls in zip(fc.values(), linestyle.values())]\n", " labels = [\"islet\", \"tumor design 1\", \"tumor design 2\", \"mesentery\"]\n", " ax.legend(handles, labels, prop={\"size\":10})\n", " ax.set_ylabel(\"Density\", fontsize=12)\n", " ax.set_xlabel(f\"{param[p]['label']} ($\\mu$m)\", fontsize=12)\n", " ax.set_xscale('log')\n", " ax.text(\n", " 0.008, top,\n", " # '(' + chr(i + 97) + ')',\n", " '(' + param[p]['figlabel'] + ')',\n", " transform=ax.transAxes,\n", " size=12\n", " )\n", " plt.show()\n", " f = os.path.join(RESULTS_PLOTS_DIR, \"morphology.svg\")\n", " fig.savefig(f, transparent=True, dpi=600, bbox_inches=\"tight\")\n", "\n", " def draw_graph3d_vedo(self):\n", " \"\"\"\n", " :return:\n", " \"\"\"\n", " # model , label : font size\n", "\n", " font_size = {\"test9_default\": 35, \"test10_default\": 35, \"test11_default\": 5, \"test24_default\": 125}\n", " settings.windowSplittingPosition = 0.4\n", " vp = Plotter(shape='3/1', interactive=True, sharecam=False)\n", "\n", " for i, sheet in enumerate(test_cases):\n", " input = get_graph(sheet=sheet)\n", " source = int(input['source'])\n", " target = int(input['target'])\n", " t = {source: 'In', target: 'Out'} # terminals inlet exit\n", " c = {source: (0, 0, 0), target: (0, 0, 0)}\n", "\n", " G = create_graph(output=input)\n", " nodes = sorted(G.nodes())\n", "\n", " # create labels\n", " l = []\n", " color = []\n", " for n in nodes:\n", " if n in t.keys():\n", " l.append(t[n])\n", " color.append(c[n])\n", " else:\n", " l.append('')\n", " color.append((128, 128, 128))\n", "\n", " # create vtk polydata\n", " pos = dict(zip(input['nodes'], input['xyz']))\n", " lines = [(pos[t], pos[h]) for t, h in input['ed_ls']]\n", " pts = Points(input['xyz'], r=8, c=color).rotateX(0)\n", " edg = Lines(lines).lw(5).rotateX(0).c(\"gray\")\n", "\n", " scale = font_size[sheet]\n", " labs = pts.labels(l, scale=scale, font='VictorMono', c='k').addPos(0.1, 1, 0.5) #FIXME: revert to (0.05,0,0.5)\n", " f = os.path.join(RESULTS_PLOTS_DIR, f'{sheet}.png')\n", "\n", " if True:\n", " plt = Plotter(N=1) # , size=(800, 600)\n", " if sheet == \"test11_default\":\n", " axes = Axes(\n", " edg, xtitle='y (\\mum)', ytitle='z (\\mum)', ztitle='x (\\mum)', xyGrid=True, zxGrid=True, yzGrid=True,\n", " xLabelSize=0.022, yLabelSize=0.022, zLabelSize=0.022, xTitleSize=0.03, yTitleSize=0.03,\n", " zTitleSize=0.03\n", " )\n", " elif sheet == \"test9_default\":\n", " axes = Axes(\n", " edg, xyGrid=True, zxGrid=True, yzGrid=True, xtitle='x (\\mum)', ytitle='y (\\mum)', ztitle='z (\\mum)',\n", " xLabelSize=0.022, yLabelSize=0.022, zLabelSize=0.022, xTitleSize=0.03, yTitleSize=0.03,\n", " zTitleSize=0.03, xTitleOffset=-0.01, yTitleOffset=-0.05,\n", " xTitlePosition=1.05, yTitlePosition=1.05,\n", " )\n", " elif sheet == \"test10_default\":\n", " axes = Axes(\n", " edg, xyGrid=True, zxGrid=True, yzGrid=True, xtitle='x (\\mum)', ytitle='y (\\mum)',\n", " ztitle='z (\\mum)',\n", " xLabelSize=0.022, yLabelSize=0.022, zLabelSize=0.022, xTitleSize=0.03, yTitleSize=0.03,\n", " zTitleSize=0.03, xTitleOffset=-0.01, yTitleOffset=-0.05,\n", " xTitlePosition=1.05, yTitlePosition=1.1,\n", " )\n", " elif sheet == \"test24_default\":\n", " axes = Axes(\n", " edg,\n", " xyGrid=True, zxGrid=True, yzGrid=True, xTitlePosition=0.15, yTitlePosition=1.02,\n", " xtitle='x (\\mum)', ytitle='y (\\mum)', ztitle='z (\\mum)',\n", " yLabelRotation=-90, xLabelRotation=-90, xTitleRotation=0, yTitleRotation=-180, xTitleOffset=-0.06, yTitleOffset=0.001,\n", " xTitleSize=0.022, yTitleSize=0.022, zTitleSize=0.022\n", " )\n", " plt.show(\n", " pts,\n", " edg,\n", " labs,\n", " axes,\n", " bg2='w',\n", " interactive=True,\n", " # zoom=1.2,\n", " at=0,\n", " # roll=90\n", " ).screenshot(f)\n", " interactive()\n", " if False:\n", " axes = Axes(edg, xyGrid=True, zxGrid=True, yzGrid=True) # zLabelRotation=-45)\n", " vp.show(\n", " pts,\n", " edg,\n", " labs,\n", " axes,\n", " bg2='w',\n", " at=i\n", " )\n", "\n", " def plot_morphological_params(self):\n", " \"\"\"\n", " generates histogram plots of length and diameter distribution\n", " :return:\n", " Reference:\n", " ----------\n", " enumerate subplot: https://stackoverflow.com/questions/22508590/enumerate-plots-in-matplotlib-figure\n", " enumerate subplot: https://stackoverflow.com/a/25544329/8281509\n", " position text label: https://stackoverflow.com/questions/8482588/putting-text-in-top-left-corner-of-matplotlib-plot\n", " \"\"\"\n", "\n", " data = {}\n", " for sheet in test_cases:\n", " data[sheet] = get_graph(sheet=sheet)\n", "\n", " fig, ((ax1, ax2, ax3, ax4), (ax5, ax6, ax7, ax8)) = plt.subplots(nrows=2, ncols=4, figsize=(16, 12))\n", " fig.subplots_adjust(wspace=0.3, hspace=0.3)\n", " axes = (ax1, ax2, ax3, ax4, ax5, ax6, ax7, ax8)\n", " param = {'d': {'bw': 5, 'label': 'Diameter'}, 'l': {'bw': 10, 'label': 'Length'}}\n", "\n", " k = 0\n", "\n", " for p in param.keys():\n", " for idx, test_case in enumerate(data):\n", " bw = param[p]['bw']\n", " label = param[p]['label']\n", " ax = axes[k]\n", " if idx == 0:\n", " ax.set_ylabel('Frequency', fontsize=12)\n", "\n", " ax.set_xlabel(f\"{label} ($\\mu$m)\", fontsize=12)\n", " ax.text(\n", " right, top,\n", " '('+fig_labels[test_case]+')',\n", " # string.ascii_uppercase[idx],\n", " transform=ax.transAxes,\n", " # weight='bold',\n", " size=20\n", " )\n", " plot_histogram(ax, data=data[test_case][p], binwidth=bw)\n", " k = k + 1\n", " plt.show()\n", " f = os.path.join(RESULTS_PLOTS_DIR, \"morphology.svg\")\n", " fig.savefig(f, transparent=True, dpi=600, bbox_inches=\"tight\")\n", "\n", " def plot_qdot(self):\n", " weights = []\n", " for test_case in test_cases:\n", " fpath = os.path.join(RESULTS_SIMGRAPH_DIR, f\"{test_case}_bc1_v1_c0.mat\")\n", " scan = get_static_data(fs=fpath, interpolate=False)\n", " for idx, key in enumerate(scan):\n", " data = scan[key]\n", " node_attr, edge_attr = data['node_attr'], data['edge_attr']\n", " weights.append(edge_attr['qdot'])\n", "\n", " fig = draw_weighted_connectivity_matrix(weights=weights, test_cases=test_cases)\n", "\n", " if pb:\n", " f = \"qdot_pb\"\n", " else:\n", " f = \"qdot\"\n", "\n", " fpath = os.path.join(RESULTS_PLOTS_DIR, f'{f}.svg')\n", " fig.savefig(fpath, transparent=True, bbox_inches='tight', dpi=600)\n", " plt.show()\n", "\n", " def Fig1(self):\n", " \"\"\" static properties for single run \"\"\"\n", " for test_case in test_cases:\n", " fs = os.path.join(RESULTS_SIMGRAPH_DIR, f\"{test_case}_bc1_v1_c0.mat\")\n", " fc = os.path.join(RESULTS_COMSOL_DIR, f\"{test_case}_bc1_v1_c0.mat\")\n", " print(fs)\n", " # pressure\n", " pressure_data = get_plot_data(\n", " fpath_s=fs,\n", " fpath_c=fc,\n", " measure=\"pressure\",\n", " sort_by_distance=False\n", " )\n", "\n", " # velocity\n", " velocity_data = get_plot_data(\n", " fpath_s=fs,\n", " fpath_c=fc,\n", " measure=\"velocity\",\n", " sort_by_distance=False\n", " )\n", "\n", " fig = plot_static_data(scan=[pressure_data, velocity_data])\n", " f = os.path.join(RESULTS_PLOTS_DIR, f\"{test_case}_pressure_velocity_scatter.svg\")\n", " fig.savefig(f, transparent=False, dpi=600, bbox_inches=\"tight\")\n", "\n", " def Fig2(self):\n", " \"\"\" static properties for pscan run \"\"\"\n", " print(RESULTS_PSCAN_FILE)\n", " pressure_data = get_plot_data(\n", " fpath_s=RESULTS_PSCAN_FILE,\n", " fpath_c=None,\n", " measure=\"pressure\",\n", " sort_by_distance=False\n", " )\n", "\n", " # velocity\n", " velocity_data = get_plot_data(\n", " fpath_s=RESULTS_PSCAN_FILE,\n", " fpath_c=None,\n", " measure=\"velocity\",\n", " sort_by_distance=False\n", " )\n", "\n", " fig = plot_static_data(scan=[pressure_data, velocity_data])\n", " f = os.path.join(RESULTS_PLOTS_DIR, f\"{test}_pressure_velocity_scatter_pscan.svg\")\n", " fig.savefig(f, transparent=True, dpi=600, bbox_inches=\"tight\")\n", "\n", " def Fig3(self):\n", " \"\"\" plot 2D histogram of axial peclet numbers for species = glc_ext \"\"\"\n", " fig = plot_peclet_distribution(\n", " fpath=RESULTS_PSCAN_FILE,\n", " pdata=[20, 200]\n", " )\n", " f = os.path.join(RESULTS_PLOTS_DIR, f\"{test}_pe_pscan.svg\")\n", " fig.savefig(\n", " f,\n", " transparent=True,\n", " dpi=600,\n", " bbox_inches=\"tight\"\n", " )\n", "\n", " def Fig4(self):\n", " \"\"\" plot conc of pscan at different peclet numbers and different nodal ppositions of the vasculature\"\"\"\n", "\n", " # symbols = ['+', '^', 's', 'p', 'v', 'd', 'o', 'x', \"_\", \"1\"]\n", " # widths = [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5]\n", " # marker = OrderedDict(zip(pressure_scan, symbols))\n", " # linewidth = OrderedDict(zip(pressure_scan, widths))\n", " # gimg = os.path.join(RESULTS_PLOTS_DIR, f'{test}.png')\n", "\n", " # ax = plt.subplot(111)\n", " fig, ax = plt.subplots()\n", " time, scan = get_pscan_dynamic_data(f=RESULTS_PSCAN_FILE, resample=False)\n", " # nodes to plot (inlet, center, outlet) (34, 4 , 39)\n", " nodes = [27, 5, 38] #5>11\n", " color = OrderedDict(zip(nodes, ['r', 'g', 'b']))\n", " print(scan.keys())\n", " linestyle = OrderedDict(zip(pressure_scan, ['dashed', 'solid']))\n", "\n", " # --------------------------------------------------------------------------------------------------------------\n", " # main plot\n", " # --------------------------------------------------------------------------------------------------------------\n", "\n", " data = {}\n", " for idx, s in enumerate(scan):\n", " species = scan[s][\"glc_ext\"]\n", " data[s] = pd.DataFrame(species)\n", " for key in data.keys():\n", " if key in pressure_scan:\n", " for n in nodes:\n", " ax.plot(\n", " time, data[key][n],\n", " color=color[n],\n", " # marker='*',#marker[key],\n", " # markersize=5,\n", " label=key,\n", " linestyle=linestyle[key],\n", " # linewidth=linewidth[key]\n", " )\n", " ax.set_xlabel('Time (s)', fontsize=12)\n", " ax.set_ylabel('Concentration (mM)', fontsize=12)\n", "\n", " # add inset\n", " # arr_img = plt.imread(gimg)\n", " # im = OffsetImage(arr_img, zoom=72. / 120)\n", " # ab = AnnotationBbox(im, (0.75, 0.25), xycoords='axes fraction', frameon=False)\n", " # ax.add_artist(ab)\n", "\n", " # --------------------------------------------------------------------------------------------------------------\n", " # inset plot\n", " # --------------------------------------------------------------------------------------------------------------\n", "\n", " rect = [.45, 0.24, .5, .4]\n", " bbox = Bbox.from_bounds(*rect)\n", " inax = fig.add_axes(bbox)\n", " inax.axis('on')\n", " inax.set_facecolor('none')\n", " # inax.grid('on', color='k')\n", "\n", " scan = get_pscan_static_data(f=RESULTS_PSCAN_FILE, interpolate=False)\n", " bins = 10\n", " fc = dict(zip(pressure_scan, [(0, 0, 0, 0.5), (0, 0, 0, 0.5), (0, 0, 0, 0.5)]))\n", " linestyle = dict(zip(pressure_scan, ['dashed', 'solid', 'dotted']))\n", "\n", " for i, key in enumerate(pressure_scan):\n", " data = scan[key]\n", " node_attr, edge_attr = data['node_attr'], data['edge_attr']\n", " x = edge_attr['pe_glc_ext']\n", " inax.hist(\n", " x,\n", " bins=bins,\n", " stacked=True,\n", " density=True,\n", " fc=fc[key],\n", " # fill=False,\n", " histtype='step',\n", " linestyle=linestyle[key],\n", " color='k', # fc[key],\n", " linewidth=2\n", " )\n", "\n", " inax.set_ylabel(\"Density\", fontsize=10)\n", " inax.set_xlabel(\"Axial peclet number\", fontsize=10)\n", "\n", " fig.show()\n", " f = os.path.join(RESULTS_PLOTS_DIR, f\"{test}_conc_pscan.svg\")\n", " ax.figure.savefig(\n", " f,\n", " transparent=True,\n", " dpi=600,\n", " bbox_inches=\"tight\"\n", " )\n", " # f = os.path.join(RESULTS_PLOTS_DIR, f\"{test}_conc_pscan.svg\")\n", " # ax.figure.savefig(f, transparent=True, dpi=600, bbox_inches=\"tight\")\n", " # plt.show()\n", "\n", " def Fig5(self):\n", " \"\"\"plot fluxes of pscan\"\"\"\n", " fig = plt.figure()\n", " data = io_utils_py._get_file(RESULTS_PSCAN_FILE, spreadsheet=False)\n", " glcim, lacex, delp = [], [], []\n", " for s in data['variable']:\n", " glcim.append(s.glcim_net)\n", " lacex.append(s.lacex_net)\n", " delp.append(s.delP)\n", "\n", " plt.plot(\n", " delp,\n", " glcim,\n", " color='k'\n", " )\n", " plt.xlabel('Time (s)')\n", " plt.ylabel('Concentration (mM)')\n", " plt.show()\n", " f = os.path.join(RESULTS_PLOTS_DIR, f\"{test}_flux_pscan.svg\")\n", " fig.savefig(f, transparent=True, dpi=600, bbox_inches=\"tight\")\n", "\n", " def Fig6(self):\n", " \"\"\"\n", " plot conc for comparing simgraph vs. comsol\n", " inset figure is loaded from a file\n", " references:\n", " ----------\n", " https://stackoverflow.com/questions/66540026/insert-an-svg-image-in-matplotlib-figure\n", " https://stackoverflow.com/questions/12161559/plot-svg-within-matplotlib-figure-embedded-in-wxpython\n", " \"\"\"\n", " gimg = os.path.join(RESULTS_PLOTS_DIR, f'{test}_inset.png')\n", " fs = os.path.join(RESULTS_SIMGRAPH_DIR, f\"{test}_bc1_v1_c0.mat\")\n", " fc = os.path.join(RESULTS_COMSOL_DIR, f\"{test}_bc1_v1_c0.mat\")\n", " time, scan = get_dynamic_data(fs=fs, fc=fc)\n", "\n", " # --------------------------------------------------------------------------------------------------------------\n", " # settings and labels\n", " # --------------------------------------------------------------------------------------------------------------\n", " linestyle = {'simgraph': 'solid', 'comsol': 'dashed'}\n", " if test == \"test11_default\":\n", " nodes = [27, 5, 38]\n", " colors = ['r', 'g', 'b']\n", " elif test == \"test24_default\":\n", " # nodes = [79, 19, 531, 13, 606, 577, 151, 311] # neighbors of 13 (53,606,577)\n", " # colors = ['r', 'c', 'y', 'k', 'brown', 'g', 'm', 'b']\n", " nodes = [79, 19, 13, 151, 311] # highlight nodes\n", " colors = ['r', 'c', 'g', 'm', 'b']\n", " color = OrderedDict(zip(nodes, colors))\n", "\n", " # --------------------------------------------------------------------------------------------------------------\n", " # main plot\n", " # --------------------------------------------------------------------------------------------------------------\n", " ax = plt.subplot(111)\n", " for key in scan.keys():\n", " species = pd.DataFrame(scan[key][\"glc_ext\"])\n", " print(species)\n", " for n in nodes:\n", " ax.plot(\n", " time, species[n-1],\n", " color=color[n],\n", " markersize=5,\n", " label=key,\n", " linestyle=linestyle[key],\n", " )\n", " ax.set_xlabel('Time (s)', fontsize=10)\n", " ax.set_ylabel('Concentration (mM)', fontsize=10)\n", "\n", " # --------------------------------------------------------------------------------------------------------------\n", " # inset plot\n", " # --------------------------------------------------------------------------------------------------------------\n", " if True:\n", " draw_graph3d_vedo(\n", " point_r=10,\n", " lw=2,\n", " points=False,\n", " sphere=True,\n", " highlight=True,\n", " line_alpha=1,\n", " fout=f\"{test}_inset.png\",\n", " default_axes=False,\n", " interact=True,\n", " sphere_r=30, #test24\n", " label_scale=100, #test24\n", " label_xoffset=80 #test24\n", " )\n", "\n", " # G = get_eucledian_lengths_wrt_origin(fs=RESULTS_SIMGRAPH_FILE, origin=34)\n", " # d = nx.get_node_attributes(G, \"ed_from_source\")\n", " if False:\n", " arr_img = plt.imread(gimg)\n", " im = OffsetImage(arr_img, zoom=0.45)\n", " ab = AnnotationBbox(im, (1, 0), xycoords='axes fraction', box_alignment=(1.1, -0.1), frameon=False)\n", " ax.add_artist(ab)\n", " plt.show()\n", "\n", " f = os.path.join(RESULTS_PLOTS_DIR, f\"{test}_conc_compare.svg\")\n", " ax.figure.savefig(\n", " f,\n", " transparent=True,\n", " dpi=600,\n", " bbox_inches=\"tight\"\n", " )\n", "\n", " def Fig7(self):\n", " \"\"\" plot conc fo pscan\"\"\"\n", " fs = os.path.join(RESULTS_SIMGRAPH_DIR, f\"{test}_pb_bc1_v1_c0.mat\")\n", " fc = os.path.join(RESULTS_COMSOL_DIR, f\"{test}_pb_bc1_v1_c0.mat\")\n", " time, scan = get_dynamic_data(fs=fs, fc=fc)\n", "\n", " # --------------------------------------------------------------------------------------------------------------\n", " # settings and labels\n", " # --------------------------------------------------------------------------------------------------------------\n", " linestyle = {'simgraph': 'solid', 'comsol': 'dashed'}\n", " nodes = [27, 11, 38]\n", " colors = ['r', 'g', 'b']\n", " color = OrderedDict(zip(nodes, colors))\n", "\n", " # --------------------------------------------------------------------------------------------------------------\n", " # inset plot\n", " # --------------------------------------------------------------------------------------------------------------\n", "\n", " vplt = draw_graph3d_vedo(\n", " point_r=10,\n", " lw=2,\n", " points=False,\n", " sphere=True,\n", " highlight=True,\n", " line_alpha=1,\n", " rotatex=0,\n", " )\n", " np_pic = Picture(screenshot(returnNumpy=True, scale=1))\n", " vplt.close()\n", "\n", " # --------------------------------------------------------------------------------------------------------------\n", " # main plot\n", " # --------------------------------------------------------------------------------------------------------------\n", " ax = plt.subplot(111)\n", " for key in scan.keys():\n", " species = pd.DataFrame(scan[key][\"glc_ext\"])\n", " for n in nodes:\n", " ax.plot(\n", " time, species[n],\n", " color=color[n],\n", " markersize=5,\n", " label=key,\n", " linestyle=linestyle[key],\n", " )\n", " ax.set_xlabel('Time (s)', fontsize=10)\n", " ax.set_ylabel('Concentration (mM)', fontsize=10)\n", "\n", " f = os.path.join(RESULTS_PLOTS_DIR, f\"{test}_conc_compare.svg\")\n", " ax.figure.savefig(\n", " f,\n", " transparent=True,\n", " dpi=600,\n", " bbox_inches=\"tight\"\n", " )\n", "\n", " # ------------------------------------------------------------\n", " # Build plot: exponential\n", " # ------------------------------------------------------------\n", " x = np.arange(0, 4, 0.1)\n", " y = 3 * np.exp(-x)\n", "\n", " plt1 = plot(\n", " x, y,\n", " xtitle='time in seconds',\n", " ytitle='some function [a.u.]',\n", " )\n", "\n", " np_pic.scale(0.025).x(2).y(1.)\n", "\n", " show(plt1, np_pic, size=(800, 600), zoom=1.2)\n", "\n", " # ------------------------------------------------------------------------------------------------------------------\n", " # figure8\n", " # ------------------------------------------------------------------------------------------------------------------\n", " def Fig8(self):\n", " \"\"\" plot conc fo pscan using 3d plot generated in matplotlib\n", " reference:\n", " \"\"\"\n", " fs = os.path.join(RESULTS_SIMGRAPH_DIR, f\"{test}_bc1_v1_c0.mat\") # FIXME:revert to _pb for test11\n", " fc = os.path.join(RESULTS_COMSOL_DIR, f\"{test}_bc1_v1_c0.mat\")\n", " time, scan = get_dynamic_data(fs=fs, fc=fc)\n", " print(scan.keys())\n", " # --------------------------------------------------------------------------------------------------------------\n", " # settings and labels\n", " # --------------------------------------------------------------------------------------------------------------\n", " G = create_graph()\n", " terminals = [335, 352] # FIXME: revert to [34, 39] for test11 # source & target\n", " hnodes = [27, 11, 38] # FIXME: revert to [40, 89, 123, 161, 312, 352] for test24 # highlight nodes\n", " terminal_labels = dict(zip(terminals, ['In', 'Out']))\n", " hnode_color = OrderedDict(zip(hnodes + terminals, ['r', 'g', 'b', 'y', 'c', 'm', 'k', 'k'])) # FIXME: revert to ['r', 'g', 'b', 'k', 'k'] for test11\n", " hnode_alpha = dict.fromkeys(hnodes + terminals, 1)\n", " label = dict.fromkeys(G.nodes(), '')\n", " color = dict.fromkeys(G.nodes(), \"grey\")\n", " alpha = dict.fromkeys(G.nodes(), 0.3)\n", " label.update(terminal_labels)\n", " color.update(hnode_color)\n", " alpha.update(hnode_alpha)\n", " linestyle = {'simgraph': 'solid', 'comsol': 'dashed'}\n", "\n", " # --------------------------------------------------------------------------------------------------------------\n", " # main plot\n", " # --------------------------------------------------------------------------------------------------------------\n", " fig, ax = plt.subplots()\n", " for key in scan.keys():\n", " species = pd.DataFrame(scan[key][\"glc_ext\"])\n", " for n in hnodes:\n", " ax.plot(\n", " time, species[n],\n", " color=color[n],\n", " marker='*',\n", " markersize=5,\n", " label=key,\n", " linestyle=linestyle[key],\n", " )\n", " ax.set_xlabel('Time (s)', fontsize=10)\n", " ax.set_ylabel('Concentration (mM)', fontsize=10)\n", "\n", " # --------------------------------------------------------------------------------------------------------------\n", " # inset plot\n", " # --------------------------------------------------------------------------------------------------------------\n", "\n", " # rect = [.3, 0.1, .8, .7]\n", " rect = [.5, 0.24, .4, .3]\n", " bbox = Bbox.from_bounds(*rect)\n", " inax = fig.add_axes(bbox) #, projection='3d')\n", " inax.axis('on')\n", " # inax.view_init(azim=90) #, elev=0)\n", " inax.set_facecolor('none')\n", " # inax.grid('on', color='k')\n", "\n", " seen = set()\n", " for i, j in G.edges():\n", " c1 = nx.get_node_attributes(G, 'pos')[i]\n", " c2 = nx.get_node_attributes(G, 'pos')[j]\n", " x = np.stack((c1, c2))\n", " inax.plot(*x.T, color='grey')\n", "\n", " if i not in seen:\n", " inax.scatter(*x[0], color=color[i], alpha=alpha[i], s=15)\n", " inax.text(*x[0], label[i], size=10, color='k')\n", " seen.add(i)\n", " if j not in seen:\n", " inax.scatter(*x[1], color=color[j], alpha=alpha[j], s=15)\n", " inax.text(*x[1], label[j], size=10, color='k')\n", " seen.add(j)\n", "\n", " inax.set_xlabel('x')\n", " inax.set_ylabel('y')\n", " inax.set_zlabel('z')\n", " inax.set_xticklabels([])\n", " inax.set_yticklabels([])\n", " inax.set_zticklabels([])\n", " inax.tick_params(axis=\"x\", colors=\"white\")\n", " inax.tick_params(axis=\"y\", colors=\"white\")\n", " inax.tick_params(axis=\"z\", colors=\"white\")\n", " inax.set_xlabel('x', labelpad=-10, loc='right')\n", " inax.set_ylabel('y', labelpad=-10, loc='top')\n", " inax.set_zlabel('z', labelpad=-10)\n", " inax.w_xaxis._axinfo.update({'grid': {'color': 'lightgrey', 'linewidth': 0.8, 'linestyle': '-'}})\n", " inax.w_yaxis._axinfo.update({'grid': {'color': 'lightgrey', 'linewidth': 0.8, 'linestyle': '-'}})\n", " inax.w_zaxis._axinfo.update({'grid': {'color': 'lightgrey', 'linewidth': 0.8, 'linestyle': '-'}})\n", "\n", "\n", " fig.show()\n", " f = os.path.join(RESULTS_PLOTS_DIR, f\"{test}_conc_compare.svg\")\n", " ax.figure.savefig(\n", " f,\n", " transparent=True,\n", " dpi=600,\n", " bbox_inches=\"tight\"\n", " )\n", "\n", " # ------------------------------------------------------------------------------------------------------------------\n", " # figure9\n", " # ------------------------------------------------------------------------------------------------------------------\n", " def Fig9(self):\n", " \"\"\"\n", " compare velocity simgraph for a pressure scan\n", " :param cname: colormap name\n", " :param fpath_c:\n", " :param fpath_s:\n", " :param pts:\n", " :param edg:\n", " :return:\n", " \"\"\"\n", " def plot_pressure():\n", " plt1 = Plotter(interactive=False)\n", " vmin, vmax = get_bounds_static(scan=scan, param='pressure', attr='node_attr')\n", "\n", " # manually build a lookup table for mapping colors\n", " lut = buildLUT([\n", " (vmin, 'dg',),\n", " (vmax, 'red'),\n", " ],\n", " vmin=vmin, vmax=vmax,\n", " interpolate=True,\n", " belowColor='yellow',\n", " )\n", " pts1 = pts.clone().pointSize(5).pointColors(node_attr['pressure'], cmap='bwr', vmin=vmin, vmax=vmax)\n", " pts1.addScalarBar3D(\n", " title='(a) pressure (Pa)',\n", " titleFont='VictorMono',\n", " labelFont='VictorMono',\n", " c='k',\n", " )\n", " pts1.scalarbar.addPos(1, 0, 0).useBounds()\n", " edg1 = edg.clone().pointColors(edge_attr['pressure'], cmap='bwr')\n", " for a in pts1.scalarbar.unpack():\n", " if a and a.name == 'Text': a.scale(1.5)\n", "\n", " axes = Axes(edg1, xtitle='y', ytitle='z', ztitle='x', xyGrid=True, zxGrid=True, yzGrid=True) # zLabelRotation=-45)\n", " # plt1.add(pts1, edg1, axes)\n", " f = os.path.join(RESULTS_PLOTS_DIR, f'{test}_pressure.png')\n", " pltp = plt1.show(pts1, edg1, axes).screenshot(f)\n", " return pltp\n", "\n", " def plot_velocity():\n", " \"\"\" line plot of velocity: 3D visualization\"\"\"\n", " vmin, vmax = get_bounds_static(scan=scan, param='velocity', attr='edge_attr', last=False)\n", "\n", " plt2 = Plotter(interactive=False)\n", " # manually build a lookup table for mapping colors\n", " lut = buildLUT([\n", " (vmin, 'b'),\n", " ((vmin+vmax)/2, 'w'),\n", " (vmax+2, 'r'),\n", " ],\n", " vmin=vmin, vmax=vmax+2,\n", " interpolate=True,\n", " aboveColor='dr',\n", " )\n", " color = []\n", " for v in edge_attr['velocity']:\n", " rgb = [0, 0, 0]\n", " lut.GetColor(v, rgb)\n", " color.append(getColorName(rgb))\n", " edge_w = OrderedDict(zip(sorted(G.edges), color))\n", " nx.set_edge_attributes(G, edge_w, 'color')\n", "\n", " nodes = pts.points()\n", " arrsv = []\n", " for i, j, attr in G.edges(data=True):\n", " dv = nodes[j - 1] - nodes[i - 1]\n", " cn = (nodes[j - 1] + nodes[i - 1]) / 2\n", " v = dv / mag(dv) * 3.5\n", " ar = Arrow(cn - v, cn + v, s=.1, c=G[i][j][\"color\"]) # attr['color']\n", " arrsv.append(ar)\n", " edg2 = edg.clone().cmap(lut, edge_attr['velocity'], on='cells').c('grey')\n", " edg2.addScalarBar3D(\n", " title='(b) velocity (\\mum/s)',\n", " titleFont='VictorMono',\n", " labelFont='VictorMono',\n", " c='k',\n", " )\n", " edg2.scalarbar.addPos(0.1, 0, 0).useBounds()\n", " for a in edg2.scalarbar.unpack():\n", " if a and a.name == 'Text': a.scale(1.4)\n", "\n", " axes = Axes(edg2, xtitle='y', ytitle='z', ztitle='x', xyGrid=True, zxGrid=True, yzGrid=True) # zLabelRotation=-45)\n", " # plt2.add(pts, edg2, arrsv, axes)\n", " f = os.path.join(RESULTS_PLOTS_DIR, f'{test}_velocity.png')\n", " pts2 = pts.clone().pointSize(5).c('gray')\n", " pltv = plt2.show(pts2, edg2, arrsv, axes).screenshot(f)\n", "\n", " return pltv\n", "\n", "\n", " input = get_graph(sheet=None)\n", " # create vtk polydata\n", " pos = dict(zip(input['nodes'], input['xyz']))\n", " lines = [(pos[t], pos[h]) for t, h in input['ed_ls']]\n", "\n", " # define points and lines for vtk\n", " pts = Points(input['xyz'], r=10).lighting('off').rotateX(-90)\n", " edg = Lines(lines).lw(5).rotateX(-90)\n", "\n", " scan = get_static_data(fs=RESULTS_SIMGRAPH_FILE, fc=None)\n", " data = scan[\"simgraph\"]\n", " node_attr, edge_attr = data['node_attr'], data['edge_attr']\n", "\n", " pltp = plot_pressure()\n", " pltv = plot_velocity()\n", " f = os.path.join(RESULTS_PLOTS_DIR, f'{test}_p_v.png')\n", " plt12 = Plotter(N=2, interactive=True, sharecam=True)\n", " plt12.show(pltp.actors, at=0, zoom=1.2)\n", " plt12.show(pltv.actors, at=1, zoom=1.2).screenshot(f)\n", "\n", " def Fig10(self):\n", " \"\"\" creates methods figure displaying the cylindrical and spherical volume elements\"\"\"\n", " if False:\n", " G, input, H = load_G_H()\n", " G_nodes = list(G.nodes())\n", " source = int(input['source'])\n", " target = int(input['target'])\n", "\n", " # node shapes and attributes\n", " s_nodes = [n for n in G_nodes if n not in [source, target]] #sphere\n", " c_nodes = list(set(H.nodes) - set(s_nodes))\n", " pos = nx.get_node_attributes(H, 'pos')\n", " node_r = nx.get_node_attributes(H, 'node_r')\n", " node_l = nx.get_node_attributes(H, 'node_l')\n", " # create vtk object\n", " nx_lines = [(pos[t], pos[h]) for t, h in G.edges()]\n", " edg = Lines(nx_lines).lw(5)\n", "\n", " pts = []\n", " for n in H.nodes():\n", " if n in s_nodes:\n", " pts.append(Sphere(r=node_r[n]*0.5, alpha=0.8).pos(pos[n]))\n", " if n in c_nodes:\n", " cyl1 = Cylinder(r=node_r[n]*0.3, height=node_l[n]*.5, cap=False, alpha=1, c='grey').pos(pos[n])\n", " cyl2 = Cylinder(r=node_r[n]*0.5, height=node_l[n]*0.5, cap=False, alpha=0.2, c='grey').pos(pos[n])\n", " cyl = merge(cyl1, cyl2)\n", " pts.append(cyl)\n", " show(pts, edg, axes=True, interactive=True, bg='w', title='plot').screenshot(\"vol\")\n", "\n", " G = nx.OrderedGraph()\n", " G.add_edges_from(ebunch_to_add=[(0, 12), (12, 1), (1, 2), (2, 3), (2, 4), (3, 10), (10, 5), (4, 11), (11, 6)])\n", " # nx_pos = {1: [0, 0, 0], 2: [12.4, 0, 0], 3: [23.14, 6.2, 0], 4: [23.14, -6.2, 0]} #30 degree elevation\n", " # nx_pos = {1: [0, 0, 0], 2: [12.4, 0, 0], 3: [21.17, 8.77, 0], 4: [21.17, -8.77, 0]} #45 degree elevation\n", " nx_pos = {0: [-12.4, 0, 0], 1: [0, 0, 0], 2: [12.4, 0, 0], 3: [18.60, 10.74, 0], 4: [18.60, -10.74, 0], 5: [24.8, 21.47, 0], 6: [24.8, -21.47, 0],\n", " 7: [9.3, 0, 0], 8: [13.95, 2.69, 0], 9: [13.95, -2.69, 0], 10: [21.7, 16.10, 0], 11: [21.7, -16.1, 0],\n", " 12: [-6.2, 0, 0]} #60 degree elevation\n", " # pos = {1: [[-6.2, 0, 0], [6.2, 0, 0]], 2: [12.4, 0, 0], 3: [[17.77, 3.1, 0], [28.5, 9.3, 0]], 4: [[17.77, -3.1, 0], [28.5, -9.3, 0]]} #30 degree elevation\n", " # pos = {1: [[-6.2, 0, 0], [6.2, 0, 0]], 2: [12.4, 0, 0], 3: [[16.78, 4.38, 0], [25.9, 13.15, 0]], 4: [[16.78, -4.38, 0], [25.9, -13.15, 0]]} #45 degree elevation\n", " pos1 = {1: [[-6.2, 0, 0], [6.2, 0, 0]], 2: [12.4, 0, 0], 3: [[15.5, 5.37, 0], [21.7, 16.10, 0]], 4: [[15.5, -5.37, 0], [21.7, -16.1, 0]],\n", " 0: [[-18.6, 0, 0], [-6.2, 0, 0]], 5: [[21.7, 16.10, 0], [27.9, 26.85, 0]], 6: [[21.7, -16.10, 0], [27.9, -26.85, 0]]} #60 degree elevation\n", " r1 = {0: 2.6, 1: 2.6, 2: 4.55, 3: 2.45, 4: 2.6, 5: 2.45, 6: 2.6, 7: 2.6, 8: 2.45, 9: 2.6} #FIXME: revert 2: 4.55\n", " r2 = {1: 8.8, 3: 8.65, 4: 8.8}\n", " color = {1: 'lightsteelblue', 2: 'r', 3: 'darkseagreen', 4: 'grey', 7: 'lightsteelblue', 8: 'darkseagreen',\n", " 9: 'grey'}\n", " nx.set_node_attributes(G, pos1, 'pos')\n", " nx.set_node_attributes(G, r1, 'r1')\n", " nx.set_node_attributes(G, r2, 'r2')\n", "\n", " plt = Plotter(shape=(1, 4), interactive=True, sharecam=0)\n", " # --------------------------------------------------------------------------------------------------------------\n", " # subfigure 0\n", " # --------------------------------------------------------------------------------------------------------------\n", " pos2 = {1: [[-12.4, 0, 0], [12.4, 0, 0]], 3: [[12.4, 0, 0], [24.8, 21.47, 0]], 4: [[12.4, 0, 0], [24.8, -21.47, 0]]}\n", " nx_lines = [(nx_pos[t], nx_pos[h]) for t, h in G.edges()]\n", " edg2 = Lines(nx_lines).lw(2)\n", "\n", " pts2 = []\n", " for n in G.nodes():\n", " if n in [0, 2, 5, 6]:\n", " pts2.append(Sphere(r=0.5, pos=nx_pos[n], c='k'))\n", " elif n in [1, 3, 4]:\n", " pts2.append(Cylinder(r=r1[n], height=12.4, cap=True, pos=pos2[n], c=color[n]).alpha(0.2))\n", "\n", " plt.show(pts2, edg2, axes=False, interactive=True, bg='w', at=0)\n", "\n", " # --------------------------------------------------------------------------------------------------------------\n", " # subfigure 1\n", " # --------------------------------------------------------------------------------------------------------------\n", " pos3 = {7: [[6.2, 0, 0], [12.4, 0, 0]], 8: [[12.4, 0, 0], [15.5, 5.37, 0]], 9: [[12.4, 0, 0], [15.5, -5.37, 0]]}\n", " pts3, edg3 = [], []\n", " for i, j in G.edges():\n", " p1 = nx_pos[i]\n", " p2 = nx_pos[j]\n", " if (i, j) in [(0, 12), (10, 5), (11, 6)]:\n", " edg3.append(DashedLine([p1, p2], spacing=0.3).lw(2))\n", " else:\n", " edg3.append(Line([p1, p2]).lw(2))\n", "\n", " for n in range(0, 10):\n", " if n in [7, 8, 9]:\n", " pts3.append(Sphere(r=0).pos(nx_pos[n]).alpha(0.4))\n", " pts3.append(Cylinder(r=r1[n], height=6.2, cap=True, pos=pos3[n], c=color[n]).alpha(0.4))\n", " elif n in [0, 1, 3, 4, 5, 6]:\n", " pts3.append(Sphere(r=0.5, c='k').pos(nx_pos[n]).alpha(1))\n", " if n in [1, 3, 4]:\n", " pts3.append(Cylinder(r=r1[n], height=12.4, cap=True, pos=pos1[n], c=color[n]).alpha(0.2))\n", " elif n == 2:\n", " pts3.append(Sphere(r=0.5, c='k').pos(nx_pos[n]).alpha(1))\n", "\n", " plt.show(pts3, edg3, axes=False, interactive=True, bg='w', at=1)\n", "\n", " # --------------------------------------------------------------------------------------------------------------\n", " # subfigure 2\n", " # --------------------------------------------------------------------------------------------------------------\n", " pts4, edg4 = [], []\n", " sph1 = []\n", " for i, j in G.edges():\n", " p1 = nx_pos[i]\n", " p2 = nx_pos[j]\n", " if (i, j) in [(0, 12), (10, 5), (11, 6)]:\n", " edg4.append(DashedLine([p1, p2], spacing=0.3).lw(2))\n", " else:\n", " edg4.append(Line([p1, p2]).lw(2))\n", "\n", " for n in G.nodes():\n", " if n == 2:\n", " pts4.append(Sphere(r=r1[n]).pos(pos1[n]).alpha(0.4))\n", " elif n in [1, 3, 4]:\n", " pts4.append(Cylinder(r=r1[n], height=12.4, cap=True, pos=pos1[n], c=color[n]).alpha(0.2))\n", "\n", " if n in [0, 1, 2, 3, 4, 5, 6]:\n", " sph1.append(Sphere(r=0.5, pos=nx_pos[n], c='k'))\n", " else:\n", " sph1.append(Sphere(r=0, pos=nx_pos[n], c='k'))\n", "\n", " # sph1 = [Sphere(r=0.5, pos=nx_pos[n], c='k') for n in G.nodes]\n", " # axes = Axes(edg, xyGrid=True, zxGrid=True, yzGrid=True)\n", " plt.show(pts4, edg4, sph1, axes=False, interactive=True, bg='w', at=2)\n", "\n", " # --------------------------------------------------------------------------------------------------------------\n", " # subfigure 3\n", " # --------------------------------------------------------------------------------------------------------------\n", " pts1, edg1 = [], []\n", " for i, j in G.edges():\n", " p1 = nx_pos[i]\n", " p2 = nx_pos[j]\n", " if (i, j) in [(0, 12), (10, 5), (11, 6)]:\n", " edg1.append(DashedLine([p1, p2], spacing=0.3).lw(2))\n", " else:\n", " edg1.append(Line([p1, p2]).lw(2))\n", "\n", " for n in G.nodes():\n", " if n == 2:\n", " pts1.append(Sphere(r=r1[n]).pos(pos1[n]).alpha(0.4))\n", " elif n in [1, 3, 4]:\n", " cylb = Cylinder(r=r1[n], height=12.4, cap=True, pos=pos1[n], c=color[n]).alpha(0.2)\n", " cylt = Cylinder(r=r2[n], height=12.4, cap=True, pos=pos1[n], c=color[n]).alpha(0.2)\n", " cyl = merge(cylb, cylt)\n", " pts1.append(cyl)\n", "\n", " if n in [0, 1, 2, 3, 4, 5, 6]:\n", " sph1.append(Sphere(r=0.5, pos=nx_pos[n], c='k'))\n", " else:\n", " sph1.append(Sphere(r=0, pos=nx_pos[n], c='k'))\n", "\n", " # sph1 = [Sphere(r=0.5, pos=nx_pos[n], c='k') for n in G.nodes]\n", " # axes = Axes(edg, xyGrid=True, zxGrid=True, yzGrid=True)\n", " plt.show(pts1, edg1, sph1, axes=False, interactive=True, bg='w', at=3).screenshot(\"method\")\n", "\n", "\n", " def Fig11(self):\n", " \"\"\" create interpolated volume for glc, lac concentration data\n", " Example:\n", " -------\n", " G = nx.gnm_random_graph(n=30, m=55, seed=1)\n", " nxpos = nx.spring_layout(G, dim=3)\n", " nxpts = [nxpos[pt] for pt in sorted(nxpos)]\n", " # node values\n", " vals = range(10, 10 + len(nxpts))\n", " print(vals)\n", " nx_lines, vals_segments = [], []\n", " for i, j in G.edges():\n", " nx_lines.append([nxpos[i], nxpos[j]])\n", " vals_segments += [vals[i], vals[j]]\n", " nx_pts = Points(nxpts, r=12)\n", " nx_edg = Lines(nx_lines).lw(2)\n", " nx_pts.cmap('YlGn', vals).addScalarBar()\n", " nx_edg.cmap('YlGn', vals_segments)\n", " labs = nx_pts.labels(vals, scale=0.03, precision=2).c('w')\n", " vol = interpolateToVolume(nx_pts, kernel='shepard', N=2, dims=(50, 50, 50))\n", " slices = []\n", " for i in range(0, 51, 10):\n", " print(i)\n", " sl = vol.xSlice(i).lighting('off').alpha(0.2).cmap('YlGn')\n", " slices.append(sl)\n", " slices.append(sl.isolines(vmin=10, vmax=40))\n", " show(nx_pts, nx_edg, labs, slices, bg='w', axes=9)\n", " # show(labs, slices, nx_pts.scalarbar, bg='w', axes=9, )\n", " \"\"\"\n", " input = get_graph(sheet=\"test11\") #sheet=\"test11_default\"\n", " pos = dict(zip(input['nodes'], input['xyz']))\n", " lines = [(pos[t], pos[h]) for t, h in input['ed_ls']]\n", " pts = Points(input['xyz'], r=10).lighting('off').rotateX(-90)\n", " edg = Lines(lines).lw(5).rotateX(-90)\n", "\n", " mets = {\"glc_ext\": 0, \"glc\": 1, \"lac\": 3, \"lac_ext\": 2}\n", " color = {\"glc_ext\": 'Reds', \"glc\": 'Reds', \"lac_ext\": 'Blues', \"lac\": 'Blues'}\n", "\n", " data = {}\n", " time, scan = get_dynamic_data(fs=RESULTS_SIMGRAPH_FILE)\n", " bounds = graph_utils_py.get_bounds_pscan(scan=scan, common=True)\n", " print(bounds)\n", " x0, x1, y0, y1, z0, z1 = pts.bounds()\n", " for species, species_s in scan['simgraph'].items():\n", " data[species] = dict(zip(time, species_s))\n", "\n", " t = 30.0 #30.0\n", " species = \"glc\"\n", " txt = Text(\n", " f'[{species} ] (mM) t = {t}s',\n", " font='VictorMono',\n", " justify='bottom-center',\n", " s=5, # 40\n", " c='k'\n", " )\n", " txt.pos((x0 + x1) / 2, y0, z1)\n", "\n", " # --------------------------------------------------------------------------------------------------------------\n", " # cell species are plotted at points\n", " # --------------------------------------------------------------------------------------------------------------\n", " def cell_species(species, show=False):\n", " pts1 = pts.clone().cmap(color[species],\n", " data[species][t],\n", " vmin=bounds[species]['simgraph']['lb'],\n", " vmax=bounds[species]['simgraph']['ub'])\n", " pts1.addScalarBar3D(\n", " titleFont='VictorMono',\n", " labelFont='VictorMono',\n", " c='k'\n", " )\n", " pts1.scalarbar.addPos(0.5, 0, 0).useBounds()\n", " for a in pts1.scalarbar.unpack():\n", " if a and a.name == 'Text': a.scale(2.5)\n", "\n", " edg1 = edg.clone()\n", "\n", " # --------------------------------------------------------------------------------------------------------------\n", " # method1\n", " if False:\n", " vol = interpolateToVolume(pts1, kernel='shepard', N=2, dims=(50, 50, 50))\n", " vol.c(color[species]).mode(0).alphaUnit(1).printInfo()\n", " if show: show(vol, pts1, edg1, bg='w', axes=9)\n", " return vol\n", "\n", " # method2\n", " if False:\n", " vol = interpolateToVolume(pts1, radius=5, kernel='linear', N=2)\n", " if species == \"glc\":\n", " iso = vol.isosurface([2, 3, 4, 5]).alpha(0.1).cmap(color[species]) # [1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5.0]\n", " elif species == \"lac\":\n", " iso = vol.isosurface([1, 2, 3, 4]).alpha(0.2).cmap(color[species])\n", " if show: show(iso, pts1, edg1, bg='w', axes=True)\n", " return iso\n", "\n", " # method3\n", " if True:\n", " vol = interpolateToVolume(pts1, radius=5, kernel='linear', N=2, dims=(50, 50, 50))\n", " slices = []\n", " for i in range(0, 51, 5):\n", " sl = vol.xSlice(i).lighting('off').alpha(0.5).cmap(color[species])\n", " slices.append(sl)\n", " # slices.append(sl.isolines(vmin=0, vmax=5))\n", " if show: show(pts1, edg1, slices, bg='w', axes=9)\n", " return slices\n", "\n", " # --------------------------------------------------------------------------------------------------------------\n", " # ext species are plotted at points and interpolated to segments\n", " # --------------------------------------------------------------------------------------------------------------\n", "\n", " def ext_species(species):\n", " pts1 = pts.clone().cmap(color[species],\n", " data[species][t],\n", " vmin=bounds[species]['simgraph']['lb'],\n", " vmax=bounds[species]['simgraph']['ub'])\n", " pts1.pointSize(0)\n", " pts1.addScalarBar3D(\n", " titleFont='VictorMono',\n", " labelFont='VictorMono',\n", " c='k'\n", " )\n", " pts1.scalarbar.addPos(0.1, 0, 0).useBounds()\n", " for a in pts1.scalarbar.unpack():\n", " if a and a.name == 'Text': a.scale(2.5)\n", "\n", " edg1 = edg.clone().cmap(color[species],\n", " get_vals_segment(point_data=data[species][t], sheet=\"test11\"),\n", " vmin=bounds[species]['simgraph']['lb'],\n", " vmax=bounds[species]['simgraph']['ub'])\n", " return pts1, edg1\n", " # show(pts1, edg1, axes=True, bg2='white', at=0, interactive=True)\n", "\n", " vol_glc = cell_species(species=\"glc\")\n", " pts_eglc, edg_eglc = ext_species(species=\"glc_ext\")\n", " vol_lac = cell_species(species=\"lac\")\n", " pts_elac, edg_elac = ext_species(species=\"lac_ext\")\n", "\n", " plt12 = Plotter(N=2, interactive=True)\n", " plt12.show(pts_eglc, edg_eglc, vol_glc, at=0) #vol_glc\n", " plt12.show(pts_elac, edg_elac, vol_lac, at=1).screenshot(\"vol2\") #vol_lac\n", "\n", " def Fig12(self):\n", " \"\"\"\n", " :param self:\n", " :return:\n", " reference:\n", " ---------\n", " https://github.com/marcomusy/vedo/blob/master/examples/volumetric/probeLine2.py\n", " \"\"\"\n", " input = get_graph(sheet=\"test11\")\n", " pos = dict(zip(input['nodes'], input['xyz']))\n", " lines = [(pos[t], pos[h]) for t, h in input['ed_ls']]\n", " pts = Points(input['xyz'], r=10).lighting('off')\n", " edg = Lines(lines).lw(5)\n", "\n", " color = {\"glc_ext\": 'Reds', \"glc\": 'Reds', \"lac_ext\": 'Blues', \"lac\": 'Blues'}\n", " time, scan = get_dynamic_data(fs=RESULTS_SIMGRAPH_FILE)\n", " bounds = graph_utils_py.get_bounds_pscan(scan=scan, common=True)\n", "\n", " data = {}\n", " for species, species_s in scan['simgraph'].items():\n", " data[species] = dict(zip(time, species_s))\n", " ts = [0.0, 1.2, 4.2, 15, 30, 60, 120, 300]\n", "\n", " y = {}\n", " for t in ts:\n", " species = \"glc\"\n", " pts1 = pts.clone().cmap(color[species],\n", " data[species][t],\n", " vmin=bounds[species]['simgraph']['lb'],\n", " vmax=bounds[species]['simgraph']['ub'])\n", " pts1.addScalarBar3D(\n", " titleFont='VictorMono',\n", " labelFont='VictorMono',\n", " c='k'\n", " )\n", " pts1.scalarbar.addPos(0.1, 0, 0).useBounds()\n", " for a in pts1.scalarbar.unpack():\n", " if a and a.name == 'Text': a.scale(2.5)\n", "\n", " vol = interpolateToVolume(pts1, radius=5, kernel='linear', N=2, dims=(50, 50, 50))\n", " p1, p2 = pos[35], pos[38]\n", " pl = probeLine(vol, p1, p2, res=10).lineWidth(4).color('r')\n", "\n", " x = pl.points()\n", " y.update({t: pl.getPointArray()})\n", "\n", " ax = plt.subplot(111)\n", " for t in ts:\n", " ax.plot(\n", " x[:, 0], y[t],\n", " # color=color[n],\n", " markersize=5,\n", " )\n", " ax.set_xlabel('Time (s)', fontsize=10)\n", " ax.set_ylabel('Concentration (mM)', fontsize=10)\n", " plt.show()\n", "\n", " show(pts, edg, pl, bg='w')\n", "\n", " def Fig13(self):\n", " \"\"\"plots steady-state time as node values and interpolates\"\"\"\n", "\n", " def plot_risetime(key):\n", " data = scan[key]\n", " node_attr, edge_attr = data['node_attr'], data['edge_attr']\n", "\n", " vplt = Plotter(interactive=False)\n", " # vmin, vmax = get_bounds_static(scan=scan, param='tss_glc_ext', attr='node_attr')\n", "\n", " pts1 = pts.clone().pointSize(10).pointColors(node_attr['tss_glc_ext'], cmap='Spectral')\n", " pts1.addScalarBar3D(\n", " title='(a) Rise time (s)',\n", " titleFont='VictorMono',\n", " labelFont='VictorMono',\n", " c='k',\n", " )\n", " pts1.scalarbar.addPos(1, 0, 0).useBounds()\n", "\n", " label = [str(n) for n in input['nodes']]\n", " labs = pts.labels(label, scale=5, font='VictorMono', c='k').addPos(0.05, 0, 0.5)\n", "\n", " edg1 = edg.clone().pointColors(edge_attr['tss_glc_ext'], cmap='Spectral')\n", " for a in pts1.scalarbar.unpack():\n", " if a and a.name == 'Text': a.scale(1.5)\n", "\n", " if test == \"test11_default\":\n", " axes = Axes(edg1, xtitle='y', ytitle='z', ztitle='x', xyGrid=True, zxGrid=True, yzGrid=True)\n", " else:\n", " axes = True\n", "\n", " vplt_return = vplt.show(pts1, edg1, labs, axes)\n", " return vplt_return\n", "\n", " input = get_graph(sheet=\"test24\")\n", " # create vtk polydata\n", " pos = dict(zip(input['nodes'], input['xyz']))\n", " lines = [(pos[t], pos[h]) for t, h in input['ed_ls']]\n", " pts = Points(input['xyz'], r=10).lighting('off').rotateX(0) #-90 for test11\n", " edg = Lines(lines).lw(5).rotateX(0)\n", "\n", " scan = get_static_data(fs=RESULTS_SIMGRAPH_FILE, fc=RESULTS_COMSOL_FILE, twitch=True)\n", " vplt_sim = plot_risetime(key=\"simgraph\")\n", " vplt_com = plot_risetime(key=\"comsol\")\n", "\n", " f = os.path.join(RESULTS_PLOTS_DIR, f'{test}_risetime.png')\n", " plt12 = Plotter(N=2, interactive=True, sharecam=True)\n", " plt12.show(vplt_sim.actors, at=0, zoom=1.2)\n", " plt12.show(vplt_com.actors, at=1, zoom=1.2).screenshot(f)\n", "\n", " def Fig14(self):\n", " \"\"\"plots rise time as node values and interpolates\"\"\"\n", "\n", " def risetime(sheet):\n", " input = get_graph(sheet=sheet)\n", " pos = dict(zip(input['nodes'], input['xyz']))\n", " lines = [(pos[t], pos[h]) for t, h in input['ed_ls']]\n", " pts = Points(input['xyz'], r=10).lighting('off').rotateX(0)\n", " edg = Lines(lines).lw(5).rotateX(0)\n", "\n", " data = scan[sheet]\n", " node_attr, edge_attr = data['node_attr'], data['edge_attr']\n", "\n", " pts1 = pts.clone().pointSize(0).cmap(lut, node_attr['rt_glc_ext'])\n", " if sheet == \"test9\":\n", " pts1.addScalarBar3D(\n", " title='Rise time (s)',\n", " titleFont='VictorMono',\n", " labelFont='VictorMono',\n", " c='k',\n", " )\n", " pts1.scalarbar.addPos(40, 0, 0).useBounds()\n", " for a in pts1.scalarbar.unpack():\n", " if a and a.name == 'Text': a.scale(1.5)\n", "\n", " label = [str(n) for n in input['nodes']]\n", " labs = pts.labels(label, scale=5, font='VictorMono', c='k').addPos(0.05, 0, 0.5)\n", "\n", " edg1 = edg.clone().cmap(lut, edge_attr['rt_glc_ext'])\n", "\n", " f = os.path.join(RESULTS_PLOTS_DIR, f'{sheet}_risetime.png')\n", " vplt = Plotter(interactive=True)\n", " vplt_return = vplt.show(pts1, edg1, axes=True, zoom=zoom[sheet]).screenshot(f)\n", " return vplt_return\n", "\n", " scan = {}\n", " tag = {\"test10\": \"design2\", \"test9\": \"design1\"}\n", " zoom = {\"test10\": 1, \"test9\": 1}\n", " for s in tag.keys():\n", " fs = os.path.join(RESULTS_SIMGRAPH_DIR, f\"{s}_default_bc1_v1_c0.mat\")\n", " r = get_static_data(fs=fs, twitch=True, sheet=s)\n", " scan[s] = r[\"simgraph\"]\n", "\n", " vmin, vmax = get_bounds_static(scan=scan, param='rt_glc_ext', attr='node_attr') # vmin=0, vmax=2034.1\n", "\n", " # manually build a lookup table for mapping colors\n", " lut = buildLUT([\n", " (vmin, 'db'),\n", " (15, 'b'),\n", " (35, 'lb'),\n", " (54, 'w'),\n", " (vmax, 'r'),\n", " ],\n", " vmin=vmin, vmax=vmax,\n", " interpolate=True\n", " )\n", "\n", " plt1 = risetime(sheet=\"test10\")\n", " plt2 = risetime(sheet=\"test9\")\n", " fpath = os.path.join(RESULTS_PLOTS_DIR, 'risetime.png')\n", " plt12 = Plotter(N=2, interactive=True, sharecam=False)\n", " plt12.show(plt1.actors, at=0, zoom=1, axes=True)\n", " plt12.show(plt2.actors, at=1, zoom=1, axes=True).screenshot(fpath)\n", "\n", " def Fig15(self):\n", " \"\"\" creates time-series plots of blood and cell species for islet vasculature\"\"\"\n", " settings.showRendererFrame = False\n", "\n", " input = get_graph(sheet=\"test11\")\n", " nodes = input['nodes']\n", " pos = dict(zip(nodes, input['xyz']))\n", " lines = [(pos[t], pos[h]) for t, h in input['ed_ls']]\n", " pts = Points(input['xyz'], r=10).lighting('off').rotateX(-90)\n", " edg = Lines(lines).lw(5).rotateX(-90)\n", " rxn_nodes = [1, 11, 15, 17, 18, 24, 28, 30, 33, 37, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51,\n", " 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73,\n", " 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95,\n", " 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112]\n", " sph_r = OrderedDict.fromkeys(nodes, 0)\n", " sph_r.update(OrderedDict.fromkeys(rxn_nodes, 3))\n", "\n", " color = {\"glc_ext\": 'Reds', \"glc\": 'Reds', \"lac_ext\": 'Blues', \"lac\": 'Blues'}\n", " label = {\"glc_ext\": 'eglucose', \"glc\": 'glucose',\n", " \"lac_ext\": 'elactate', \"lac\": 'lactate'}\n", " data = {}\n", " time, scan = get_dynamic_data(fs=RESULTS_SIMGRAPH_FILE)\n", " ts = [0.6, 6.0, 30.0]\n", " plt = Plotter(shape=(len(ts), nspecies), interactive=0, sharecam=1, bg2='white') #\n", "\n", " # ------------------------------------------------------------------------------------------------------------------\n", " x0, x1, y0, y1, z0, z1 = pts.bounds()\n", " bounds = graph_utils_py.get_bounds_pscan(scan=scan, common=False)\n", "\n", " at = 0\n", " for idx, key in enumerate(scan):\n", " for i in range(nspecies):\n", " s = SPECIES[i]\n", " c = scan[key][s]\n", " data[s] = dict(zip(time, c))\n", "\n", " for t in ts:\n", " for species in SPECIES:\n", " txt = Text(\n", " f't: {t}s',\n", " font='VictorMono',\n", " justify='bottom-left',\n", " s=5,\n", " c='k'\n", " )\n", " txt.pos((x0 + x1) / 2, y0, z1)\n", " vmin = bounds[species]['simgraph']['lb']\n", " vmax = bounds[species]['simgraph']['ub']\n", "\n", " if species in [\"glc\", \"lac\"]:\n", " temp = []\n", " for value in data[species][t]:\n", " value_c = colorMap(\n", " value,\n", " name=color[species],\n", " vmin=vmin,\n", " vmax=vmax\n", " )\n", " temp.append(value_c)\n", " sph_c = OrderedDict(zip(nodes, temp))\n", "\n", " # pts\n", " pts1 = pts.clone().cmap(\n", " color[species],\n", " data[species][t],\n", " vmin=vmin,\n", " vmax=vmax\n", " )\n", "\n", " sph1 = [Sphere(r=sph_r[n]).color(sph_c[n]).pos(pos[n]).alpha(1).rotateX(-90) for n in nodes]\n", " if ts.index(t) <= len(ts)-1:\n", " pts1.addScalarBar3D(\n", " title=f'[{label[species]}] (mM)',\n", " titleFont='VictorMono',\n", " labelFont='VictorMono',\n", " c='k'\n", " )\n", " pts1.scalarbar.addPos(15, 0, 0).useBounds()\n", " for a in pts1.scalarbar.unpack():\n", " if a and a.name == 'Text': a.scale(2.5)\n", "\n", " # edg\n", " edg1 = edg.clone()\n", "\n", " else:\n", " pts1 = pts.clone().pointSize(2).cmap(\n", " color[species],\n", " data[species][t],\n", " vmin=vmin,\n", " vmax=vmax\n", " )\n", " edg1 = edg.clone().cmap(\n", " color[species],\n", " get_vals_segment(point_data=data[species][t], sheet=\"test11\"),\n", " vmin=vmin,\n", " vmax=vmax\n", " )\n", "\n", " if ts.index(t) <= len(ts)-1:\n", " edg1.addScalarBar3D(\n", " title=f'[{label[species]}] (mM)',\n", " titleFont='VictorMono',\n", " labelFont='VictorMono',\n", " c='k'\n", " )\n", " edg1.scalarbar.addPos(15, 0, 0).useBounds()\n", " for a in edg1.scalarbar.unpack():\n", " if a and a.name == 'Text':\n", " a.scale(2)\n", "\n", " axes = Axes(edg1, xtitle='y', ytitle='z', ztitle='x', xyGrid=True, zxGrid=False, yzGrid=False)\n", " if species in [\"glc\", \"lac\"]:\n", " plt.show(\n", " # pts1,\n", " sph1,\n", " edg1,\n", " pts1.scalarbar,\n", " txt,\n", " axes,\n", " at=at,\n", " zoom=1,\n", " )\n", " else:\n", " plt.show(\n", " pts1,\n", " edg1,\n", " txt,\n", " axes,\n", " at=at,\n", " zoom=1,\n", " )\n", "\n", " at = at + 1\n", " interactive()\n", " plt.screenshot(os.path.join(RESULTS_PLOTS_DIR, f'{test}_conc_bvcell.png'))\n", "\n", " def Fig16(self):\n", " \"\"\"\n", " plot glucose and lactate net exchange flux\n", " Reference:\n", " https://stackoverflow.com/questions/14762181/adding-a-y-axis-label-to-secondary-y-axis-in-matplotlib\n", " https://stackoverflow.com/questions/58490203/matplotlib-plot-a-line-with-open-markers-where-the-line-is-not-visible-within\n", " ---------\n", " Example:\n", " x = np.arange(0, 10, 0.1)\n", " y1 = 0.05 * x ** 2\n", " y2 = -1 * y1\n", " \"\"\"\n", "\n", " data = io_utils_py._get_file(file=RESULTS_GSCAN_FILE, spreadsheet=False)\n", "\n", " dose, glcim, lacex = [], [], []\n", " for s in data['variable']:\n", " dose.append(s.dose)\n", " glcim.append(s.glcim_net)\n", " lacex.append(abs(s.lacex_net))\n", "\n", " fig, ax = plt.subplots()\n", " # fig, ax1 = plt.subplots()\n", " # ax2 = ax1.twinx()\n", " ax.plot(dose, glcim, 'r', marker='o', markerfacecolor='none', markersize=4)\n", " ax.plot(dose, lacex, 'b', marker='s', markerfacecolor='none', markersize=4)\n", "\n", " ax.set_xlabel('Glucose dose (mM)', fontsize=12)\n", " ax.set_ylabel('Net exchange flux (pmole/min)', fontsize=12)\n", " # ax2.set_ylabel('Lactate release (pmole/min)', fontsize=12)\n", "\n", " fig.show()\n", " f = os.path.join(RESULTS_PLOTS_DIR, f\"{test}_flux.svg\")\n", " ax.figure.savefig(\n", " f,\n", " transparent=True,\n", " dpi=600,\n", " bbox_inches=\"tight\"\n", " )\n", "\n", " def Fig17(self):\n", " \"\"\"\n", " compare conc simgraph vs. comsol for glc_ext evolution in islet vasculature\n", " :param pts:\n", " :param edg:\n", " :return:\n", " ref: https://stackoverflow.com/questions/44947505/how-to-make-a-movie-out-of-images-in-python\n", " Note to self: If the colormap appears black/blank for cell species, check input file\n", " \"\"\"\n", "\n", " def plot_timeseries(scan):\n", "\n", " images = []\n", " data = {}\n", " for key, species_s in scan.items():\n", " data_s = dict(zip(time, species_s[\"glc_ext\"]))\n", " data[key] = data_s\n", "\n", " # ts = [0.0, 1.2, 4.2]\n", " bounds = graph_utils_py.get_bounds_pscan(scan=scan, common=True)\n", " x0, x1, y0, y1, z0, z1 = pts.bounds()\n", " plt = Plotter(N=2, interactive=False)\n", "\n", " for t in time:\n", " for key in scan.keys():\n", " txt = Text3D(\n", " f'[glc_ext ] (mM) t = {round(t,2)}s',\n", " font='VictorMono',\n", " justify='bottom-center',\n", " s=5,\n", " c='k'\n", " )\n", " txt.pos((x0 + x1) / 2, y0, z1)\n", "\n", " # pts\n", " pts1 = pts.clone()\n", " pts1.pointSize(1)\n", "\n", " # edg\n", " edg1 = edg.clone().cmap(\n", " \"Reds\",\n", " get_vals_segment(point_data=data[key][t], sheet=\"test11\"),\n", " vmin=bounds[\"glc_ext\"][key]['lb'],\n", " vmax=bounds[\"glc_ext\"][key]['ub']\n", " )\n", " edg1.addScalarBar3D(\n", " title=key,\n", " titleFont='VictorMono',\n", " labelFont='VictorMono',\n", " c='k',\n", " titleXOffset=3.5\n", " )\n", " edg1.scalarbar.addPos(0.1, 0, 0).useBounds()\n", " for a in edg1.scalarbar.unpack():\n", " if a and a.name == 'Text': a.scale(2.5)\n", "\n", " if key == \"simgraph\":\n", " plt.show(\n", " pts1,\n", " edg1,\n", " txt,\n", " axes=True,\n", " bg2='white',\n", " at=0,\n", " # zoom=1.5,\n", " )\n", " elif key == \"comsol\":\n", " plt.show(\n", " pts1,\n", " edg1,\n", " txt,\n", " axes=True,\n", " bg2='white',\n", " at=1,\n", " # zoom=1.5,\n", " )\n", "\n", " fpath = os.path.join(RESULTS_PLOTS_DIR, f'{test}_conc_glc_ext_{round(t,2)}s.png')\n", " images.append(fpath)\n", " plt.screenshot(fpath)\n", "\n", " return images\n", "\n", " input = get_graph(sheet=\"test11\")\n", " pos = dict(zip(input['nodes'], input['xyz']))\n", " lines = [(pos[t], pos[h]) for t, h in input['ed_ls']]\n", " pts = Points(input['xyz'], r=10).lighting('off').rotateX(-90)\n", " edg = Lines(lines).lw(5).rotateX(-90)\n", "\n", " time, scan = get_dynamic_data(fs=RESULTS_SIMGRAPH_FILE, fc=RESULTS_COMSOL_FILE)\n", "\n", " # create video\n", " images = plot_timeseries(scan=scan)\n", " print(images)\n", " video_name = os.path.join(RESULTS_PLOTS_DIR, \"islet_compare.mp4\")\n", " frame = cv2.imread(os.path.join(RESULTS_PLOTS_DIR, images[0]))\n", " height, width, layers = frame.shape\n", " fourcc = 0\n", " fps = 1\n", " video = cv2.VideoWriter(video_name, fourcc, fps, (width, height))\n", " for image in images:\n", " fpath = os.path.join(RESULTS_PLOTS_DIR, image)\n", " video.write(cv2.imread(fpath))\n", " os.remove(fpath) # remove the image file\n", " cv2.destroyAllWindows()\n", " video.release()\n", "\n", " def Fig18(self):\n", " \"\"\"\n", " compare conc simgraph vs. comsol for glc_ext evolution in mesentery vasculature\n", " :param pts:\n", " :param edg:\n", " :return:\n", " ref: https://stackoverflow.com/questions/44947505/how-to-make-a-movie-out-of-images-in-python\n", " Note to self: If the colormap appears black/blank for cell species, check input file\n", " \"\"\"\n", "\n", " def plot_timeseries(scan, sheet):\n", "\n", " images = []\n", " data = {}\n", " for key, species_s in scan.items():\n", " data_s = dict(zip(time, species_s[\"glc_ext\"]))\n", " data[key] = data_s\n", "\n", " # ts = [0.0, 1.2, 4.2]\n", " bounds = graph_utils_py.get_bounds_pscan(scan=scan, common=True)\n", " x0, x1, y0, y1, z0, z1 = pts.bounds()\n", " plt = Plotter(N=2, interactive=False)\n", "\n", "\n", " for t in time:\n", " for key in scan.keys():\n", " txt = Text3D(\n", " f'[glc_ext ] (mM) t = {round(t,2)}s',\n", " font='VictorMono',\n", " justify='bottom-center',\n", " s=150,\n", " c='k'\n", " )\n", " txt.pos((x0 + x1) / 2, y0, 5000)\n", "\n", " # pts\n", " pts1 = pts.clone()\n", " pts1.pointSize(1)\n", "\n", " # edg\n", " edg1 = edg.clone().cmap(\n", " lut,\n", " get_vals_segment(point_data=data[key][t], sheet=sheet),\n", " vmin=bounds[\"glc_ext\"][key]['lb'],\n", " vmax=bounds[\"glc_ext\"][key]['ub']\n", " )\n", " edg1.addScalarBar3D(\n", " title=key,\n", " titleFont='VictorMono',\n", " labelFont='VictorMono',\n", " c='k',\n", " titleXOffset=2.5\n", " )\n", " edg1.scalarbar.addPos(0.1, 0, 0).useBounds()\n", " for a in edg1.scalarbar.unpack():\n", " if a and a.name == 'Text': a.scale(1)\n", "\n", " if key == \"simgraph\":\n", " plt.show(\n", " pts1,\n", " edg1,\n", " txt,\n", " axes=True,\n", " bg2='white',\n", " at=0,\n", " zoom=1.3,\n", " )\n", " elif key == \"comsol\":\n", " plt.show(\n", " pts1,\n", " edg1,\n", " txt,\n", " axes=True,\n", " bg2='white',\n", " at=1,\n", " zoom=1.3,\n", " )\n", "\n", " fpath = os.path.join(RESULTS_PLOTS_DIR, f'{test}_conc_glc_ext_{round(t,2)}s.png')\n", " images.append(fpath)\n", " plt.screenshot(fpath)\n", "\n", " return images\n", "\n", " input = get_graph(sheet=\"test24\")\n", " pos = dict(zip(input['nodes'], input['xyz']))\n", " lines = [(pos[t], pos[h]) for t, h in input['ed_ls']]\n", " pts = Points(input['xyz'], r=10).lighting('off')\n", " edg = Lines(lines).lw(5)\n", "\n", " # ------------------------------------------------------------------------------------------------------------------\n", " # manually build a lookup table for mapping colors\n", " # ------------------------------------------------------------------------------------------------------------------\n", " lut = buildLUT([\n", " (0, 'w'),\n", " (5, 'red'),\n", " ],\n", " vmin=0,\n", " vmax=5,\n", " interpolate=True,\n", " aboveColor='darkred',\n", " belowColor='w'\n", " )\n", " # ------------------------------------------------------------------------------------------------------------------\n", "\n", " time, scan = get_dynamic_data(fs=RESULTS_SIMGRAPH_FILE, fc=RESULTS_COMSOL_FILE)\n", "\n", " # create video\n", " images = plot_timeseries(scan=scan, sheet=\"test24\")\n", " video_name = os.path.join(RESULTS_PLOTS_DIR, \"mesentery_compare.mp4\")\n", " frame = cv2.imread(os.path.join(RESULTS_PLOTS_DIR, images[0]))\n", " height, width, layers = frame.shape\n", " fourcc = 0\n", " fps = 1\n", " video = cv2.VideoWriter(video_name, fourcc, fps, (width, height))\n", " for image in images:\n", " fpath = os.path.join(RESULTS_PLOTS_DIR, image)\n", " video.write(cv2.imread(fpath))\n", " os.remove(fpath) # remove the image file\n", " cv2.destroyAllWindows()\n", " video.release()\n", "\n", "\n", " def Fig19(self):\n", " \"\"\"\n", " compare conc simgraph vs. comsol for glc_ext evolution in tumor vasculature\n", " :param pts:\n", " :param edg:\n", " :return:\n", " ref: https://stackoverflow.com/questions/44947505/how-to-make-a-movie-out-of-images-in-python\n", " Note to self: If the colormap appears black/blank for cell species, check input file\n", " \"\"\"\n", "\n", " def plot_timeseries(scan, sheet, t):\n", "\n", " input = get_graph(sheet=sheet)\n", " pos = dict(zip(input['nodes'], input['xyz']))\n", " lines = [(pos[t], pos[h]) for t, h in input['ed_ls']]\n", " pts = Points(input['xyz'], r=10).lighting('off').rotateX(0)\n", " edg = Lines(lines).lw(5).rotateX(0)\n", "\n", " data = dict(zip(time, scan[\"simgraph\"][\"glc_ext\"]))\n", "\n", " bounds = graph_utils_py.get_bounds_pscan(scan=scan, common=True)\n", " x0, x1, y0, y1, z0, z1 = pts.bounds()\n", " plt = Plotter(N=2, interactive=False)\n", " label = {\"test9\": \"Tumor design1\", \"test10\": \"Tumor design2\"}\n", "\n", " for key in scan.keys():\n", " txt = Text3D(\n", " f'[glc_ext ] (mM) t = {round(t,2)}s',\n", " font='VictorMono',\n", " justify='bottom-center',\n", " s=25,\n", " c='k'\n", " )\n", " txt.pos((x0 + x1) / 2, y0, z1+500)\n", "\n", " # pts\n", " pts1 = pts.clone()\n", " pts1 = pts.clone().cmap(\n", " \"Reds\",\n", " data[t],\n", " vmin=bounds[\"glc_ext\"][key]['lb'],\n", " vmax=bounds[\"glc_ext\"][key]['ub']\n", " )\n", " pts1.pointSize(1)\n", " # edg\n", " edg1 = edg.clone().cmap(\n", " \"Reds\",\n", " get_vals_segment(point_data=data[t], sheet=sheet),\n", " vmin=bounds[\"glc_ext\"][key]['lb'],\n", " vmax=bounds[\"glc_ext\"][key]['ub']\n", " )\n", " edg1.addScalarBar3D(\n", " title=label[sheet],\n", " titleFont='VictorMono',\n", " labelFont='VictorMono',\n", " c='k',\n", " titleXOffset=2.5\n", " )\n", "\n", " edg1.scalarbar.addPos(0.1, 0, 0).useBounds()\n", " for a in edg1.scalarbar.unpack():\n", " if a and a.name == 'Text': a.scale(1.5)\n", "\n", " vplt_return = plt.add(\n", " [pts1,\n", " edg1,\n", " txt],\n", " # axes=True,\n", " # bg2='white',\n", " at=0,\n", " # zoom=1.5,\n", " )\n", "\n", " return vplt_return\n", "\n", " scan = {}\n", " time, scan[\"design1\"] = get_dynamic_data(fs=os.path.join(RESULTS_SIMGRAPH_DIR, \"test9_default_bc1_v1_c0.mat\"))\n", " time, scan[\"design2\"] = get_dynamic_data(fs=os.path.join(RESULTS_SIMGRAPH_DIR, \"test10_default_bc1_v1_c0.mat\"))\n", " images = []\n", " for t in time:\n", " plt1 = plot_timeseries(scan=scan[\"design1\"], sheet=\"test9\", t=t)\n", " plt2 = plot_timeseries(scan=scan[\"design2\"], sheet=\"test10\", t=t)\n", " fpath = os.path.join(RESULTS_PLOTS_DIR, f'tumor_conc_glc_ext_{round(t,2)}s.png')\n", " plt12 = Plotter(N=2, interactive=False, sharecam=True)\n", " plt12.show(plt1.actors, at=0, zoom=1.2, axes=True, bg2='gray')\n", " plt12.show(plt2.actors, at=1, zoom=1.2, axes=True, bg2='gray').screenshot(fpath)\n", " images.append(fpath)\n", "\n", " # create video\n", " video_name = os.path.join(RESULTS_PLOTS_DIR, \"tumor_compare.mp4\")\n", " frame = cv2.imread(os.path.join(RESULTS_PLOTS_DIR, images[0]))\n", " height, width, layers = frame.shape\n", " fourcc = 0\n", " fps = 1\n", " video = cv2.VideoWriter(video_name, fourcc, fps, (width, height))\n", " for image in images:\n", " fpath = os.path.join(RESULTS_PLOTS_DIR, image)\n", " video.write(cv2.imread(fpath))\n", " os.remove(fpath) # remove the image file\n", " cv2.destroyAllWindows()\n", " video.release()\n", "\n", "\n", "if __name__ == '__main__':\n", "\n", " # ------------------------------------------------------------------------------------------------------------------\n", " # Generate geometry\n", " # ------------------------------------------------------------------------------------------------------------------\n", " # G = create_graph(directed=True, attributes=True)\n", " # G = get_eucledian_lengths_wrt_origin(origin=34, edge_ed=False)\n", " # G, H = mesh_network(G=G)\n", " # draw_graph3d_vedo(Gs=[G], label_scale=5, highlight=False, default_label=True)\n", " # exit()\n", " manuscript_figures = ManuscriptFigures()\n", " # manuscript_figures.figures()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.11.1" } }, "nbformat": 4, "nbformat_minor": 5 }