Photron - Data processing

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)

Check if the cameras are available

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

Load configuration file

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']

Geometry data

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')  

Remapping: px $\xrightarrow{}$ mm

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'])])

Visualization of the LOS

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)

Source: .mp4

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()

Plasma detection

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
    
    # disable camera plasma detection for now
    #if len(plasmaFrame)==0: 
    #    print(f'Plasma detection via camera {Position} failed')
    
    if fEnd >= nFrames:
        fEnd = nFrames-1
    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

Plasma Position

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]

Height: 1280 
Width: 1 
img_all shape (1280, 1004, 3)
newsize (8032, 1280)

<ipython-input-1-349bdac94605>:5: RuntimeWarning: invalid value encountered in scalar divide
  r = [np.sum(data[:,i]*cams_mm)/np.sum(data[:,i]) for i in range(data.shape[1])]

Height: 1280 
Width: 1 
img_all shape (1280, 1004, 3)
newsize (8032, 1280)

<ipython-input-1-349bdac94605>:5: RuntimeWarning: invalid value encountered in scalar divide
  r = [np.sum(data[:,i]*cams_mm)/np.sum(data[:,i]) for i in range(data.shape[1])]

Data:

• Radial plasma position
• Vertical plasma position

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)