import os, glob, requests, cv2
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy import ndimage
from PIL import Image
from io import BytesIO
from IPython.display import Markdown
from IPython.display import Image as DispImage
from scipy.interpolate import interp1d
from matplotlib.ticker import (MultipleLocator, FormatStrFormatter, AutoMinorLocator)
ds=np.DataSource()
class StopExecution(Exception):
def _render_traceback_(self):
pass
def remove_old(name):
if os.path.exists(name):
os.remove(name)
files = ['analysis.html', 'icon-fig.png','plasma_film_Radial.png','plasma_film_Vertical.png']
for file in files:
remove_old(file)
if os.path.exists('Results/'):
file = glob.glob('Results/*')
for f in file:
os.remove(f)
shot_no =0
Vertical = True; Radial = True
if not os.path.exists('Camera_Radial/Data.mp4'):
Radial = False
if not os.path.exists('Camera_Vertical/Data.mp4'):
Vertical = False
if not Vertical and not Radial: raise StopExecution
camsConf = requests.get(f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/FastCameras/OutputCamsConfig.json').json()
FPS = {}
for cams in camsConf:
FPS[cams['name']] = cams['recordRate']
To accurately remap the pixels to the position at the center of the vacuum chamber, we will utilize the data obtained from the calibration of the cameras. For comprehensive details on calibration procedures, please refer to [BP, J. Chlum] or [PRPL, J.Chlum & M.Odlozilik].
# Load geometric info on camera's line of sight
CamsRad_LOS = np.load('CamsRad_chords_coord.npz')
CamsVert_LOS = np.load('CamsVert_chords_coord.npz')
shiftx = -0.4 # shift in x axis; to position the chamber at the center with x = 0
def interp_LOS(i_chord, cams, index_coord, data_coord, value):
f_LOS = interp1d(cams[index_coord][i_chord,:], cams[data_coord][i_chord,:])
return float(f_LOS(value))
camsRad_centerx = np.array([interp_LOS(i_chord, CamsRad_LOS, 'z', 'x', 0) for i_chord in range(CamsRad_LOS['nchords'])])
camsRad_centerx += shiftx
camsVert_centerz = np.array([interp_LOS(i_chord, CamsVert_LOS, 'x', 'z', 0.4) for i_chord in range(CamsVert_LOS['nchords'])])
def plot_toroid(elements, R, r):
u = np.linspace(0, 2*np.pi, elements)
v = np.linspace(0, 2*np.pi, elements)
uu, vv = np.meshgrid(u, v)
xx = (R+r*np.cos(vv))*np.cos(uu)
yy = (R+r*np.cos(vv))*np.sin(uu)
zz = r*np.sin(vv)
coord = [xx, yy, zz]
return coord
# Create additional toroidal structures
elements = 100
coord_plasma = plot_toroid(elements, 0.4, 0.085)
coord_chamber = plot_toroid(elements, 0.4, 0.095)
# Parameters to plot only a small part of the torus
a = 10; b = 2
fig = plt.figure(figsize = (10,10))
ax = fig.add_subplot(221, projection='3d')
ax2 = fig.add_subplot(222)
# Plot boundary defined by limiter
ax.plot_surface(coord_plasma[0][:,:a]+shiftx, coord_plasma[1][:,:a], coord_plasma[2][:,:a], antialiased=True, color='orange', edgecolors='k', lw = 0.05, alpha = 0.4)
ax.plot_surface(coord_plasma[0][:,-b:]+shiftx, coord_plasma[1][:,-b:], coord_plasma[2][:,-b:], antialiased=True, color='orange', edgecolors='k', lw = 0.05, alpha = 0.4)
limiter = plt.Circle((0, 0), 0.085, color='tab:orange', fill = False); ax2.add_patch(limiter)
# Plot chamber
ax.plot_surface(coord_chamber[0][:,:a]+shiftx, coord_chamber[1][:,:a], coord_chamber[2][:,:a], antialiased=True, color='grey', edgecolors='k', lw = 0.05, alpha = 0.4)
ax.plot_surface(coord_chamber[0][:,-b:]+shiftx, coord_chamber[1][:,-b:], coord_chamber[2][:,-b:], antialiased=True, color='grey', edgecolors='k', lw = 0.05, alpha = 0.4)
ax.set(xlabel = ('x [m]'), ylabel = ('y [m]'), zlabel = ('z [m]'))
chamber = plt.Circle((0, 0), 0.095, color = 'grey', fill = False); ax2.add_patch(chamber);
sc_unit = 1 # default is [m]
n_chords = 50 # Plot every 50th chord
# Note that chord with index 0 corresponds to first chord at LFS and top, respectively
# Plot chords
for cams, color in zip([CamsRad_LOS,CamsVert_LOS], ['tab:red', 'tab:blue']):
for x,y,z in zip(cams['x'][::n_chords,:], cams['y'][::n_chords,:], cams['z'][::n_chords,:]):
ax.plot3D(x*sc_unit+shiftx, y*sc_unit, z*sc_unit, color = color)
ax2.plot(x*sc_unit+shiftx, z*sc_unit, color = color, alpha = 0.5, lw =1)
# Plot cuts at the center of the chamber
ax.plot3D(camsRad_centerx, np.zeros_like(camsRad_centerx),np.zeros_like(camsRad_centerx), c = 'k', lw = 1)
ax.plot3D(np.zeros_like(camsRad_centerx),np.zeros_like(camsRad_centerx),camsVert_centerz, c = 'k', lw = 1)
ax2.plot(camsRad_centerx, np.zeros_like(camsRad_centerx), c = 'k', lw = 2)
ax2.plot(np.zeros_like(camsRad_centerx),camsVert_centerz, c = 'k', lw = 2)
# Set axes
ax.axes.set_xlim3d(left=-0.4, right=0.4)
ax.axes.set_ylim3d(bottom=-0.4, top=0.4)
ax.axes.set_zlim3d(bottom=-0.4, top=0.4)
edge = 0.2
ax2.set(xlim = (-edge,edge), ylim = (-edge,edge),xlabel = 'r [m]', ylabel = 'z [m]', aspect = 'equal');
plt.tight_layout(pad = 6)
def get_frame(file, Position='Radial'):
global nFrames
cv2.destroyAllWindows()
video = cv2.VideoCapture(file)
current=0
try:
os.mkdir('Camera_'+Position+'/Frames')
except FileExistsError:
filelist = glob.glob(os.path.join('Camera_'+Position+'/Frames/', "*"))
for f in filelist:
os.remove(f)
while (video.isOpened()):
ret, frame=video.read()
if ret:
name = 'Camera_' + Position + '/Frames/%i.png'%current
idx_fr = np.shape(frame)[0]//2
if Position == 'Radial':
cv2.imwrite(name, frame[idx_fr-1:idx_fr,:,:]) #frame[height,width]
else:
cv2.imwrite(name, frame[idx_fr-1:idx_fr,:,:]) #frame[height,width]
current += 1
else:
break
nFrames=current
video.release()
cv2.destroyAllWindows()
ret, frame=video.read()
return
get_frame('Camera_Radial/Data.mp4','Radial')
get_frame('Camera_Vertical/Data.mp4','Vertical')
def img_shape(img_name,Position='Radial'):
img=cv2.imread('Camera_'+Position+'/Frames/%s.png'%img_name)
height=img.shape[0]
width=img.shape[1]
return height, width
img_shape('0')
(1, 1280)
def make_image(FramStart,FramEnd, Position = 'Radial'):
nFrames=FramEnd-FramStart
frames=[]
width, height = img_shape('0',Position) #width is new height
img_all = np.zeros((height, (FramEnd-FramStart)*width,3), np.uint8)
i=0
for frame in range(FramStart,FramEnd+1):
if Position == 'Radial':
img=cv2.rotate(cv2.imread('Camera_'+Position+f'/Frames/{frame}.png'),cv2.ROTATE_90_CLOCKWISE)
elif Position == 'Vertical':
img=cv2.rotate(cv2.imread('Camera_'+Position+f'/Frames/{frame}.png'),cv2.ROTATE_90_CLOCKWISE)
img_all[:height,i*width:(i+1)*width,:3]=img
i+=1
print('Height:', height, '\nWidth:', width,'\nimg_all shape',img_all.shape)
cv2.imwrite('plasma-film_' + Position + '.png',img_all)
cv2.destroyAllWindows()
def plasma_detection(Position):
plasmaFrame=[]
for frame in range(nFrames):
img=plt.imread('Camera_' + Position + f'/Frames/{frame}.png')
rad=np.sum(img)
if rad>100:
plasmaFrame.append(frame)
plasma_start=float(requests.get(f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/PlasmaDetection/Results/t_plasma_start').text)
plasma_end=float(requests.get(f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/PlasmaDetection/Results/t_plasma_end').text)
FPMS = FPS[f'Camera {Position}']*1e-3
fStart = int(plasma_start*FPMS);fEnd=int(plasma_end*FPMS)+1
if len(plasmaFrame)==0:
print(f'Plasma detection via camera {Position} failed')
return fStart,fEnd
if abs(fStart-plasmaFrame[0])>2: fStart=plasmaFrame[0]; print('plasma start detected via camera')
if abs(fEnd-plasmaFrame[-1])>2: fEnd=plasmaFrame[-1]; print('plasma end detected via camera')
return fStart, fEnd
The plasma position is determined by calculating the center of mass of the signal: $$ \Delta_r[t_i] = \frac{\sum\limits_j^{1280}sig[j]\cdot r_j}{\sum\limits_j^{1280}sig[j]} $$ Where the previously calculated geometry information $r_j$ is employed to accurately map pixels to the central position of the chamber.\ Note that the geometry data assumes a fixed frame height (width for raw images) of 1280 pixels.
def plasma_position(shot_no, plasma_start, plasma_end, Position,symb, cams_mm):
duration = float(plasma_end)-float(plasma_start)
data=np.sum(plt.imread('plasma-film_'+Position+'.png'),axis=2)
r = [np.sum(data[:,i]*cams_mm)/np.sum(data[:,i]) for i in range(data.shape[1])]
camera_time = np.linspace(float(plasma_start), float(plasma_end), len(r))
plasma_position_data = pd.Series(r, index = camera_time)
return plasma_position_data
def PositionAndImg(camsRad_centerx, camsVert_centerz):
if Radial and Vertical: Position = ['Radial', 'Vertical']; nRow = 2
elif Radial: Position = ['Radial']; symb = ['r'];nRow = 1
elif Vertical: Position = ['Vertical'];nRow = 1
fig, ax = plt.subplots(nRow, 1, figsize = (10,9), sharex = True)
if nRow == 1: ax = [ax]
PositionDict = {'Radial': 'r', 'Vertical': 'v'}
for i, (Position, symb), cams_mm in zip(range(nRow), PositionDict.items(), [camsRad_centerx, camsVert_centerz]):
fStart, fEnd = plasma_detection(Position)
t_frame = 1/(FPS[f'Camera {Position}']*1e-3)
CamStart = fStart * t_frame; CamEnd = fEnd * t_frame
make_image(fStart, fEnd, Position)
plasma_position_camera = plasma_position(shot_no, CamStart, CamEnd, Position, symb, cams_mm)
FONT = 'DejaVu Sans'
img = Image.open('plasma-film_'+Position+'.png')
# number of time steps (in figure) covered by one frame
# TODO: should be determined by FPS,number of frames (~ discharge duration);
# tak aby i pri ruzne FPS, ale stejne delce vyboje, byl obrazek stejne velky (+neprilis maly, neprilis velky)
k = 8
newsize = (img.size[0]*k, img.size[1]) #resize image
print('newsize',newsize)
img = img.resize(newsize)
cv2.imwrite('plasma-film_' + Position + '2.png',np.asarray(img))
plasma_position_camera.plot(ax = ax[i], label = 'Fast Camera: '+ Position, lw=1)
ax[i].set_ylabel('$\Delta$'+symb+' [mm]',fontname = FONT, fontweight = 'bold', fontsize = 11)
ax[i].set_ylim(-95,95)
loclegend='best'
leg = ax[i].legend(loc = loclegend, shadow = True, fancybox=False)
leg.get_frame().set_linewidth(1)
leg.get_frame().set_edgecolor('k')
for text in leg.get_texts():
plt.setp(text, fontname=FONT, fontsize = 8)
for line, text in zip(leg.get_lines(), leg.get_texts()):
text.set_color(line.get_color())
plt.xticks(fontname=FONT, fontweight = 'bold', fontsize = 10)
plt.yticks(fontname=FONT, fontweight = 'bold', fontsize = 10)
ax[i].tick_params(which = 'major', direction='out', length=6, width=1.5)
ax[i].tick_params(which = 'minor', direction='out', length=3, width=1)
ax[i].spines['top'].set_visible(False)
ax[i].spines['right'].set_visible(False)
for axis in ['bottom','left']:
ax[i].spines[axis].set_linewidth(1.5)
ax[i].xaxis.set_minor_locator(AutoMinorLocator(2))
ax[i].yaxis.set_minor_locator(AutoMinorLocator(2))
ax[i].grid(which = 'major', c = 'gray', linewidth = 0.5, linestyle = 'solid')
ax[i].grid(which = 'minor', c = 'gray', linewidth = 0.3, linestyle = 'dashed')
for j in [-85,85,0]: ax[i].axhline(y=j, color='k', ls='--', lw=1, alpha=0.4)
fig.savefig('icon-fig.png')
# Save Data
savedata = 'Camera_'+Position+'/Camera'+Position+'Position'
plasma_position_camera.to_csv(savedata)
# Compute error -> TODO: zpresnit
FOV = abs(cams_mm[-1] - cams_mm[0])
PxLength = (FOV/cams_mm.size) # mm
err = PxLength*10
ax[i].fill_between(plasma_position_camera.index, plasma_position_camera - err, plasma_position_camera + err, alpha = 0.3)
ax[-1].set_xlabel('Time [ms]',fontname = FONT, fontweight = 'bold', fontsize = 11)
PositionAndImg(camsRad_centerx*1e3, camsVert_centerz*1e3) #center coords of chords in [mm]
plasma start detected via camera plasma end detected via camera Height: 1280 Width: 1 img_all shape (1280, 269, 3) newsize (2152, 1280) plasma start detected via camera plasma end detected via camera Height: 1280 Width: 1 img_all shape (1280, 267, 3) newsize (2136, 1280)
def icon_fig():
vert=plt.imread('plasma-film_Vertical2.png')
rad=plt.imread('plasma-film_Radial2.png')
maxlen=min(vert.shape[1],rad.shape[1])
stacked=np.vstack((rad[:,:maxlen,:],vert[:,:maxlen,:]))
img = cv2.convertScaleAbs(stacked, alpha=(255.0))
font = cv2.FONT_HERSHEY_SIMPLEX; org = [(int(rad.shape[0]*0.15), int(vert.shape[0]*0.15)),((150, int(rad.shape[0]*1.1)))]; fontScale = 4;color = (0, 0,255); thickness = 2
img = cv2.putText(img, 'Radial Position', org[0], font, fontScale, color, thickness, cv2.LINE_AA)
img = cv2.putText(img, 'Vertical Position', org[1], font, fontScale, color, thickness, cv2.LINE_AA)
y_org = [int(rad.shape[0]*0.15), int(rad.shape[0]*0.95), int(rad.shape[0]*1.1), int((rad.shape[0]+vert.shape[0])*0.95)]; name = ['LFS','HFS','top', 'bottom']
for n, y in zip(name, y_org):
img = cv2.putText(img, n, (maxlen-600, y), font, fontScale, color, thickness, cv2.LINE_AA)
cv2.imwrite('ScreenShotAll.png',img)
if Vertical and Radial:
icon_fig()
# show resulting image
fig, ax = plt.subplots(dpi=200)
im = plt.imread('ScreenShotAll.png')
ax.axhline(np.argmin(abs(camsRad_centerx)), c='tab:red', alpha = 0.6)
ax.axhline(np.argmin(abs(camsVert_centerz))+np.shape(im)[0]//2, c='tab:red', alpha = 0.6)
plt.imshow(im)