{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# **GOLEM MFR tomography**\n", "This notebook loads the camera data from two certical and radial fast cameras on tokamak GOLEM and then performs tomopgraphic inversion using MFR algorithm.\n", "Geometry matrix data are needed - notebook called geometry_matrix.ipynb creates all needed.\n", " \n", "\n", "*Required python libraries for this notebook: numpy, matplotlib,, opencv, tomotok, IPython.display and scipy.ndimage*\n", "________________________________________________" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Code:\n", "Necessary imports:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "tags": [] }, "outputs": [], "source": [ "from typing import Optional\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import cv2 #opencv\n", "import requests\n", "import os, glob\n", "from pathlib import Path\n", "from scipy.ndimage import gaussian_filter\n", "from matplotlib.ticker import FuncFormatter\n", "\n", "from tomotok.core.derivative import standard_anisotropic_derivative_matrices, all_direction_derivative_matrices\n", "#from tomotok.core.derivative import compute_iso_dmats, compute_aniso_dmats\n", " \n", "from tomotok.core.geometry import RegularGrid\n", "from tomotok.core.geometry.io import load_sparse_gmat\n", "from tomotok.core.inversions import CholmodMfr, Mfr, Fixt, CholmodFixt\n", "\n", "from IPython.display import Image\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "\n", "shot_no = \n", "channel = 0 #color channel to use (0,1,2) = (B,G,R)\n", "\n", "camsConf = requests.get(f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/FastCameras/OutputCamsConfig.json').json()\n", "FPS = camsConf[0]['recordRate'] #camera frame rate\n", "FPMS = FPS * 1e-3 #camera frame rate [1/ms]\n", "\n", "current_folder = Path(f'/home/golem/FastCameras/{shot_no}')\n", "\n", "result_folder = Path(current_folder / 'results')\n", "Path(result_folder).mkdir(exist_ok=True)\n", "Path(result_folder / \"data\").mkdir(exist_ok=True)\n", "Path(result_folder / \"graphs\").mkdir(exist_ok=True)\n", "Path(result_folder / \"stats\").mkdir(exist_ok=True)\n", "\n", "frames_folder = Path(current_folder / 'frames')\n", "Path(frames_folder).mkdir(exist_ok=True)\n", "Path(frames_folder / \"vertical\").mkdir(exist_ok=True)\n", "Path(frames_folder / \"radial\").mkdir(exist_ok=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### **Defined functions:**\n", "Function creates mask for poloidal plane dataset:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "tags": [] }, "outputs": [], "source": [ "def create_mask(grid: RegularGrid, limiter_mask: bool, plasmaradius: Optional[float] = None):\n", " \"\"\"Creates mask, of the same size as grid, of True and False vlaues to define where to compute the reconstructions.\n", "\n", " Args:\n", " grid (RegularGrid): Reconstruction grid.\n", " limiter_mask (bool): True if you want the results up to limiter, False if you want to see results to LCFS.\n", " plasmaradius (float, optional): Radius of plasma (and of the reconstruction domain). If limiter_mask is False, this parameter can be None.\n", "\n", " Returns:\n", " RegularGrid class ndarray of the given grid size, with True values where the reconstrucition will be computed (inside limiter/plasma radius).\n", " \"\"\"\n", " \n", " plasmacentre = 0.40 #minor tokamak radius\n", " chamberradius = 0.1 #chamber radius\n", " \n", " #Limiter and chamber radii definitions\n", " t = np.linspace(0,2*np.pi,100) #t is just a parameter for the curves/circles\n", "\n", " if limiter_mask==True:\n", " #limiter curve definition\n", " r = plasmacentre + plasmaradius*np.sin(t)\n", " z = plasmaradius*np.cos(t)\n", " bdm = grid.is_inside(r, z)\n", " else:\n", " #chamber curve definition\n", " R = plasmacentre + chamberradius*np.sin(t)\n", " Z = chamberradius*np.cos(t)\n", " bdm = grid.is_inside(R, Z)\n", " return bdm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Function that extracts images from loaded video data:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def imagesfromvideo(minrange:int,maxrange:int):\n", " \"\"\"Extracts single frames from both GOLEM camera data and saves them\n", "\n", " Args:\n", " ver_vid (cv2.VideoCapture): vertical camera data\n", " rad_vid (cv2.VideoCapture): radial camera data\n", " Returns:\n", " \"\"\"\n", " #loads from golem shot website - .mp4\n", " capver=cv2.VideoCapture(f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/FastCameras/Camera_Vertical/Data.mp4')\n", " caprad=cv2.VideoCapture(f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/FastCameras/Camera_Radial/Data.mp4')\n", " i = minrange \n", " \n", " capver.set(1,i)\n", " caprad.set(1,i)\n", " \n", " retver, framever = capver.read(1) \n", " retrad, framerad = caprad.read(1) \n", " while retver and retrad and i < maxrange:\n", " retver, framever = capver.read(1) \n", " retrad, framerad = caprad.read(1) \n", " \n", " #Saving frames\n", " namever = str(frames_folder / \"vertical\" / str(str(i)+'.png'))\n", " cv2.imwrite(namever, framever) \n", " namerad = str(frames_folder / \"radial\" / str(str(i)+'.png'))\n", " cv2.imwrite(namerad, framerad) \n", " i=i+1\n", " capver.release()\n", " caprad.release()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "tags": [] }, "outputs": [], "source": [ "def datafromimage(framenumber:int,channel: int):\n", " \"\"\"Extracts signal (frame given by frame number) from both GOLEM camera frames. \n", "\n", " Args:\n", " framenumber (int): Number of frame which will be extracted\n", " channel (int): color channel to be used 012 BGR\n", "\n", "\n", " Returns:\n", " sig: Signal from given frame in ndarray form.\n", " \"\"\"\n", " namever = str(frames_folder / \"vertical\" / str(str(framenumber)+'.png'))\n", " namerad = str(frames_folder / \"radial\" / str(str(framenumber)+'.png'))\n", "\n", " ogimagerad = cv2.imread(namerad,cv2.IMREAD_UNCHANGED )[:,:,channel] #loading radial image \n", " ogimagever = cv2.imread(namever,cv2.IMREAD_UNCHANGED )[:,:,channel] #loading vertical image\n", "\n", " #THIS MUST BE CHANGED DEPENDING ON THE GEOMETRY MATRIX USED - MUST CONSTRUCT IT FOR EACH ROW DIFFERENTLTY\n", " rowrad = 511 #which radial camera detector row to use for reconstruction\n", " rowver = 511 #which vertical camera detector row to use for reconstruction\n", "\n", " #Discarding obstructed detectors\n", " leftstoprad = 0 #leftmost R camera detector to be used (useful if cameras' view is obstructed on edges, 0 if all detectors see the plasma)\n", " rightstoprad = 1280 #rightmost R camera detector to be used (1280 if all)\n", " totalspanrad = rightstoprad - leftstoprad\n", "\n", " leftstopver = 0 #leftmost V camera detector to be used (useful if cameras' view is obstructed on edges, 0 if all detectors see the plasma)\n", " rightstopver = 1280 #rightmost V camera detector to be used (1280 if all)\n", " totalspanver = rightstopver - leftstopver\n", "\n", " fullheight = 1024 #height of full camera picture in pixels\n", " fullwidth = 1280 #width of full camera picture in pixels\n", " imgheight = np.shape(ogimagerad)[0] #height of loaded image, \n", " imgheighthalf = int((fullheight-imgheight)/2) #where to insert loaded image into 1024 by 1280 array for easier handling of distortion\n", "\n", " \n", " imagerad = np.zeros((fullheight-imgheight, fullwidth))\n", " imagever = np.zeros((fullheight-imgheight, fullwidth))\n", " imagerad = np.insert(imagerad, imgheighthalf, ogimagerad, axis=0)\n", " imagever = np.insert(imagever, imgheighthalf, ogimagever, axis=0)\n", " \n", " #imagever = np.flip(imagever)\n", " #imagerad = np.flip(imagerad)\n", "\n", " imagerad = imagerad[rowrad:rowrad+1,leftstoprad:rightstoprad] #taking signal from only one row of camera detectors\n", " imagever = imagever[rowver:rowver+1,leftstopver:rightstopver]\n", " imagever = np.squeeze(imagever,axis=0)\n", " imagerad = np.squeeze(imagerad,axis=0)\n", " \n", " #imagever[390:450] = 1.05*gaussian_filter(imagever, 30)[390:450]\n", " \n", " imagerad = imagerad\n", " #imagerad[970:] = gaussian_filter(imagerad, 30)[970:]\n", " \n", " sig = np.concatenate((imagever, imagerad), axis=0)\n", " \n", " \n", " return sig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Function where the tomographic inversion is performed:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "tags": [] }, "outputs": [], "source": [ "def datasolve(framenumber: int, sig, grid, gmat, magflux, mask, fixed:bool, alpha=[0.7,0.7,0.7], initial_guess = None):\n", " \"\"\" Tomographic reconstruction using MFR of single frame data, which saves result\n", "\n", " Args:\n", " framenumber (int): framenumber\n", " sig (_type_): camera signal to be reconstructed, most likely from datafromimage()\n", " grid (_type_): reconstruction grid\n", " magflux (_type_): magnetic flux surfaces\n", " mask (_type_): mask for the reconstruction grid, reconstruction only within limiter right now \n", " fixed (bool): if the function should use the fixed parameter set True (could speed up the reconstruction)\n", "\n", " Returns:\n", " _type_: np.ndarray of reconstructed plane\n", " \"\"\" \n", " \n", " data = sig.reshape(1, -1) # data should have shape (#timeslices, #channels/pixels)\n", " \n", " errors = .02*data + 5 #error estimation. This works, but should be improved and studied further. Has very large impact on results\n", "\n", "\n", " dmats = standard_anisotropic_derivative_matrices(grid, flux=magflux, mask=mask) #derivative matrix with mask\n", " aniso = [1, 0.1, 1, 0.1] #parameter dictating anisotropy of smoothing perp/parallel to mag flux surfaces\n", "\n", " #Alternative for older Tomotok versions:\n", " #dmats = compute_aniso_dmats(grid, magflux=magflux, mask=mask) #derivative matrix with mask\n", " #aniso = 2.2 #parameter dictating anisotropy of smoothing perp/parallel to mag flux surfaces\n", " #initial_guess = None\n", "\n", " gmat_sel = gmat[:, mask.flatten()] #geometric matrix with mask\n", "\n", "\n", " if fixed == True:\n", " \n", " #Solver that uses a fixed regularization parameter and is therefore faster (about 10 times). \n", " try:\n", " solver = Fixt() #CholmodFixt()\n", " print(\"using CholmodFixt\")\n", " except ImportError:\n", " solver = Fixt()\n", " print(\"using Fixt\")\n", " \n", " #res, chi = solver(data, gmat_sel, dmats, errors, aniso=aniso, parameters = alpha)\n", " res, chi = solver(data, gmat_sel, dmats, errors, derivative_weights=aniso, parameters = alpha, initial_guess = initial_guess)\n", " print('chi = {:.2f}'.format(chi[-1]['chi'][-1]))\n", "\n", " else: \n", " #Solvers that use full Minimum Fischer Regularization. Slower but \"better\".\n", " try: \n", " solver = Mfr() # CholmodMfr() requires scikit.sparse\n", " print(\"using CholmodMFR\")\n", " except ImportError:\n", " solver = Mfr() # slowlier version based on scipy.sparse\n", " print(\"using MFR\")\n", " \n", " #res, chi = solver(data, gmat_sel, dmats, errors, aniso=aniso)\n", " res, chi = solver(data, gmat_sel, dmats, errors, derivative_weights=aniso, initial_guess = initial_guess)\n", " \n", " print('last reg parameter logalpha = {:.2f}'.format(chi[-1]['logalpha'][-1]))\n", " print('chi = {:.2f}'.format(chi[-1]['chi'][-1]))\n", "\n", "\n", " print(np.shape(res))\n", " #reshaping result to original grid shape\n", " res2 = np.zeros((1, grid.shape[0]*grid.shape[1]))\n", " res2[0, mask.flatten()] = res\n", " res2 = res2.reshape(1, *grid.shape)\n", " np.save(result_folder / 'data' / str(\"res_\"+ str(framenumber)), res2, allow_pickle=True, fix_imports=True)\n", " np.save(result_folder / 'stats' / str(\"stats_\"+ str(framenumber)), chi, allow_pickle=True, fix_imports=True)\n", " return res2, chi\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### **Main part of the code:**\n", "Here some parameters are defined, like plasma radius etc. to set boundaries for example" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "tags": [] }, "outputs": [], "source": [ "# Define some \"parameters\" used in tokamak GOLEM\n", "\n", "#Parameters:\n", "plasmaradius=0.085 #Radius of plasma (and of the reconstruction domain), 0.085 for limiter. \n", "padding = 0.015 #padding to work with slightly bigger grid for better visualization\n", "#GOLEM parameters:\n", "chamberradius = 0.1 #chamber radius\n", "plasmacentre = 0.40 #minor tokamak radius\n", "\n", "#Curves for limiter contures etc.:\n", "t = np.linspace(0,2*np.pi,100) #t is just a parameter for the curves/circles\n", "#Plasma (determined by limiter)\n", "r = plasmacentre + plasmaradius*np.sin(t)\n", "z = plasmaradius*np.cos(t)\n", "#Chamber wall\n", "R = plasmacentre + chamberradius*np.sin(t)\n", "Z = chamberradius*np.cos(t)\n", "\n", "#Grid minimum and maximum locations in poloidal plane of the GOLEM tokamak for better plotting:\n", "gridrmin = plasmacentre - plasmaradius - padding\n", "gridrmax = plasmacentre + plasmaradius + padding\n", "gridzmin = - plasmaradius - padding\n", "gridzmax = plasmaradius + padding\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Load geometry matrix and grid:**\n", "\n", "The geometry matrix describes the geometry of the field of view of the cameras in relation to the reconstruction grid. The geometry matrix element [i,j] describes how much the emissivity in the i-th region contributes to the signal at the j-th detector. Here, this is done using the line of sight approximation in which the [i,j] element is proportional to the length of the line of sight of the j-th detector intersecting the i-th region.\n", "\n", "The geometry matrix is specific to the cameras used and to the grid size. If any of these are changed, a new geometry matrix has to be generated.\n", "\n", "This geometry matrix was created using Tomotok from Calcam camera calibrations. It is specifically for the UX100 color cameras and their 511th pixel row and a 50x50 reconstruction grid. \n" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "tags": [] }, "outputs": [], "source": [ "# Loading sparse_gmat (class with gmat and grid) and extracting data\n", "gmat_load=load_sparse_gmat(f'/home/golem/FastCameras/{shot_no}/gmat/gmat_40x40') #loading data\n", "gmat=gmat_load[0] # Defining geometry matrix gmat\n", "grid=gmat_load[1] # Defining grid\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A graph showing the cameras' field of view using the sum of geometry matrix elements. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "if True: #set True if you want to plot the sum of geometry matrix elements\n", " plt.rcParams['font.size'] = '15'\n", "\n", " f, a = plt.subplots(figsize=(5,5),dpi=100)\n", " \n", " \n", " # img = a.imshow(res.reshape(grid.shape), origin='lower', extent=grid.extent) # nodes must be squares\n", " img = a.pcolorfast(grid.r_border, grid.z_border, gmat.toarray().sum(0).reshape(grid.shape))\n", " plt.plot(r,z, 'w', lw=1)\n", " plt.plot(R,Z, 'w', lw=1)\n", " a.set_aspect(1)\n", " a.set_xlabel('R [m]')\n", " a.set_ylabel('z [m]')\n", "\n", " plt.title('Geometry matrix')\n", "\n", " plt.xticks([gridrmin,gridrmax,0.4,0.45,0.35])\n", " plt.yticks([gridzmin,gridzmax,0.,0.05,-0.05])\n", " #plt.xticks(rotation = 90)\n", "\n", " cbar = f.colorbar(img, ax=a, fraction=0.0452, pad=0.05)\n", " cbar.ax.set_ylabel('[-]')\n", "\n", " #imagename='geometric_matrix.png'\n", " #plt.savefig(str(current_folder / \"imagename\"),bbox_inches=\"tight\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Create mask:**\n", "\n", "The mask is created to limit the reconstruction domain to only the inside of the tokamak chamber. " ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "tags": [] }, "outputs": [], "source": [ "# Calling mask function for defined grid\n", "bdm=create_mask(grid, limiter_mask=False, plasmaradius=plasmaradius)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "tags": [] }, "outputs": [], "source": [ "if False: #set True if you want to plot the mask\n", " plt.figure()\n", " plt.imshow(bdm)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Generating Psi:**\n", "\n", "The magnetic surfaces are used for anisotropical smoothing of the results. \n", "\n", "Currently, the Solovyev solution of the Grad-Shafranov equation is used." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "tags": [] }, "outputs": [], "source": [ "#Psi generation\n", "nr = grid.nr\n", "nz = grid.nz\n", "span = plasmaradius+padding\n", "offsetx = 0.01 #offset of centre in x direction \n", "offsety = 0 #offset of centre in y direction \n", "x = np.linspace(-span, span, nr)\n", "y = np.linspace(-span, span, nz)\n", "mx, my = np.meshgrid(x, y)\n", "\n", "#Simplest possible Psi - concentric circles\n", "#ipsi = np.sqrt((mx -offsetx)*( mx-offsetx) + (my-offsety) *(my-offsety))\n", "\n", "#Simplest ~physical Psi - Solovyev solution of G-S equation\n", "R0 = 0.4\n", "a = 0.085\n", "eps = a/R0\n", "tau = 1\n", "sigma = 1\n", "ipsi = ( ((mx)/a) - (eps/2)*(1-((mx)/a)**2))**2 + (1-eps**2/4)*(my/a)**2\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "if True: #set True if you want to plot the \"Psi\" - smoothing of the result\n", " plt.rcParams['font.size'] = '15'\n", " f, a = plt.subplots(figsize=(5,5),dpi=100)\n", " # img = a.imshow(res.reshape(grid.shape), origin='lower', extent=grid.extent) # nodes must be squares\n", " img = a.pcolorfast(grid.r_border, grid.z_border, ipsi)\n", " plt.plot(r,z, 'w', lw=1)\n", " plt.plot(R,Z, 'k', lw=1)\n", "\n", " a.set_aspect(1)\n", " a.set_xlabel('R [m]')\n", " a.set_ylabel('z [m]')\n", "\n", " plt.title('')\n", "\n", " plt.xticks([gridrmin,gridrmax,0.4,0.45,0.35])\n", " plt.yticks([gridzmin,gridzmax,0.,0.05,-0.05])\n", " #plt.xticks(rotation = 90)\n", "\n", " cbar = f.colorbar(img, ax=a, fraction=0.0452, pad=0.05,format = \"{x:.2}\")\n", " cbar.ax.set_ylabel('Psi [-]')\n", "\n", " #imagename='mag_flux.png'\n", " #plt.savefig(str(current_folder / \"imagename\"),bbox_inches=\"tight\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set the range of frames you want to compute:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "#Tomographic reconstructions of all frames selected\n", "plasma_start=float(requests.get(f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/PlasmaDetection/Results/t_plasma_start').text)\n", "plasma_end=float(requests.get(f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/PlasmaDetection/Results/t_plasma_end').text)\n", "\n", "minrange = int(plasma_start*FPMS)\n", "maxrange = int(plasma_end*FPMS)+1\n", "\n", "imagesfromvideo(minrange,maxrange)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "tags": [] }, "outputs": [], "source": [ "if False: #set True if you want to test and see the image data\n", " #Frame numbers to reconstruct\n", " frame = 278 #starting frame number\n", " data = datafromimage(frame,0)\n", " plt.xlabel('Detector number')\n", " plt.ylabel('Signal')\n", " plt.plot(data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute and immediately plot all tomography results:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "tags": [] }, "outputs": [], "source": [ "alpha = (0,0,0)\n", "\n", "for i in range(minrange,maxrange):\n", " if i%50 == 0:\n", " print(\"Tomo reconstruction number \"+str(i))\n", " \n", " data = datafromimage(i,channel=channel)\n", " \n", " if i == minrange:\n", " res, stats = datasolve(i, data, grid, gmat, ipsi, bdm,fixed=False)\n", " logalpha = stats[0]['logalpha']\n", " alpha = (10**logalpha[0], 10**logalpha[1], 10**logalpha[2])\n", " chi = stats[0]['chi'][-1]\n", " else:\n", " initial_guess = np.load(result_folder / 'data' / str(\"res_\"+ str(i-1)+ \".npy\"))\n", " initial_guess = initial_guess[0, bdm]\n", " \n", " res, stats = datasolve(i, data, grid, gmat, ipsi, bdm, fixed=True, alpha = alpha, initial_guess=initial_guess )\n", " \n", " chi = stats[0]['chi'][-1]\n", " if abs(chi - 1) > 0.1:\n", " res, stats = datasolve(i, data, grid, gmat, ipsi, bdm, fixed=False, initial_guess=initial_guess)\n", " logalpha = stats[0]['logalpha']\n", " alpha = (10**logalpha[0], 10**logalpha[1], 10**logalpha[2])\n", " chi = stats[0]['chi'][-1]\n", "\n", " \n", " plt.figure()\n", " f, a = plt.subplots(figsize=(5,5),dpi=50)\n", " img = a.pcolorfast(grid.r_border, grid.z_border, np.load(result_folder / 'data' / str(\"res_\"+ str(i)+ \".npy\"), mmap_mode=None, allow_pickle=False, fix_imports=True)[0], snap=True)\n", " plt.plot(r,z, 'w', lw=2)\n", " plt.plot(R,Z, 'k', lw=2)\n", " a.set_aspect(1)\n", " a.set_xlabel('R [m]')\n", " a.set_ylabel('z [m]')\n", " \n", " time = (i)/FPMS\n", "\n", " plt.title('$t_d$ = '+\"{:.3f}\".format(time)+' ms')\n", "\n", " plt.xticks([gridrmin,gridrmax,0.4,0.45,0.35])\n", " plt.yticks([gridzmin,gridzmax,0.,0.05,-0.05])\n", " #plt.xticks(rotation = 90)\n", " \n", " fmt = lambda x, pos: '{:1.1e}'.format(x)\n", " cbar = f.colorbar(img, ax=a, fraction=0.0452, pad=0.05,format = FuncFormatter(fmt))\n", " cbar.ax.set_ylabel('[-]')\n", "\n", " imagename= result_folder / 'graphs' / str(\"res_\"+ str(i)+ \".png\")\n", " plt.savefig(imagename,bbox_inches=\"tight\",transparent=True)\n", " plt.close('all') \n", " #Magflux based on centre of gravity of previous frame, should be improved\n", " #offsetx = np.sum(res * mx)/np.sum(res)\n", " #offsety = np.sum(res * my)/np.sum(res)\n", " #ipsi = np.sqrt((mx -offsetx)*( mx-offsetx) + (my-offsety) *(my-offsety))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Replot all results with the same color scale:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "#Plot all reconstructions\n", "\n", "#Finding maximum in all results_iso - to have a stationary colorbar\n", "maxres= 0\n", "maxres = np.zeros(maxrange-minrange)\n", "\n", "for i in range(minrange,maxrange):\n", " res2 = np.load(result_folder / 'data' / str(\"res_\"+ str(i)+ \".npy\"), mmap_mode=None, allow_pickle=True, fix_imports=True)\n", " maxres[i-int(minrange)] = np.max(res2)\n", " \n", "maxres = gaussian_filter(maxres, 20)\n", " \n", "plt.rcParams['font.size'] = '15'\n", "\n", "\n", "for i in range(minrange,maxrange):\n", " f, a = plt.subplots(figsize=(5,5),dpi=100)\n", " img = a.pcolorfast(grid.r_border, grid.z_border, np.load(result_folder / 'data' / str(\"res_\"+ str(i)+ \".npy\"), mmap_mode=None, allow_pickle=True, fix_imports=True)[0], vmax=maxres[i-minrange])\n", " plt.plot(r,z, 'w', lw=1)\n", " plt.plot(R,Z, 'k', lw=1)\n", " a.set_aspect(1)\n", " a.set_xlabel('R [m]')\n", " a.set_ylabel('z [m]')\n", " \n", " time = (i)/FPMS\n", "\n", " plt.title('$t_d$ = '+\"{:.3f}\".format(time)+' ms')\n", "\n", " plt.xticks([gridrmin,gridrmax,0.4,0.45,0.35])\n", " plt.yticks([gridzmin,gridzmax,0.,0.05,-0.05])\n", " #plt.xticks(rotation = 90)\n", " \n", " fmt = lambda x, pos: '{:1.1e}'.format(x)\n", " cbar = f.colorbar(img, ax=a, fraction=0.0452, pad=0.05,format = FuncFormatter(fmt))\n", " cbar.ax.set_ylabel('[-]')\n", "\n", " imagename=result_folder / 'graphs' / str(\"res_\"+ str(i)+ \".png\")\n", " plt.savefig(imagename,bbox_inches=\"tight\",transparent=True)\n", " plt.close('all')\n", " \n", " if i%50 == 0:\n", " print(i)\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "if False:\n", " index = 476\n", " \n", " to_plot = np.load(result_folder / 'data' / str(\"res_\"+ str(index)+ \".npy\"), mmap_mode=None, allow_pickle=True, fix_imports=True)[0]\n", " to_plot = to_plot[bdm]\n", "\n", " to_plot2 = np.zeros((1, grid.shape[0]*grid.shape[1]))\n", " to_plot2[0, bdm.flatten()] = to_plot\n", " to_plot2 = to_plot2.reshape(1, *grid.shape)\n", "\n", "\n", " f, a = plt.subplots(figsize=(5,5),dpi=100)\n", " img = a.pcolorfast(grid.r_border, grid.z_border, to_plot2[0])\n", " plt.plot(r,z, 'w', lw=1)\n", " plt.plot(R,Z, 'k', lw=1)\n", " a.set_aspect(1)\n", " a.set_xlabel('R [m]')\n", " a.set_ylabel('z [m]')\n", "\n", " time = (i)/FPMS\n", "\n", " plt.title('$t_d$ = '+\"{:.3f}\".format(time)+' ms')\n", "\n", " plt.xticks([gridrmin,gridrmax,0.4,0.45,0.35])\n", " plt.yticks([gridzmin,gridzmax,0.,0.05,-0.05])\n", " #plt.xticks(rotation = 90)\n", "\n", " cbar = f.colorbar(img, ax=a, fraction=0.0452, pad=0.05)\n", " cbar.ax.set_ylabel('[-]')\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Combine all graphs into an animation:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "img_array = []\n", "\n", "#create iconfig for homepage\n", "i = int(0.5*(maxrange-minrange)+minrange)\n", "img = cv2.imread(str(result_folder / 'graphs' / str(\"res_\"+ str(i)+ \".png\")))\n", "cv2.imwrite(str(current_folder / 'tomo_iconfig.png'), img) \n", "cv2.imwrite(str(current_folder / 'tomo_ScreenShotAll.png'), img) \n", "\n", "for i in range(minrange,maxrange):\n", " img = cv2.imread(str(result_folder / 'graphs' / str(\"res_\"+ str(i)+ \".png\")))\n", " img_array.append(img)\n", "\n", "size = (np.shape(img)[1],np.shape(img)[0])\n", "\n", "out = cv2.VideoWriter(str(current_folder / 'TomographyResult.avi'),cv2.VideoWriter_fourcc(*'DIVX'), 10, size)\n", "\n", "for i in range(len(img_array)):\n", " out.write(img_array[i])\n", "out.release()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Hollowness parameter computation for use in database" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "def quality_parameters(data):\n", " \"\"\"\n", " Calculates a quality parameter based on the symmetry of values in a 2D array with respect to a specified center point.\n", "\n", " The function evaluates symmetry by comparing the maximum values along vertical (above vs. below), \n", " horizontal (left vs. right), and diagonal directions (top-left vs. bottom-right, top-right vs. bottom-left)\n", " passing through a specified center point. The quality parameter increases when symmetry is high and decreases \n", " when asymmetry is present.\n", "\n", " Args:\n", " data (np.ndarray): A 2D numpy array representing the data.\n", "\n", " Returns:\n", " tuple: A tuple containing the calculated quality parameter, symmetry score, and hollowness score.\n", "\n", " Raises:\n", " IndexError: If the data array is not large enough to include the center point.\n", " ZeroDivisionError: If the value at the center point is zero.\n", " \"\"\"\n", " # Define the center point\n", " center = (21, 20)\n", "\n", " # Check the dimensions of the input array\n", " if data.shape[0] <= center[0] or data.shape[1] <= center[1]:\n", " raise IndexError(\"Input data array must be at least of \\size (22, 21).\")\n", " \n", " # Retrieve the value at the center point\n", " center_value = data[center[0], center[1]]\n", " \n", " # Check if the center value is zero to avoid division by zero, handle gracefully\n", " if center_value == 0:\n", " raise ZeroDivisionError(\"The value at the center point is zero.\")\n", "\n", " # Maximum values above and below in the column\n", " max_above = np.max(data[:center[0], center[1]]) if center[0] > 0 else 0\n", " max_below = np.max(data[center[0]+1:, center[1]]) if center[0] < data.shape[0] - 1 else 0\n", " \n", " # Maximum values left and right in the row\n", " max_left = np.max(data[center[0], :center[1]]) if center[1] > 0 else 0\n", " max_right = np.max(data[center[0], center[1]+1:]) if center[1] < data.shape[1] - 1 else 0\n", "\n", " # Diagonal comparison: top-left vs. bottom-right and top-right vs. bottom-left\n", " # Top-left quadrant\n", " top_left = data[:center[0], :center[1]]\n", " max_top_left = np.max(top_left) if top_left.size > 0 else 0\n", " \n", " # Bottom-right quadrant\n", " bottom_right = data[center[0]+1:, center[1]+1:]\n", " max_bottom_right = np.max(bottom_right) if bottom_right.size > 0 else 0\n", " \n", " # Top-right quadrant\n", " top_right = data[:center[0], center[1]+1:]\n", " max_top_right = np.max(top_right) if top_right.size > 0 else 0\n", " \n", " # Bottom-left quadrant\n", " bottom_left = data[center[0]+1:, :center[1]]\n", " max_bottom_left = np.max(bottom_left) if bottom_left.size > 0 else 0\n", "\n", " # Symmetry measures (inverse of asymmetry):\n", " vertical_symmetry = 1 - abs(max_above - max_below) / max(max_above, max_below, 1) # Penalize vertical asymmetry\n", " horizontal_symmetry = 1 - abs(max_left - max_right) / max(max_left, max_right, 1) # Penalize horizontal asymmetry\n", " diagonal_symmetry_1 = 1 - abs(max_top_left - max_bottom_right) / max(max_top_left, max_bottom_right, 1) # Top-left vs. bottom-right\n", " diagonal_symmetry_2 = 1 - abs(max_top_right - max_bottom_left) / max(max_top_right, max_bottom_left, 1) # Top-right vs. bottom-left\n", " \n", " # Total symmetry score (higher if symmetrical):\n", " total_symmetry = (vertical_symmetry + horizontal_symmetry + diagonal_symmetry_1 + diagonal_symmetry_2) / 4\n", "\n", " # Calculate hollowness based on the relation of total maximums to the center value\n", " total_max = (max_above + max_below + max_left + max_right + max_top_left + max_bottom_right + max_top_right + max_bottom_left) / 8\n", " total_hollowness = total_max / center_value\n", "\n", " # Calculate the quality parameter, adjusted by the symmetry score\n", " parameter = total_hollowness * total_symmetry\n", " \n", " return parameter, total_symmetry, total_hollowness" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Initialize lists to hold your parameters\n", "frames= []\n", "parameters = []\n", "symmetries = []\n", "hollownesses = []\n", "\n", "# Loop through the range\n", "for i in range(minrange, maxrange):\n", " # Load the data\n", " data = np.load(result_folder / 'data' / str(\"res_\"+ str(i)+ \".npy\"), mmap_mode=None, allow_pickle=True, fix_imports=True)[0]\n", " \n", " # Extract the parameters\n", " parameter, total_symmetry, total_hollowness = quality_parameters(data)\n", " \n", " # Append the parameters to the lists\n", " frames.append(int(i))\n", " parameters.append(parameter)\n", " symmetries.append(total_symmetry)\n", " hollownesses.append(total_hollowness)\n", "\n", "# Convert lists to a NumPy array for easier saving\n", "results_array = np.column_stack((frames, parameters, symmetries, hollownesses))\n", "\n", "\n", "total_frames=int(maxrange-minrange)\n", "# Calculate mean values\n", "mean_values = [np.mean(parameters), np.mean(symmetries), np.mean(hollownesses)]\n", "# Define column names\n", "header = 'Frame,Parameter,Total_Symmetry,Total_Hollowness'\n", "\n", "# Save the array to a text file\n", "with open(result_folder / 'quality_parameters.txt', 'w') as f:\n", " np.savetxt(f, results_array, header=header, comments='', delimiter=',', fmt='%f')\n", " f.write('\\nTotal frames and Mean Values:\\n')\n", " f.write(f'Total_frames,Parameter_Mean,Total_Symmetry_Mean,Total_Hollowness_Mean\\n')\n", " f.write(f'{total_frames},{mean_values[0]},{mean_values[1]},{mean_values[2]}\\n')" ] } ], "metadata": { "interpreter": { "hash": "823b433449232ffa024257f9a4d323559a5ef8fa1cde7207ec0907d532cc106a" }, "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.7.3" } }, "nbformat": 4, "nbformat_minor": 4 }