diff --git a/README.md b/README.md index bb624da..aff4724 100644 --- a/README.md +++ b/README.md @@ -148,25 +148,28 @@ First, clone this repo! In /src you will find: -* PySpectrometer2-Picam2-v1.0.py (PySpectrometer for Raspberry Pi) -* PySpectrometer2-USB-v1.0.py (USB version of this program (This is for USB Cameras See end of Readme)). -* specFunctions.py (A library of functions including: Wavelength to RGB, SavGol filter from Scipy, Peak detect from peakutils, readcal and writecal. +* PySpectrometer2.py (The main program. See USB camera usage at the end of the Readme) +* specFunctions.py (A library of functions including: Wavelength to RGB, SavGol filter from Scipy, Peak detect from peakutils, readcal and writecal) ## Dependencies Run: **sudo apt-get install python3-opencv** **Also note, this build is designed for Raspberry Pi OS Bullseye, and will only work with the new libcamera based python library (picamera2)** -It will **not** work with older versions of Raspbery Pi OS. You **will** however be able to use PySpectrometer2-USB-v1.0.py with an external USB camera with other Operating Systems. +It will **not** work with older versions of Raspbery Pi OS. You **will** however be able to use with an external USB camera with other Operating Systems. -To run the program, first make it executable by running: **chmod +x PySpectrometer2-Picam2-v1.0.py** +To run the program, first make it executable by running: **chmod +x PySpectrometer2.py** -Run by typing: **./PySpectrometer2-Picam2-v1.0.py** +Run by typing: **./PySpectrometer2.py --X** +Where --X is one of: +* --picam (if using a picamera) +* --webcam (if using the first available camera on your computer) +* --device NUMBER (to specify the device) -Note to also display the waterfall display, run with the option: **./PySpectrometer2-Picam2-v1.0.py --waterfall** +Note to also display the waterfall display, run with the option: **--waterfall** -To run in fullscreen mode (perform calibration in standard mode first), run with the option: **./PySpectrometer2-Picam2-v1.0.py --fullscreen** +To run in fullscreen mode (perform calibration in standard mode first), run with the option: **--fullscreen** When first started, the spectrometer is in an uncalibrated state! You must therefore perform the calibration procedure, but at this stage you should be able to focus and align your camera with your spectroscope using the preview window. Is is expected that red light is on the right, and blue-violet on the left. An excellent choice for this is daylight, as well defined Fraunhoffer lines are indicative of good camera focus. @@ -259,7 +262,7 @@ For an external USB camera, first find the device by issuing: Once you have determined this, you can run the program. For example if your camera is /dev/video3 and you require a framerate of 15fps you would issue: -**./PySpectrometer2-USB-v1.0.py --device 3 --fps 15** +**./PySpectrometer2.py --device 3 --fps 15** If you want fine control over camera settings use guvcview: **sudo apt-get install guvcview** diff --git a/src/PySpectrometer2-Picam2-v1.0.py b/src/PySpectrometer2-Picam2-v1.0.py deleted file mode 100644 index c0965ef..0000000 --- a/src/PySpectrometer2-Picam2-v1.0.py +++ /dev/null @@ -1,465 +0,0 @@ -#!/usr/bin/env python3 - -''' -PySpectrometer2 Les Wright 2022 -https://www.youtube.com/leslaboratory -https://github.com/leswright1977 - -This project is a follow on from: https://github.com/leswright1977/PySpectrometer - -This is a more advanced, but more flexible version of the original program. Tk Has been dropped as the GUI to allow fullscreen mode on Raspberry Pi systems and the iterface is designed to fit 800*480 screens, which seem to be a common resolutin for RPi LCD's, paving the way for the creation of a stand alone benchtop instrument. - -Whats new: -Higher resolution (800px wide graph) -3 row pixel averaging of sensor data -Fullscreen option for the Spectrometer graph -3rd order polymonial fit of calibration data for accurate measurement. -Improved graph labelling -Labelled measurement cursors -Optional waterfall display for recording spectra changes over time. -Key Bindings for all operations - -All old features have been kept, including peak hold, peak detect, Savitsky Golay filter, and the ability to save graphs as png and data as CSV. - -For instructions please consult the readme! -''' - - -import cv2 -import time -import numpy as np -from specFunctions import wavelength_to_rgb,savitzky_golay,peakIndexes,readcal,writecal,background,generateGraticule -import base64 -import argparse -from picamera2 import Picamera2 - -parser = argparse.ArgumentParser() -group = parser.add_mutually_exclusive_group() -group.add_argument("--fullscreen", help="Fullscreen (Native 800*480)",action="store_true") -group.add_argument("--waterfall", help="Enable Waterfall (Windowed only)",action="store_true") -args = parser.parse_args() -dispFullscreen = False -dispWaterfall = False -if args.fullscreen: - print("Fullscreen Spectrometer enabled") - dispFullscreen = True -if args.waterfall: - print("Waterfall display enabled") - dispWaterfall = True - - - -frameWidth = 800 -frameHeight = 600 - -picam2 = Picamera2() -#need to spend more time at: https://datasheets.raspberrypi.com/camera/picamera2-manual.pdf -#but this will do for now! -#min and max microseconds per frame gives framerate. -#30fps (33333, 33333) -#25fps (40000, 40000) - -picamGain = 10.0 - -video_config = picam2.create_video_configuration(main={"format": 'RGB888', "size": (frameWidth, frameHeight)}, controls={"FrameDurationLimits": (33333, 33333)}) -picam2.configure(video_config) -picam2.start() - -#Change analog gain -#picam2.set_controls({"AnalogueGain": 10.0}) #Default 1 -#picam2.set_controls({"Brightness": 0.2}) #Default 0 range -1.0 to +1.0 -#picam2.set_controls({"Contrast": 1.8}) #Default 1 range 0.0-32.0 - - - -title1 = 'PySpectrometer 2 - Spectrograph' -title2 = 'PySpectrometer 2 - Waterfall' -stackHeight = 320+80+80 #height of the displayed CV window (graph+preview+messages) - -if dispWaterfall == True: - #watefall first so spectrum is on top - cv2.namedWindow(title2,cv2.WINDOW_GUI_NORMAL) - cv2.resizeWindow(title2,frameWidth,stackHeight) - cv2.moveWindow(title2,200,200); - -if dispFullscreen == True: - cv2.namedWindow(title1,cv2.WND_PROP_FULLSCREEN) - cv2.setWindowProperty(title1,cv2.WND_PROP_FULLSCREEN,cv2.WINDOW_FULLSCREEN) -else: - cv2.namedWindow(title1,cv2.WINDOW_GUI_NORMAL) - cv2.resizeWindow(title1,frameWidth,stackHeight) - cv2.moveWindow(title1,0,0); - -#settings for peak detect -savpoly = 7 #savgol filter polynomial max val 15 -mindist = 50 #minumum distance between peaks max val 100 -thresh = 20 #Threshold max val 100 - -calibrate = False - -clickArray = [] -cursorX = 0 -cursorY = 0 -def handle_mouse(event,x,y,flags,param): - global clickArray - global cursorX - global cursorY - mouseYOffset = 160 - if event == cv2.EVENT_MOUSEMOVE: - cursorX = x - cursorY = y - if event == cv2.EVENT_LBUTTONDOWN: - mouseX = x - mouseY = y-mouseYOffset - clickArray.append([mouseX,mouseY]) -#listen for click on plot window -cv2.setMouseCallback(title1,handle_mouse) - - -font=cv2.FONT_HERSHEY_SIMPLEX - -intensity = [0] * frameWidth #array for intensity data...full of zeroes - -holdpeaks = False #are we holding peaks? -measure = False #are we measuring? -recPixels = False #are we measuring pixels and recording clicks? - - -#messages -msg1 = "" -saveMsg = "No data saved" - -#blank image for Waterfall -waterfall = np.zeros([320,frameWidth,3],dtype=np.uint8) -waterfall.fill(0) #fill black - -#Go grab the computed calibration data -caldata = readcal(frameWidth) -wavelengthData = caldata[0] -calmsg1 = caldata[1] -calmsg2 = caldata[2] -calmsg3 = caldata[3] - -#generate the craticule data -graticuleData = generateGraticule(wavelengthData) -tens = (graticuleData[0]) -fifties = (graticuleData[1]) - -def snapshot(savedata): - now = time.strftime("%Y%m%d--%H%M%S") - timenow = time.strftime("%H:%M:%S") - imdata1 = savedata[0] - graphdata = savedata[1] - if dispWaterfall == True: - imdata2 = savedata[2] - cv2.imwrite("waterfall-" + now + ".png",imdata2) - cv2.imwrite("spectrum-" + now + ".png",imdata1) - #print(graphdata[0]) #wavelengths - #print(graphdata[1]) #intensities - f = open("Spectrum-"+now+'.csv','w') - f.write('Wavelength,Intensity\r\n') - for x in zip(graphdata[0],graphdata[1]): - f.write(str(x[0])+','+str(x[1])+'\r\n') - f.close() - message = "Last Save: "+timenow - return(message) - - -while True: - # Capture frame-by-frame - frame = picam2.capture_array() - y=int((frameHeight/2)-40) #origin of the vertical crop - #y=200 #origin of the vert crop - x=0 #origin of the horiz crop - h=80 #height of the crop - w=frameWidth #width of the crop - cropped = frame[y:y+h, x:x+w] - bwimage = cv2.cvtColor(cropped,cv2.COLOR_BGR2GRAY) - rows,cols = bwimage.shape - halfway =int(rows/2) - #show our line on the original image - #now a 3px wide region - cv2.line(cropped,(0,halfway-2),(frameWidth,halfway-2),(255,255,255),1) - cv2.line(cropped,(0,halfway+2),(frameWidth,halfway+2),(255,255,255),1) - - #banner image - decoded_data = base64.b64decode(background) - np_data = np.frombuffer(decoded_data,np.uint8) - img = cv2.imdecode(np_data,3) - messages = img - - #blank image for Graph - graph = np.zeros([320,frameWidth,3],dtype=np.uint8) - graph.fill(255) #fill white - - #Display a graticule calibrated with cal data - textoffset = 12 - #vertial lines every whole 10nm - for position in tens: - cv2.line(graph,(position,15),(position,320),(200,200,200),1) - - #vertical lines every whole 50nm - for positiondata in fifties: - cv2.line(graph,(positiondata[0],15),(positiondata[0],320),(0,0,0),1) - cv2.putText(graph,str(positiondata[1])+'nm',(positiondata[0]-textoffset,12),font,0.4,(0,0,0),1, cv2.LINE_AA) - - #horizontal lines - for i in range (320): - if i>=64: - if i%64==0: #suppress the first line then draw the rest... - cv2.line(graph,(0,i),(frameWidth,i),(100,100,100),1) - - #Now process the intensity data and display it - #intensity = [] - for i in range(cols): - #data = bwimage[halfway,i] #pull the pixel data from the halfway mark - #print(type(data)) #numpy.uint8 - #average the data of 3 rows of pixels: - dataminus1 = bwimage[halfway-1,i] - datazero = bwimage[halfway,i] #pull the pixel data from the halfway mark - dataplus1 = bwimage[halfway+1,i] - data = (int(dataminus1)+int(datazero)+int(dataplus1))/3 - data = np.uint8(data) - - - if holdpeaks == True: - if data > intensity[i]: - intensity[i] = data - else: - intensity[i] = data - - if dispWaterfall == True: - #waterfall.... - #data is smoothed at this point!!!!!! - #create an empty array for the data - wdata = np.zeros([1,frameWidth,3],dtype=np.uint8) - index=0 - for i in intensity: - rgb = wavelength_to_rgb(round(wavelengthData[index]))#derive the color from the wavelenthData array - luminosity = intensity[index]/255 - b = int(round(rgb[0]*luminosity)) - g = int(round(rgb[1]*luminosity)) - r = int(round(rgb[2]*luminosity)) - #print(b,g,r) - #wdata[0,index]=(r,g,b) #fix me!!! how do we deal with this data?? - wdata[0,index]=(r,g,b) - index+=1 - #bright and contrast of final image - contrast = 2.5 - brightness =10 - wdata = cv2.addWeighted( wdata, contrast, wdata, 0, brightness) - waterfall = np.insert(waterfall, 0, wdata, axis=0) #insert line to beginning of array - waterfall = waterfall[:-1].copy() #remove last element from array - - hsv = cv2.cvtColor(waterfall, cv2.COLOR_BGR2HSV) - - - #Draw the intensity data :-) - #first filter if not holding peaks! - - if holdpeaks == False: - intensity = savitzky_golay(intensity,17,savpoly) - intensity = np.array(intensity) - intensity = intensity.astype(int) - holdmsg = "Holdpeaks OFF" - else: - holdmsg = "Holdpeaks ON" - - - #now draw the intensity data.... - index=0 - for i in intensity: - rgb = wavelength_to_rgb(round(wavelengthData[index]))#derive the color from the wvalenthData array - r = rgb[0] - g = rgb[1] - b = rgb[2] - #or some reason origin is top left. - cv2.line(graph, (index,320), (index,320-i), (b,g,r), 1) - cv2.line(graph, (index,319-i), (index,320-i), (0,0,0), 1,cv2.LINE_AA) - index+=1 - - - #find peaks and label them - textoffset = 12 - thresh = int(thresh) #make sure the data is int. - indexes = peakIndexes(intensity, thres=thresh/max(intensity), min_dist=mindist) - #print(indexes) - for i in indexes: - height = intensity[i] - height = 310-height - wavelength = round(wavelengthData[i],1) - cv2.rectangle(graph,((i-textoffset)-2,height),((i-textoffset)+60,height-15),(0,255,255),-1) - cv2.rectangle(graph,((i-textoffset)-2,height),((i-textoffset)+60,height-15),(0,0,0),1) - cv2.putText(graph,str(wavelength)+'nm',(i-textoffset,height-3),font,0.4,(0,0,0),1, cv2.LINE_AA) - #flagpoles - cv2.line(graph,(i,height),(i,height+10),(0,0,0),1) - - - if measure == True: - #show the cursor! - cv2.line(graph,(cursorX,cursorY-140),(cursorX,cursorY-180),(0,0,0),1) - cv2.line(graph,(cursorX-20,cursorY-160),(cursorX+20,cursorY-160),(0,0,0),1) - cv2.putText(graph,str(round(wavelengthData[cursorX],2))+'nm',(cursorX+5,cursorY-165),font,0.4,(0,0,0),1, cv2.LINE_AA) - - if recPixels == True: - #display the points - cv2.line(graph,(cursorX,cursorY-140),(cursorX,cursorY-180),(0,0,0),1) - cv2.line(graph,(cursorX-20,cursorY-160),(cursorX+20,cursorY-160),(0,0,0),1) - cv2.putText(graph,str(cursorX)+'px',(cursorX+5,cursorY-165),font,0.4,(0,0,0),1, cv2.LINE_AA) - else: - #also make sure the click array stays empty - clickArray = [] - - if clickArray: - for data in clickArray: - mouseX=data[0] - mouseY=data[1] - cv2.circle(graph,(mouseX,mouseY),5,(0,0,0),-1) - #we can display text :-) so we can work out wavelength from x-pos and display it ultimately - cv2.putText(graph,str(mouseX),(mouseX+5,mouseY),cv2.FONT_HERSHEY_SIMPLEX,0.4,(0,0,0)) - - - - - #stack the images and display the spectrum - spectrum_vertical = np.vstack((messages,cropped, graph)) - #dividing lines... - cv2.line(spectrum_vertical,(0,80),(frameWidth,80),(255,255,255),1) - cv2.line(spectrum_vertical,(0,160),(frameWidth,160),(255,255,255),1) - #print the messages - cv2.putText(spectrum_vertical,calmsg1,(490,15),font,0.4,(0,255,255),1, cv2.LINE_AA) - cv2.putText(spectrum_vertical,calmsg3,(490,33),font,0.4,(0,255,255),1, cv2.LINE_AA) - cv2.putText(spectrum_vertical,saveMsg,(490,51),font,0.4,(0,255,255),1, cv2.LINE_AA) - cv2.putText(spectrum_vertical,"Gain: "+str(picamGain),(490,69),font,0.4,(0,255,255),1, cv2.LINE_AA) - #Second column - cv2.putText(spectrum_vertical,holdmsg,(640,15),font,0.4,(0,255,255),1, cv2.LINE_AA) - cv2.putText(spectrum_vertical,"Savgol Filter: "+str(savpoly),(640,33),font,0.4,(0,255,255),1, cv2.LINE_AA) - cv2.putText(spectrum_vertical,"Label Peak Width: "+str(mindist),(640,51),font,0.4,(0,255,255),1, cv2.LINE_AA) - cv2.putText(spectrum_vertical,"Label Threshold: "+str(thresh),(640,69),font,0.4,(0,255,255),1, cv2.LINE_AA) - cv2.imshow(title1,spectrum_vertical) - - if dispWaterfall == True: - #stack the images and display the waterfall - waterfall_vertical = np.vstack((messages,cropped, waterfall)) - #dividing lines... - cv2.line(waterfall_vertical,(0,80),(frameWidth,80),(255,255,255),1) - cv2.line(waterfall_vertical,(0,160),(frameWidth,160),(255,255,255),1) - #Draw this stuff over the top of the image! - #Display a graticule calibrated with cal data - textoffset = 12 - - #vertical lines every whole 50nm - for positiondata in fifties: - for i in range(162,480): - if i%20 == 0: - cv2.line(waterfall_vertical,(positiondata[0],i),(positiondata[0],i+1),(0,0,0),2) - cv2.line(waterfall_vertical,(positiondata[0],i),(positiondata[0],i+1),(255,255,255),1) - cv2.putText(waterfall_vertical,str(positiondata[1])+'nm',(positiondata[0]-textoffset,475),font,0.4,(0,0,0),2, cv2.LINE_AA) - cv2.putText(waterfall_vertical,str(positiondata[1])+'nm',(positiondata[0]-textoffset,475),font,0.4,(255,255,255),1, cv2.LINE_AA) - - cv2.putText(waterfall_vertical,calmsg1,(490,15),font,0.4,(0,255,255),1, cv2.LINE_AA) - cv2.putText(waterfall_vertical,calmsg3,(490,33),font,0.4,(0,255,255),1, cv2.LINE_AA) - cv2.putText(waterfall_vertical,saveMsg,(490,51),font,0.4,(0,255,255),1, cv2.LINE_AA) - cv2.putText(waterfall_vertical,"Gain: "+str(picamGain),(490,69),font,0.4,(0,255,255),1, cv2.LINE_AA) - - cv2.putText(waterfall_vertical,holdmsg,(640,15),font,0.4,(0,255,255),1, cv2.LINE_AA) - - cv2.imshow(title2,waterfall_vertical) - - - keyPress = cv2.waitKey(1) - if keyPress == ord('q'): - break - elif keyPress == ord('h'): - if holdpeaks == False: - holdpeaks = True - elif holdpeaks == True: - holdpeaks = False - elif keyPress == ord("s"): - #package up the data! - graphdata = [] - graphdata.append(wavelengthData) - graphdata.append(intensity) - if dispWaterfall == True: - savedata = [] - savedata.append(spectrum_vertical) - savedata.append(graphdata) - savedata.append(waterfall_vertical) - else: - savedata = [] - savedata.append(spectrum_vertical) - savedata.append(graphdata) - saveMsg = snapshot(savedata) - elif keyPress == ord("c"): - calcomplete = writecal(clickArray) - if calcomplete: - #overwrite wavelength data - #Go grab the computed calibration data - caldata = readcal(frameWidth) - wavelengthData = caldata[0] - calmsg1 = caldata[1] - calmsg2 = caldata[2] - calmsg3 = caldata[3] - #overwrite graticule data - graticuleData = generateGraticule(wavelengthData) - tens = (graticuleData[0]) - fifties = (graticuleData[1]) - elif keyPress == ord("x"): - clickArray = [] - elif keyPress == ord("m"): - recPixels = False #turn off recpixels! - if measure == False: - measure = True - elif measure == True: - measure = False - elif keyPress == ord("p"): - measure = False #turn off measure! - if recPixels == False: - recPixels = True - elif recPixels == True: - recPixels = False - elif keyPress == ord("o"):#sav up - savpoly+=1 - if savpoly >=15: - savpoly=15 - elif keyPress == ord("l"):#sav down - savpoly-=1 - if savpoly <=0: - savpoly=0 - elif keyPress == ord("i"):#Peak width up - mindist+=1 - if mindist >=100: - mindist=100 - elif keyPress == ord("k"):#Peak Width down - mindist-=1 - if mindist <=0: - mindist=0 - elif keyPress == ord("u"):#label thresh up - thresh+=1 - if thresh >=100: - thresh=100 - elif keyPress == ord("j"):#label thresh down - thresh-=1 - if thresh <=0: - thresh=0 - - elif keyPress == ord("t"):#Gain up! - picamGain += 1 - if picamGain >=50: - picamGain = 50.0 - picam2.set_controls({"AnalogueGain": picamGain}) - print("Camera Gain: "+str(picamGain)) - elif keyPress == ord("g"):#Gain down - picamGain -= 1 - if picamGain <=0: - picamGain = 0.0 - picam2.set_controls({"AnalogueGain": picamGain}) - print("Camera Gain: "+str(picamGain)) - - - - -#Everything done -cv2.destroyAllWindows() - - diff --git a/src/PySpectrometer2-USB-v1.0.py b/src/PySpectrometer2-USB-v1.0.py deleted file mode 100644 index 8bb6adc..0000000 --- a/src/PySpectrometer2-USB-v1.0.py +++ /dev/null @@ -1,456 +0,0 @@ -#!/usr/bin/env python3 - -''' -PySpectrometer2 Les Wright 2022 -https://www.youtube.com/leslaboratory -https://github.com/leswright1977 - -This project is a follow on from: https://github.com/leswright1977/PySpectrometer - -This is a more advanced, but more flexible version of the original program. Tk Has been dropped as the GUI to allow fullscreen mode on Raspberry Pi systems and the iterface is designed to fit 800*480 screens, which seem to be a common resolutin for RPi LCD's, paving the way for the creation of a stand alone benchtop instrument. - -Whats new: -Higher resolution (800px wide graph) -3 row pixel averaging of sensor data -Fullscreen option for the Spectrometer graph -3rd order polymonial fit of calibration data for accurate measurement. -Improved graph labelling -Labelled measurement cursors -Optional waterfall display for recording spectra changes over time. -Key Bindings for all operations - -All old features have been kept, including peak hold, peak detect, Savitsky Golay filter, and the ability to save graphs as png and data as CSV. - -For instructions please consult the readme! - -''' - - -import cv2 -import time -import numpy as np -from specFunctions import wavelength_to_rgb,savitzky_golay,peakIndexes,readcal,writecal,background,generateGraticule -import base64 -import argparse - -parser = argparse.ArgumentParser() -parser.add_argument("--device", type=int, default=0, help="Video Device number e.g. 0, use v4l2-ctl --list-devices") -parser.add_argument("--fps", type=int, default=30, help="Frame Rate e.g. 30") -group = parser.add_mutually_exclusive_group() -group.add_argument("--fullscreen", help="Fullscreen (Native 800*480)",action="store_true") -group.add_argument("--waterfall", help="Enable Waterfall (Windowed only)",action="store_true") -args = parser.parse_args() -dispFullscreen = False -dispWaterfall = False -if args.fullscreen: - print("Fullscreen Spectrometer enabled") - dispFullscreen = True -if args.waterfall: - print("Waterfall display enabled") - dispWaterfall = True - -if args.device: - dev = args.device -else: - dev = 0 - -if args.fps: - fps = args.fps -else: - fps = 30 - -frameWidth = 800 -frameHeight = 600 - - -#init video -cap = cv2.VideoCapture('/dev/video'+str(dev), cv2.CAP_V4L) -#cap = cv2.VideoCapture(0) -print("[info] W, H, FPS") -cap.set(cv2.CAP_PROP_FRAME_WIDTH,frameWidth) -cap.set(cv2.CAP_PROP_FRAME_HEIGHT,frameHeight) -cap.set(cv2.CAP_PROP_FPS,fps) -print(cap.get(cv2.CAP_PROP_FRAME_WIDTH)) -print(cap.get(cv2.CAP_PROP_FRAME_HEIGHT)) -print(cap.get(cv2.CAP_PROP_FPS)) -cfps = (cap.get(cv2.CAP_PROP_FPS)) - - -title1 = 'PySpectrometer 2 - Spectrograph' -title2 = 'PySpectrometer 2 - Waterfall' -stackHeight = 320+80+80 #height of the displayed CV window (graph+preview+messages) - -if dispWaterfall == True: - #watefall first so spectrum is on top - cv2.namedWindow(title2,cv2.WINDOW_GUI_NORMAL) - cv2.resizeWindow(title2,frameWidth,stackHeight) - cv2.moveWindow(title2,200,200); - -if dispFullscreen == True: - cv2.namedWindow(title1,cv2.WND_PROP_FULLSCREEN) - cv2.setWindowProperty(title1,cv2.WND_PROP_FULLSCREEN,cv2.WINDOW_FULLSCREEN) -else: - cv2.namedWindow(title1,cv2.WINDOW_GUI_NORMAL) - cv2.resizeWindow(title1,frameWidth,stackHeight) - cv2.moveWindow(title1,0,0); - -#settings for peak detect -savpoly = 7 #savgol filter polynomial max val 15 -mindist = 50 #minumum distance between peaks max val 100 -thresh = 20 #Threshold max val 100 - -calibrate = False - -clickArray = [] -cursorX = 0 -cursorY = 0 -def handle_mouse(event,x,y,flags,param): - global clickArray - global cursorX - global cursorY - mouseYOffset = 160 - if event == cv2.EVENT_MOUSEMOVE: - cursorX = x - cursorY = y - if event == cv2.EVENT_LBUTTONDOWN: - mouseX = x - mouseY = y-mouseYOffset - clickArray.append([mouseX,mouseY]) -#listen for click on plot window -cv2.setMouseCallback(title1,handle_mouse) - - -font=cv2.FONT_HERSHEY_SIMPLEX - -intensity = [0] * frameWidth #array for intensity data...full of zeroes - -holdpeaks = False #are we holding peaks? -measure = False #are we measuring? -recPixels = False #are we measuring pixels and recording clicks? - - -#messages -msg1 = "" -saveMsg = "No data saved" - -#blank image for Waterfall -waterfall = np.zeros([320,frameWidth,3],dtype=np.uint8) -waterfall.fill(0) #fill black - -#Go grab the computed calibration data -caldata = readcal(frameWidth) -wavelengthData = caldata[0] -calmsg1 = caldata[1] -calmsg2 = caldata[2] -calmsg3 = caldata[3] - -#generate the craticule data -graticuleData = generateGraticule(wavelengthData) -tens = (graticuleData[0]) -fifties = (graticuleData[1]) - -def snapshot(savedata): - now = time.strftime("%Y%m%d--%H%M%S") - timenow = time.strftime("%H:%M:%S") - imdata1 = savedata[0] - graphdata = savedata[1] - if dispWaterfall == True: - imdata2 = savedata[2] - cv2.imwrite("waterfall-" + now + ".png",imdata2) - cv2.imwrite("spectrum-" + now + ".png",imdata1) - #print(graphdata[0]) #wavelengths - #print(graphdata[1]) #intensities - f = open("Spectrum-"+now+'.csv','w') - f.write('Wavelength,Intensity\r\n') - for x in zip(graphdata[0],graphdata[1]): - f.write(str(x[0])+','+str(x[1])+'\r\n') - f.close() - message = "Last Save: "+timenow - return(message) - -while(cap.isOpened()): - # Capture frame-by-frame - ret, frame = cap.read() - - if ret == True: - y=int((frameHeight/2)-40) #origin of the vertical crop - #y=200 #origin of the vert crop - x=0 #origin of the horiz crop - h=80 #height of the crop - w=frameWidth #width of the crop - cropped = frame[y:y+h, x:x+w] - bwimage = cv2.cvtColor(cropped,cv2.COLOR_BGR2GRAY) - rows,cols = bwimage.shape - halfway =int(rows/2) - #show our line on the original image - #now a 3px wide region - cv2.line(cropped,(0,halfway-2),(frameWidth,halfway-2),(255,255,255),1) - cv2.line(cropped,(0,halfway+2),(frameWidth,halfway+2),(255,255,255),1) - - #banner image - decoded_data = base64.b64decode(background) - np_data = np.frombuffer(decoded_data,np.uint8) - img = cv2.imdecode(np_data,3) - messages = img - - #blank image for Graph - graph = np.zeros([320,frameWidth,3],dtype=np.uint8) - graph.fill(255) #fill white - - #Display a graticule calibrated with cal data - textoffset = 12 - #vertial lines every whole 10nm - for position in tens: - cv2.line(graph,(position,15),(position,320),(200,200,200),1) - - #vertical lines every whole 50nm - for positiondata in fifties: - cv2.line(graph,(positiondata[0],15),(positiondata[0],320),(0,0,0),1) - cv2.putText(graph,str(positiondata[1])+'nm',(positiondata[0]-textoffset,12),font,0.4,(0,0,0),1, cv2.LINE_AA) - - #horizontal lines - for i in range (320): - if i>=64: - if i%64==0: #suppress the first line then draw the rest... - cv2.line(graph,(0,i),(frameWidth,i),(100,100,100),1) - - #Now process the intensity data and display it - #intensity = [] - for i in range(cols): - #data = bwimage[halfway,i] #pull the pixel data from the halfway mark - #print(type(data)) #numpy.uint8 - #average the data of 3 rows of pixels: - dataminus1 = bwimage[halfway-1,i] - datazero = bwimage[halfway,i] #pull the pixel data from the halfway mark - dataplus1 = bwimage[halfway+1,i] - data = (int(dataminus1)+int(datazero)+int(dataplus1))/3 - data = np.uint8(data) - - - if holdpeaks == True: - if data > intensity[i]: - intensity[i] = data - else: - intensity[i] = data - - if dispWaterfall == True: - #waterfall.... - #data is smoothed at this point!!!!!! - #create an empty array for the data - wdata = np.zeros([1,frameWidth,3],dtype=np.uint8) - index=0 - for i in intensity: - rgb = wavelength_to_rgb(round(wavelengthData[index]))#derive the color from the wavelenthData array - luminosity = intensity[index]/255 - b = int(round(rgb[0]*luminosity)) - g = int(round(rgb[1]*luminosity)) - r = int(round(rgb[2]*luminosity)) - #print(b,g,r) - #wdata[0,index]=(r,g,b) #fix me!!! how do we deal with this data?? - wdata[0,index]=(r,g,b) - index+=1 - waterfall = np.insert(waterfall, 0, wdata, axis=0) #insert line to beginning of array - waterfall = waterfall[:-1].copy() #remove last element from array - - hsv = cv2.cvtColor(waterfall, cv2.COLOR_BGR2HSV) - - - - #Draw the intensity data :-) - #first filter if not holding peaks! - - if holdpeaks == False: - intensity = savitzky_golay(intensity,17,savpoly) - intensity = np.array(intensity) - intensity = intensity.astype(int) - holdmsg = "Holdpeaks OFF" - else: - holdmsg = "Holdpeaks ON" - - - #now draw the intensity data.... - index=0 - for i in intensity: - rgb = wavelength_to_rgb(round(wavelengthData[index]))#derive the color from the wvalenthData array - r = rgb[0] - g = rgb[1] - b = rgb[2] - #or some reason origin is top left. - cv2.line(graph, (index,320), (index,320-i), (b,g,r), 1) - cv2.line(graph, (index,319-i), (index,320-i), (0,0,0), 1,cv2.LINE_AA) - index+=1 - - - #find peaks and label them - textoffset = 12 - thresh = int(thresh) #make sure the data is int. - indexes = peakIndexes(intensity, thres=thresh/max(intensity), min_dist=mindist) - #print(indexes) - for i in indexes: - height = intensity[i] - height = 310-height - wavelength = round(wavelengthData[i],1) - cv2.rectangle(graph,((i-textoffset)-2,height),((i-textoffset)+60,height-15),(0,255,255),-1) - cv2.rectangle(graph,((i-textoffset)-2,height),((i-textoffset)+60,height-15),(0,0,0),1) - cv2.putText(graph,str(wavelength)+'nm',(i-textoffset,height-3),font,0.4,(0,0,0),1, cv2.LINE_AA) - #flagpoles - cv2.line(graph,(i,height),(i,height+10),(0,0,0),1) - - - if measure == True: - #show the cursor! - cv2.line(graph,(cursorX,cursorY-140),(cursorX,cursorY-180),(0,0,0),1) - cv2.line(graph,(cursorX-20,cursorY-160),(cursorX+20,cursorY-160),(0,0,0),1) - cv2.putText(graph,str(round(wavelengthData[cursorX],2))+'nm',(cursorX+5,cursorY-165),font,0.4,(0,0,0),1, cv2.LINE_AA) - - if recPixels == True: - #display the points - cv2.line(graph,(cursorX,cursorY-140),(cursorX,cursorY-180),(0,0,0),1) - cv2.line(graph,(cursorX-20,cursorY-160),(cursorX+20,cursorY-160),(0,0,0),1) - cv2.putText(graph,str(cursorX)+'px',(cursorX+5,cursorY-165),font,0.4,(0,0,0),1, cv2.LINE_AA) - else: - #also make sure the click array stays empty - clickArray = [] - - if clickArray: - for data in clickArray: - mouseX=data[0] - mouseY=data[1] - cv2.circle(graph,(mouseX,mouseY),5,(0,0,0),-1) - #we can display text :-) so we can work out wavelength from x-pos and display it ultimately - cv2.putText(graph,str(mouseX),(mouseX+5,mouseY),cv2.FONT_HERSHEY_SIMPLEX,0.4,(0,0,0)) - - - - - #stack the images and display the spectrum - spectrum_vertical = np.vstack((messages,cropped, graph)) - #dividing lines... - cv2.line(spectrum_vertical,(0,80),(frameWidth,80),(255,255,255),1) - cv2.line(spectrum_vertical,(0,160),(frameWidth,160),(255,255,255),1) - #print the messages - cv2.putText(spectrum_vertical,calmsg1,(490,15),font,0.4,(0,255,255),1, cv2.LINE_AA) - cv2.putText(spectrum_vertical,calmsg3,(490,33),font,0.4,(0,255,255),1, cv2.LINE_AA) - cv2.putText(spectrum_vertical,"Framerate: "+str(cfps),(490,51),font,0.4,(0,255,255),1, cv2.LINE_AA) - cv2.putText(spectrum_vertical,saveMsg,(490,69),font,0.4,(0,255,255),1, cv2.LINE_AA) - #Second column - cv2.putText(spectrum_vertical,holdmsg,(640,15),font,0.4,(0,255,255),1, cv2.LINE_AA) - cv2.putText(spectrum_vertical,"Savgol Filter: "+str(savpoly),(640,33),font,0.4,(0,255,255),1, cv2.LINE_AA) - cv2.putText(spectrum_vertical,"Label Peak Width: "+str(mindist),(640,51),font,0.4,(0,255,255),1, cv2.LINE_AA) - cv2.putText(spectrum_vertical,"Label Threshold: "+str(thresh),(640,69),font,0.4,(0,255,255),1, cv2.LINE_AA) - cv2.imshow(title1,spectrum_vertical) - - if dispWaterfall == True: - #stack the images and display the waterfall - waterfall_vertical = np.vstack((messages,cropped, waterfall)) - #dividing lines... - cv2.line(waterfall_vertical,(0,80),(frameWidth,80),(255,255,255),1) - cv2.line(waterfall_vertical,(0,160),(frameWidth,160),(255,255,255),1) - #Draw this stuff over the top of the image! - #Display a graticule calibrated with cal data - textoffset = 12 - - #vertical lines every whole 50nm - for positiondata in fifties: - for i in range(162,480): - if i%20 == 0: - cv2.line(waterfall_vertical,(positiondata[0],i),(positiondata[0],i+1),(0,0,0),2) - cv2.line(waterfall_vertical,(positiondata[0],i),(positiondata[0],i+1),(255,255,255),1) - cv2.putText(waterfall_vertical,str(positiondata[1])+'nm',(positiondata[0]-textoffset,475),font,0.4,(0,0,0),2, cv2.LINE_AA) - cv2.putText(waterfall_vertical,str(positiondata[1])+'nm',(positiondata[0]-textoffset,475),font,0.4,(255,255,255),1, cv2.LINE_AA) - - cv2.putText(waterfall_vertical,calmsg1,(490,15),font,0.4,(0,255,255),1, cv2.LINE_AA) - cv2.putText(waterfall_vertical,calmsg2,(490,33),font,0.4,(0,255,255),1, cv2.LINE_AA) - cv2.putText(waterfall_vertical,calmsg3,(490,51),font,0.4,(0,255,255),1, cv2.LINE_AA) - cv2.putText(waterfall_vertical,saveMsg,(490,69),font,0.4,(0,255,255),1, cv2.LINE_AA) - - cv2.putText(waterfall_vertical,holdmsg,(640,15),font,0.4,(0,255,255),1, cv2.LINE_AA) - - cv2.imshow(title2,waterfall_vertical) - - - keyPress = cv2.waitKey(1) - if keyPress == ord('q'): - break - elif keyPress == ord('h'): - if holdpeaks == False: - holdpeaks = True - elif holdpeaks == True: - holdpeaks = False - elif keyPress == ord("s"): - #package up the data! - graphdata = [] - graphdata.append(wavelengthData) - graphdata.append(intensity) - if dispWaterfall == True: - savedata = [] - savedata.append(spectrum_vertical) - savedata.append(graphdata) - savedata.append(waterfall_vertical) - else: - savedata = [] - savedata.append(spectrum_vertical) - savedata.append(graphdata) - saveMsg = snapshot(savedata) - elif keyPress == ord("c"): - calcomplete = writecal(clickArray) - if calcomplete: - #overwrite wavelength data - #Go grab the computed calibration data - caldata = readcal(frameWidth) - wavelengthData = caldata[0] - calmsg1 = caldata[1] - calmsg2 = caldata[2] - calmsg3 = caldata[3] - #overwrite graticule data - graticuleData = generateGraticule(wavelengthData) - tens = (graticuleData[0]) - fifties = (graticuleData[1]) - elif keyPress == ord("x"): - clickArray = [] - elif keyPress == ord("m"): - recPixels = False #turn off recpixels! - if measure == False: - measure = True - elif measure == True: - measure = False - elif keyPress == ord("p"): - measure = False #turn off measure! - if recPixels == False: - recPixels = True - elif recPixels == True: - recPixels = False - elif keyPress == ord("o"):#sav up - savpoly+=1 - if savpoly >=15: - savpoly=15 - elif keyPress == ord("l"):#sav down - savpoly-=1 - if savpoly <=0: - savpoly=0 - elif keyPress == ord("i"):#Peak width up - mindist+=1 - if mindist >=100: - mindist=100 - elif keyPress == ord("k"):#Peak Width down - mindist-=1 - if mindist <=0: - mindist=0 - elif keyPress == ord("u"):#label thresh up - thresh+=1 - if thresh >=100: - thresh=100 - elif keyPress == ord("j"):#label thresh down - thresh-=1 - if thresh <=0: - thresh=0 - else: - break - - -#Everything done, release the vid -cap.release() - -cv2.destroyAllWindows() - - diff --git a/src/PySpectrometer2.py b/src/PySpectrometer2.py new file mode 100755 index 0000000..e82be4f --- /dev/null +++ b/src/PySpectrometer2.py @@ -0,0 +1,585 @@ +#!/usr/bin/env python3 + +''' +PySpectrometer2 Les Wright 2022 +https://www.youtube.com/leslaboratory +https://github.com/leswright1977 + +This project is a follow on from: https://github.com/leswright1977/PySpectrometer + +This is a more advanced, but more flexible version of the original program. Tk Has been dropped as the GUI to allow fullscreen mode on Raspberry Pi systems and the iterface is designed to fit 800*480 screens, which seem to be a common resolutin for RPi LCD's, paving the way for the creation of a stand alone benchtop instrument. + +Whats new: +Higher resolution (800px wide graph) +3 row pixel averaging of sensor data +Fullscreen option for the Spectrometer graph +3rd order polymonial fit of calibration data for accurate measurement. +Improved graph labelling +Labelled measurement cursors +Optional waterfall display for recording spectra changes over time. +Key Bindings for all operations + +All old features have been kept, including peak hold, peak detect, Savitsky Golay filter, and the ability to save graphs as png and data as CSV. + +For instructions please consult the readme! +''' + + +import cv2 +import time +import numpy as np +from specFunctions import wavelength_to_rgb,savitzky_golay,peakIndexes,readcal,writecal,background,generateGraticule,closest +import base64 +import argparse +from pprint import pprint + +parser = argparse.ArgumentParser() +parser.add_argument("--device", type=int, help="Video Device number e.g. 0, use v4l2-ctl --list-devices") +parser.add_argument("--webcam", action=argparse.BooleanOptionalAction, help="Use the first available camera") +parser.add_argument("--picam", action=argparse.BooleanOptionalAction, help="Use the picamera") +parser.add_argument("--fps", type=int, default=30, help="Frame Rate e.g. 30") +group = parser.add_mutually_exclusive_group() +group.add_argument("--fullscreen", help="Fullscreen (Native 800*480)",action="store_true") +group.add_argument("--waterfall", help="Enable Waterfall (Windowed only)",action="store_true") +args = parser.parse_args() +dispFullscreen = False +dispWaterfall = False +if args.fullscreen: + print("Fullscreen Spectrometer enabled") + dispFullscreen = True +if args.waterfall: + print("Waterfall display enabled") + dispWaterfall = True + +use_webcam = args.webcam is not None +use_device = args.device is not None +use_picamera = args.picam is not None + +fps = args.fps +dev = 0 +if args.device: + dev = args.device + +peak_intensities = [] +frameWidth = 800 +frameHeight = int(frameWidth * 0.75) +message_loc1 = frameWidth - 310 +message_loc2 = frameWidth - 160 + +picam2 = None +if use_picamera: + from libcamera import controls + from picamera2 import Picamera2 + + picam2 = Picamera2() + #need to spend more time at: https://datasheets.raspberrypi.com/camera/picamera2-manual.pdf + #but this will do for now! + #min and max microseconds per frame gives framerate. + #30fps (33333, 33333) + #25fps (40000, 40000) + + NoiseReductionMode = controls.draft.NoiseReductionModeEnum + picamGain = 10.0 + video_config = picam2.create_video_configuration( + main={ + "format": 'RGB888', + "size": (frameWidth, frameHeight) + }, + controls={ + "NoiseReductionMode": NoiseReductionMode.Fast, + "FrameDurationLimits": (33333, 33333), + "AnalogueGain": picamGain + }) + pprint(video_config["controls"]) + picam2.configure(video_config) + picam2.start() + + #Change analog gain + #picam2.set_controls({"AnalogueGain": 10.0}) #Default 1 + #picam2.set_controls({"Brightness": 0.2}) #Default 0 range -1.0 to +1.0 + #picam2.set_controls({"Contrast": 1.8}) #Default 1 range 0.0-32.0 + +cap = None +if use_device or use_webcam: + if use_device: + cap = cv2.VideoCapture("/dev/video{}".format(dev), cv2.CAP_V4L) + elif use_webcam: + cap = cv2.VideoCapture(0) + print("[info] W, H, FPS") + cap.set(cv2.CAP_PROP_FRAME_WIDTH,frameWidth) + cap.set(cv2.CAP_PROP_FRAME_HEIGHT,frameHeight) + cap.set(cv2.CAP_PROP_FPS,fps) + print(cap.get(cv2.CAP_PROP_FRAME_WIDTH)) + print(cap.get(cv2.CAP_PROP_FRAME_HEIGHT)) + print(cap.get(cv2.CAP_PROP_FPS)) + cfps = (cap.get(cv2.CAP_PROP_FPS)) + + +title1 = 'PySpectrometer 2 - Spectrograph' +title2 = 'PySpectrometer 2 - Waterfall' +stackHeight = 320+80+80 #height of the displayed CV window (graph+preview+messages) + +if dispWaterfall: + #watefall first so spectrum is on top + cv2.namedWindow(title2,cv2.WINDOW_GUI_NORMAL) + cv2.resizeWindow(title2,frameWidth,stackHeight) + cv2.moveWindow(title2,200,200); + +if dispFullscreen: + cv2.namedWindow(title1,cv2.WND_PROP_FULLSCREEN) + cv2.setWindowProperty(title1,cv2.WND_PROP_FULLSCREEN,cv2.WINDOW_FULLSCREEN) +else: + cv2.namedWindow(title1,cv2.WINDOW_GUI_NORMAL) + cv2.resizeWindow(title1,frameWidth,stackHeight) + cv2.moveWindow(title1,0,0); + +#settings for peak detect +savpoly = 7 #savgol filter polynomial max val 15 +mindist = 50 #minumum distance between peaks max val 100 +thresh = 20 #Threshold max val 100 + +calibrate = False + +clickArray = [] +cursorX = 0 +cursorY = 0 +def handle_mouse(event,x,y,flags,param): + global clickArray, cursorX, cursorY, peak_intensities, wavelengthData + mouseYOffset = 160 + if event == cv2.EVENT_MOUSEMOVE: + cursorX = x + cursorY = y + elif event == cv2.EVENT_LBUTTONDOWN: + mouseX = x + mouseY = y-mouseYOffset + if len(peak_intensities) > 0: + # use raw peaks + closest_peak_index = 0 + peak_distance = 5 + left = mouseX - peak_distance + right = mouseX + peak_distance + peak_index = 0 + peak_intensity = 0 + for index, i in enumerate(raw_intensity[left:right]): + if i > peak_intensity: + peak_intensity = i + peak_index = index+left + closest_peak_index = peak_index + + # # use labeled peaks + # closest_peak_index = closest(peak_intensities, mouseX) + + wl = wavelengthData[closest_peak_index] + clickArray.append([mouseX, mouseY, closest_peak_index, wl]) + print("added a point for peak {:.1f}nm at {}".format(wl, closest_peak_index)) + elif event == cv2.EVENT_RBUTTONDOWN: + if len(clickArray) == 0: + return + _, _, closest_peak_index, wl = clickArray[-1] + print("removing point for peak {:.1f}nm at {}".format(wl, closest_peak_index)) + clickArray = clickArray[:-1] +#listen for click on plot window +cv2.setMouseCallback(title1,handle_mouse) + + +font=cv2.FONT_HERSHEY_SIMPLEX + +raw_intensity = np.zeros(frameWidth) # array for intensity data...full of zeroes + +# messages +msg1 = "" +saveMsg = "No data saved" + +# blank image for Waterfall, filled black +waterfall = np.full([320, frameWidth, 3], fill_value=0, dtype=np.uint8) + +#Go grab the computed calibration data +caldata = readcal(frameWidth) +wavelengthData = caldata[0] + +calmsg1 = caldata[1] +calmsg2 = caldata[2] +calmsg3 = caldata[3] + +#generate the craticule data +graticuleData = generateGraticule(wavelengthData) +tens = (graticuleData[0]) +fifties = (graticuleData[1]) + +# load the banner image once +banner_image = base64.b64decode(background) +banner_image = np.frombuffer(banner_image, np.uint8) +banner_image = cv2.imdecode(banner_image, 3) + +# the background is fixed to 800x80, so we need to scale it if +# frameWidth and frameHeight have changed +banner_image_resized = np.zeros([banner_image.shape[0], frameWidth, 3], dtype=np.uint8) +w1 = banner_image.shape[1] +w2 = banner_image_resized.shape[1] +xoff = round((w2-w1)/2) +# center the background img in the newly sized background +banner_image_resized[:, xoff:xoff+w1,:] = banner_image +banner_image = banner_image_resized +# white border at the bottom +cv2.line(banner_image,(0,79),(frameWidth,79),(255,255,255),1) +messages = np.copy(banner_image) + +spectrum_vertical = None +waterfall_vertical = None +textoffset = 12 + +def build_graph_base(): + global textoffset + # blank image for Graph, filled white + result = np.full([320, frameWidth, 3], fill_value=255, dtype=np.uint8) + + # vertial lines every whole 10nm + for position in tens: + cv2.line(result,(position,15),(position,320),(200,200,200),1) + + # vertical lines every whole 50nm + for positiondata in fifties: + cv2.line(result,(positiondata[0],15),(positiondata[0],320),(0,0,0),1) + cv2.putText(result,"{}nm".format(positiondata[1]),(positiondata[0]-textoffset,12),font,0.4,(0,0,0),1, cv2.LINE_AA) + + # horizontal lines + for i in range (320): + if i >= 64 and i%64 == 0: # suppress the first line then draw the rest... + cv2.line(result,(0,i),(frameWidth,i),(100,100,100),1) + + return result + +graph_base = build_graph_base() +graph = np.copy(graph_base) + +def compute_wavelength_rgbs(wavelengthData): + result = [] + for wld in wavelengthData: + # derive the color from the wavelenthData array + rgb = wavelength_to_rgb(round(wld)) + result.append(rgb) + return result + +def compute_wavelength_rgbs_bg(array, graph): + result = np.full([320, frameWidth, 3], fill_value=255, dtype=np.uint8) + for index, wl in enumerate(array): + r, g, b = wavelength_data_rgbs[index] + result[:, index, :] = [b, g, r] + # overlay the graticule + result = cv2.bitwise_and(graph, result) + return result + +wavelength_data_rgbs = compute_wavelength_rgbs(wavelengthData) +wavelength_data_rgbs_bg = compute_wavelength_rgbs_bg(wavelength_data_rgbs, graph_base) + +def snapshot(savedata): + now = time.strftime("%Y%m%d--%H%M%S") + timenow = time.strftime("%H:%M:%S") + imdata1 = savedata[0] + graphdata = savedata[1] + if dispWaterfall: + imdata2 = savedata[2] + cv2.imwrite("waterfall-" + now + ".png",imdata2) + cv2.imwrite("spectrum-" + now + ".png",imdata1) + f = open("Spectrum-"+now+'.csv','w') + f.write('Wavelength,Intensity\r\n') + for x in zip(graphdata[0],graphdata[1]): + f.write("{},{}\r\n".format(*x)) + f.close() + message = "Last Save: "+timenow + return(message) + +if not use_picamera: + # triggers starting the cap. doing this here to avoid memory/time this takes + # from skewing measurements below + cap.isOpened() + +def runall(): + global graticuleData, tens, fifties, msg1, saveMsg, waterfall, wavelengthData, peak_intensities + global caldata, calmsg1, calmsg2, calmsg3, savpoly, mindist, thresh, messages, banner_image + global calibrate, clickArray, cursorX, cursorY, picamGain, spectrum_vertical, waterfall_vertical + global graph, graph_base, wavelength_data_rgbs, raw_intensity, picam2, cap, wavelength_data_rgbs_bg + global textoffset + + holdpeaks = False + measure = False # show cursor measurements + recPixels = False # measure pixels and record clicks + + pts = np.zeros((frameWidth+2, 2), dtype=np.int) + pts[0] = [frameWidth, 319] + pts[1] = [0, 319] + pts[2:, 0] = np.arange(frameWidth) + mask = np.zeros(graph.shape[:2], np.uint8) + while (True if use_picamera else cap.isOpened()): + # Capture frame-by-frame + frame = None + ret = True + if use_picamera: + frame = picam2.capture_array() + else: + ret, frame = cap.read() + + if not ret: + break + + x = 0 # origin of the horiz + h = 80 # height of the crop + y = int((frameHeight/2)-(h/2)) # origin of the vertical crop + w = frameWidth #width of the crop + cropped = frame[y:y+h, x:x+w] + bwimage = cv2.cvtColor(cropped,cv2.COLOR_BGR2GRAY) + rows,cols = bwimage.shape + halfway =int(rows/2) + #show our line on the original image + #now a 3px wide region + cv2.line(cropped,(0,halfway-2),(frameWidth,halfway-2),(255,255,255),1) + cv2.line(cropped,(0,halfway+2),(frameWidth,halfway+2),(255,255,255),1) + + np.copyto(messages, banner_image) + # reset the graph to the base + np.copyto(graph, graph_base) + + num_mean = 3 # average the data from this many rows of pixels + + # get the mean value for each column spanning "num_mean" rows + current_intensities = np.uint8(np.mean(bwimage[halfway-(num_mean//2):halfway+(num_mean//2)+1, :], axis=0)) + + if holdpeaks: + # find the maximums, doing so in-place + np.maximum(raw_intensity, current_intensities, casting="no", out=raw_intensity) + else: + raw_intensity = current_intensities + + np.clip(raw_intensity, 0, 255, out=raw_intensity) + + if dispWaterfall: + #data is smoothed at this point!!!!!! + #create an empty array for the data + wdata = np.zeros([1,frameWidth,3],dtype=np.uint8) + for index, i in enumerate(raw_intensity): + luminosity = i/255.0 + rgb = wavelength_data_rgbs[index] + b = int(round(rgb[0]*luminosity)) + g = int(round(rgb[1]*luminosity)) + r = int(round(rgb[2]*luminosity)) + #wdata[0,index]=(r,g,b) #fix me!!! how do we deal with this data?? + wdata[0,index] = (r,g,b) + + if use_picamera: + contrast = 2.5 + brightness =10 + wdata = cv2.addWeighted( wdata, contrast, wdata, 0, brightness) + + # rolling stream of data + for index in range(waterfall.shape[0]-1, 0, -1): + waterfall[index] = waterfall[index-1] + waterfall[0] = wdata + + intensity = None + if not holdpeaks: + intensity = savitzky_golay(raw_intensity,17,savpoly) + intensity = np.array(intensity, dtype=np.int32) + holdmsg = "Holdpeaks OFF" + else: + intensity = np.int32(raw_intensity) + holdmsg = "Holdpeaks ON" + + np.clip(intensity, 0, 255, out=intensity) + + # draw instensity data + pts[2:, 1] = 319-intensity + # mask holds where we should draw the rainbow + mask.fill(0) + cv2.drawContours(mask, [pts], -1, 1, -1, cv2.LINE_AA) + mask3d = np.dstack([mask]*3) + # draw the rainbow using the mask onto graph + np.putmask(graph, mask3d, wavelength_data_rgbs_bg) + # draw the black line + graph[pts[2:, 1], pts[2:, 0]] = 0 + + # find peaks and label them + thresh = int(thresh) #make sure the data is int. + peak_intensities = peakIndexes(intensity, thres=thresh/max(intensity), min_dist=mindist) + + for i in peak_intensities: + height = intensity[i] + height = 310-height + wavelength = round(wavelengthData[i],1) + cv2.rectangle(graph,((i-textoffset)-2,height),((i-textoffset)+60,height-15),(0,255,255),-1) + cv2.rectangle(graph,((i-textoffset)-2,height),((i-textoffset)+60,height-15),(0,0,0),1) + cv2.putText(graph,"{}nm".format(wavelength),(i-textoffset,height-3),font,0.4,(0,0,0),1, cv2.LINE_AA) + # flagpoles + cv2.line(graph,(i,height),(i,height+10),(0,0,0),1) + + if measure: + # show the cursor! + cv2.line(graph,(cursorX,cursorY-140),(cursorX,cursorY-180),(0,0,0),1) + cv2.line(graph,(cursorX-20,cursorY-160),(cursorX+20,cursorY-160),(0,0,0),1) + cv2.putText(graph,"{}nm".format(round(wavelengthData[cursorX],2)),(cursorX+5,cursorY-165),font,0.4,(0,0,0),1, cv2.LINE_AA) + + if recPixels: + # display the points + cv2.line(graph,(cursorX,cursorY-140),(cursorX,cursorY-180),(0,0,0),1) + cv2.line(graph,(cursorX-20,cursorY-160),(cursorX+20,cursorY-160),(0,0,0),1) + cv2.putText(graph,"{}px".format(cursorX),(cursorX+5,cursorY-165),font,0.4,(0,0,0),1, cv2.LINE_AA) + + for mouseX, mouseY, _, _ in clickArray: + cv2.circle(graph,(mouseX,mouseY),5,(0,0,0),-1) + #we can display text :-) so we can work out wavelength from x-pos and display it ultimately + cv2.putText(graph,str(mouseX),(mouseX+5,mouseY),cv2.FONT_HERSHEY_SIMPLEX,0.4,(0,0,0)) + + # stack the images and display the spectrum (using concatenate instead of + # vstack to reuse the array and save memory allocations/time) + if spectrum_vertical is None: + spectrum_vertical = np.concatenate((messages, cropped, graph), axis=0) + else: + np.concatenate((messages, cropped, graph), out=spectrum_vertical, axis=0) + + #print the messages + cv2.putText(spectrum_vertical,calmsg1,(message_loc1,15),font,0.4,(0,255,255),1, cv2.LINE_AA) + cv2.putText(spectrum_vertical,calmsg3,(message_loc1,33),font,0.4,(0,255,255),1, cv2.LINE_AA) + + if use_picamera: + cv2.putText(spectrum_vertical,saveMsg,(message_loc1,51),font,0.4,(0,255,255),1, cv2.LINE_AA) + cv2.putText(spectrum_vertical,"Gain: {}".format(picamGain),(message_loc1,69),font,0.4,(0,255,255),1, cv2.LINE_AA) + else: + cv2.putText(spectrum_vertical,"Framerate: {}".format(cfps),(message_loc1,51),font,0.4,(0,255,255),1, cv2.LINE_AA) + cv2.putText(spectrum_vertical,saveMsg,(message_loc1,69),font,0.4,(0,255,255),1, cv2.LINE_AA) + + #Second column + cv2.putText(spectrum_vertical,holdmsg,(message_loc2,15),font,0.4,(0,255,255),1, cv2.LINE_AA) + cv2.putText(spectrum_vertical,"Savgol Filter: {}".format(savpoly),(message_loc2,33),font,0.4,(0,255,255),1, cv2.LINE_AA) + cv2.putText(spectrum_vertical,"Label Peak Width: {}".format(mindist),(message_loc2,51),font,0.4,(0,255,255),1, cv2.LINE_AA) + cv2.putText(spectrum_vertical,"Label Threshold: {}".format(thresh),(message_loc2,69),font,0.4,(0,255,255),1, cv2.LINE_AA) + cv2.imshow(title1,spectrum_vertical) + + if dispWaterfall: + #stack the images and display the waterfall + if waterfall_vertical is None: + waterfall_vertical = np.concatenate((messages, cropped, waterfall), axis=0) + else: + np.concatenate((messages, cropped, waterfall), out=waterfall_vertical, axis=0) + + #dividing lines... + cv2.line(waterfall_vertical,(0,80),(frameWidth,80),(255,255,255),1) + cv2.line(waterfall_vertical,(0,160),(frameWidth,160),(255,255,255),1) + #Draw this stuff over the top of the image! + #Display a graticule calibrated with cal data + + #vertical lines every whole 50nm + for positiondata in fifties: + for i in range(162,480): + if i%20 == 0: + cv2.line(waterfall_vertical,(positiondata[0],i),(positiondata[0],i+1),(0,0,0),2) + cv2.line(waterfall_vertical,(positiondata[0],i),(positiondata[0],i+1),(255,255,255),1) + cv2.putText(waterfall_vertical,"{}nm".format(positiondata[1]),(positiondata[0]-textoffset,475),font,0.4,(0,0,0),2, cv2.LINE_AA) + cv2.putText(waterfall_vertical,"{}nm".format(positiondata[1]),(positiondata[0]-textoffset,475),font,0.4,(255,255,255),1, cv2.LINE_AA) + + cv2.putText(waterfall_vertical,calmsg1,(message_loc1,15),font,0.4,(0,255,255),1, cv2.LINE_AA) + + if use_picamera: + cv2.putText(waterfall_vertical,calmsg3,(message_loc1,33),font,0.4,(0,255,255),1, cv2.LINE_AA) + cv2.putText(waterfall_vertical,saveMsg,(message_loc1,51),font,0.4,(0,255,255),1, cv2.LINE_AA) + cv2.putText(waterfall_vertical,"Gain: {}".format(picamGain),(message_loc1,69),font,0.4,(0,255,255),1, cv2.LINE_AA) + else: + cv2.putText(waterfall_vertical,calmsg2,(message_loc1,33),font,0.4,(0,255,255),1, cv2.LINE_AA) + cv2.putText(waterfall_vertical,calmsg3,(message_loc1,51),font,0.4,(0,255,255),1, cv2.LINE_AA) + cv2.putText(waterfall_vertical,saveMsg,(message_loc1,69),font,0.4,(0,255,255),1, cv2.LINE_AA) + + cv2.putText(waterfall_vertical,holdmsg,(message_loc2,15),font,0.4,(0,255,255),1, cv2.LINE_AA) + + cv2.imshow(title2,waterfall_vertical) + + while True: + keyPress = cv2.waitKey(1) + if keyPress == -1: + break + elif keyPress == ord("q"): + return + elif keyPress == ord("h"): + holdpeaks = not holdpeaks + elif keyPress == ord("s"): + #package up the data! + graphdata = [] + graphdata.append(wavelengthData) + graphdata.append(intensity) + if dispWaterfall: + savedata = [] + savedata.append(spectrum_vertical) + savedata.append(graphdata) + savedata.append(waterfall_vertical) + else: + savedata = [] + savedata.append(spectrum_vertical) + savedata.append(graphdata) + saveMsg = snapshot(savedata) + elif keyPress == ord("c"): + if writecal(clickArray, frameWidth): + #overwrite wavelength data + #Go grab the computed calibration data + caldata = readcal(frameWidth) + calmsg1 = caldata[1] + calmsg2 = caldata[2] + calmsg3 = caldata[3] + #overwrite graticule data + graticuleData = generateGraticule(wavelengthData) + tens = (graticuleData[0]) + fifties = (graticuleData[1]) + graph_base = build_graph_base() + wavelengthData = caldata[0] + wavelength_data_rgbs = compute_wavelength_rgbs(wavelengthData) + wavelength_data_rgbs_bg = compute_wavelength_rgbs_bg(wavelength_data_rgbs, graph_base) + elif keyPress == ord("x"): + clickArray = [] + elif keyPress == ord("m"): + measure = not measure + recPixels = False # turn off recpixels! + clickArray = [] + elif keyPress == ord("p"): + measure = False #turn off measure! + recPixels = not recPixels + if not recPixels: + clickArray = [] + elif keyPress == ord("o"):#sav up + savpoly+=1 + if savpoly >=15: + savpoly=15 + elif keyPress == ord("l"):#sav down + savpoly-=1 + if savpoly <=0: + savpoly=0 + elif keyPress == ord("i"):#Peak width up + mindist+=1 + if mindist >=100: + mindist=100 + elif keyPress == ord("k"):#Peak Width down + mindist-=1 + if mindist <=0: + mindist=0 + elif keyPress == ord("u"):#label thresh up + thresh+=1 + if thresh >=100: + thresh=100 + elif keyPress == ord("j"):#label thresh down + thresh-=1 + if thresh <=0: + thresh=0 + elif use_picamera and keyPress == ord("t"):#Gain up! + picamGain += 1 + if picamGain >=50: + picamGain = 50.0 + picam2.set_controls({"AnalogueGain": picamGain}) + print("Camera Gain: {}".format(picamGain)) + elif use_picamera and keyPress == ord("g"):#Gain down + picamGain -= 1 + if picamGain <=0: + picamGain = 0.0 + picam2.set_controls({"AnalogueGain": picamGain}) + print("Camera Gain: {}".format(picamGain)) + +runall() + +#Everything done +if use_device or use_webcam: + cap.release() +cv2.destroyAllWindows() \ No newline at end of file diff --git a/src/specFunctions.py b/src/specFunctions.py index 72c4cee..5b1f5b6 100644 --- a/src/specFunctions.py +++ b/src/specFunctions.py @@ -30,413 +30,419 @@ import time def wavelength_to_rgb(nm): - #from: Chris Webb https://www.codedrome.com/exploring-the-visible-spectrum-in-python/ - #returns RGB vals for a given wavelength - gamma = 0.8 - max_intensity = 255 - factor = 0 - rgb = {"R": 0, "G": 0, "B": 0} - if 380 <= nm <= 439: - rgb["R"] = -(nm - 440) / (440 - 380) - rgb["G"] = 0.0 - rgb["B"] = 1.0 - elif 440 <= nm <= 489: - rgb["R"] = 0.0 - rgb["G"] = (nm - 440) / (490 - 440) - rgb["B"] = 1.0 - elif 490 <= nm <= 509: - rgb["R"] = 0.0 - rgb["G"] = 1.0 - rgb["B"] = -(nm - 510) / (510 - 490) - elif 510 <= nm <= 579: - rgb["R"] = (nm - 510) / (580 - 510) - rgb["G"] = 1.0 - rgb["B"] = 0.0 - elif 580 <= nm <= 644: - rgb["R"] = 1.0 - rgb["G"] = -(nm - 645) / (645 - 580) - rgb["B"] = 0.0 - elif 645 <= nm <= 780: - rgb["R"] = 1.0 - rgb["G"] = 0.0 - rgb["B"] = 0.0 - if 380 <= nm <= 419: - factor = 0.3 + 0.7 * (nm - 380) / (420 - 380) - elif 420 <= nm <= 700: - factor = 1.0 - elif 701 <= nm <= 780: - factor = 0.3 + 0.7 * (780 - nm) / (780 - 700) - if rgb["R"] > 0: - rgb["R"] = int(max_intensity * ((rgb["R"] * factor) ** gamma)) - else: - rgb["R"] = 0 - if rgb["G"] > 0: - rgb["G"] = int(max_intensity * ((rgb["G"] * factor) ** gamma)) - else: - rgb["G"] = 0 - if rgb["B"] > 0: - rgb["B"] = int(max_intensity * ((rgb["B"] * factor) ** gamma)) - else: - rgb["B"] = 0 - #display no color as gray - if(rgb["R"]+rgb["G"]+rgb["B"]) == 0: - rgb["R"] = 155 - rgb["G"] = 155 - rgb["B"] = 155 - return (rgb["R"], rgb["G"], rgb["B"]) + #from: Chris Webb https://www.codedrome.com/exploring-the-visible-spectrum-in-python/ + #returns RGB vals for a given wavelength + gamma = 0.8 + max_intensity = 255 + factor = 0 + rgb = {"R": 0, "G": 0, "B": 0} + if 380 <= nm <= 439: + rgb["R"] = -(nm - 440) / (440 - 380) + rgb["G"] = 0.0 + rgb["B"] = 1.0 + elif 440 <= nm <= 489: + rgb["R"] = 0.0 + rgb["G"] = (nm - 440) / (490 - 440) + rgb["B"] = 1.0 + elif 490 <= nm <= 509: + rgb["R"] = 0.0 + rgb["G"] = 1.0 + rgb["B"] = -(nm - 510) / (510 - 490) + elif 510 <= nm <= 579: + rgb["R"] = (nm - 510) / (580 - 510) + rgb["G"] = 1.0 + rgb["B"] = 0.0 + elif 580 <= nm <= 644: + rgb["R"] = 1.0 + rgb["G"] = -(nm - 645) / (645 - 580) + rgb["B"] = 0.0 + elif 645 <= nm <= 780: + rgb["R"] = 1.0 + rgb["G"] = 0.0 + rgb["B"] = 0.0 + if 380 <= nm <= 419: + factor = 0.3 + 0.7 * (nm - 380) / (420 - 380) + elif 420 <= nm <= 700: + factor = 1.0 + elif 701 <= nm <= 780: + factor = 0.3 + 0.7 * (780 - nm) / (780 - 700) + if rgb["R"] > 0: + rgb["R"] = int(max_intensity * ((rgb["R"] * factor) ** gamma)) + else: + rgb["R"] = 0 + if rgb["G"] > 0: + rgb["G"] = int(max_intensity * ((rgb["G"] * factor) ** gamma)) + else: + rgb["G"] = 0 + if rgb["B"] > 0: + rgb["B"] = int(max_intensity * ((rgb["B"] * factor) ** gamma)) + else: + rgb["B"] = 0 + #display no color as gray + if(rgb["R"]+rgb["G"]+rgb["B"]) == 0: + rgb["R"] = 155 + rgb["G"] = 155 + rgb["B"] = 155 + return (rgb["R"], rgb["G"], rgb["B"]) def savitzky_golay(y, window_size, order, deriv=0, rate=1): - #scipy - #From: https://scipy.github.io/old-wiki/pages/Cookbook/SavitzkyGolay - ''' - Copyright (c) 2001-2002 Enthought, Inc. 2003-2022, SciPy Developers. - All rights reserved. - - Redistribution and use in source and binary forms, with or without - modification, are permitted provided that the following conditions - are met: - - 1. Redistributions of source code must retain the above copyright - notice, this list of conditions and the following disclaimer. - - 2. Redistributions in binary form must reproduce the above - copyright notice, this list of conditions and the following - disclaimer in the documentation and/or other materials provided - with the distribution. - - 3. Neither the name of the copyright holder nor the names of its - contributors may be used to endorse or promote products derived - from this software without specific prior written permission. - - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS - "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT - LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR - A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT - OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, - SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT - LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, - DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY - THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT - (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE - OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - ''' - import numpy as np - from math import factorial - try: - window_size = np.abs(np.int(window_size)) - order = np.abs(np.int(order)) - except ValueError: - raise ValueError("window_size and order have to be of type int") - if window_size % 2 != 1 or window_size < 1: - raise TypeError("window_size size must be a positive odd number") - if window_size < order + 2: - raise TypeError("window_size is too small for the polynomials order") - order_range = range(order+1) - half_window = (window_size -1) // 2 - # precompute coefficients - b = np.mat([[k**i for i in order_range] for k in range(-half_window, half_window+1)]) - m = np.linalg.pinv(b).A[deriv] * rate**deriv * factorial(deriv) - # pad the signal at the extremes with - # values taken from the signal itself - firstvals = y[0] - np.abs( y[1:half_window+1][::-1] - y[0] ) - lastvals = y[-1] + np.abs(y[-half_window-1:-1][::-1] - y[-1]) - y = np.concatenate((firstvals, y, lastvals)) - return np.convolve( m[::-1], y, mode='valid') + #scipy + #From: https://scipy.github.io/old-wiki/pages/Cookbook/SavitzkyGolay + ''' + Copyright (c) 2001-2002 Enthought, Inc. 2003-2022, SciPy Developers. + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + 3. Neither the name of the copyright holder nor the names of its + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + ''' + import numpy as np + from math import factorial + try: + window_size = np.abs(np.int(window_size)) + order = np.abs(np.int(order)) + except ValueError: + raise ValueError("window_size and order have to be of type int") + if window_size % 2 != 1 or window_size < 1: + raise TypeError("window_size size must be a positive odd number") + if window_size < order + 2: + raise TypeError("window_size is too small for the polynomials order") + order_range = range(order+1) + half_window = (window_size -1) // 2 + # precompute coefficients + b = np.mat([[k**i for i in order_range] for k in range(-half_window, half_window+1)]) + m = np.linalg.pinv(b).A[deriv] * rate**deriv * factorial(deriv) + # pad the signal at the extremes with + # values taken from the signal itself + firstvals = y[0] - np.abs( y[1:half_window+1][::-1] - y[0] ) + lastvals = y[-1] + np.abs(y[-half_window-1:-1][::-1] - y[-1]) + y = np.concatenate((firstvals, y, lastvals)) + return np.convolve( m[::-1], y, mode='valid') def peakIndexes(y, thres=0.3, min_dist=1, thres_abs=False): - #from peakutils - #from https://bitbucket.org/lucashnegri/peakutils/raw/f48d65a9b55f61fb65f368b75a2c53cbce132a0c/peakutils/peak.py - ''' - The MIT License (MIT) - - Copyright (c) 2014-2022 Lucas Hermann Negri - - Permission is hereby granted, free of charge, to any person obtaining a copy - of this software and associated documentation files (the "Software"), to deal - in the Software without restriction, including without limitation the rights - to use, copy, modify, merge, publish, distribute, sublicense, and/or sell - copies of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be included in - all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR - IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE - AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER - LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, - OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN - THE SOFTWARE. - ''' - if isinstance(y, np.ndarray) and np.issubdtype(y.dtype, np.unsignedinteger): - raise ValueError("y must be signed") - - if not thres_abs: - thres = thres * (np.max(y) - np.min(y)) + np.min(y) - - min_dist = int(min_dist) - - # compute first order difference - dy = np.diff(y) - - # propagate left and right values successively to fill all plateau pixels (0-value) - zeros, = np.where(dy == 0) - - # check if the signal is totally flat - if len(zeros) == len(y) - 1: - return np.array([]) - - if len(zeros): - # compute first order difference of zero indexes - zeros_diff = np.diff(zeros) - # check when zeros are not chained together - zeros_diff_not_one, = np.add(np.where(zeros_diff != 1), 1) - # make an array of the chained zero indexes - zero_plateaus = np.split(zeros, zeros_diff_not_one) - - # fix if leftmost value in dy is zero - if zero_plateaus[0][0] == 0: - dy[zero_plateaus[0]] = dy[zero_plateaus[0][-1] + 1] - zero_plateaus.pop(0) - - # fix if rightmost value of dy is zero - if len(zero_plateaus) and zero_plateaus[-1][-1] == len(dy) - 1: - dy[zero_plateaus[-1]] = dy[zero_plateaus[-1][0] - 1] - zero_plateaus.pop(-1) - - # for each chain of zero indexes - for plateau in zero_plateaus: - median = np.median(plateau) - # set leftmost values to leftmost non zero values - dy[plateau[plateau < median]] = dy[plateau[0] - 1] - # set rightmost and middle values to rightmost non zero values - dy[plateau[plateau >= median]] = dy[plateau[-1] + 1] - - # find the peaks by using the first order difference - peaks = np.where( - (np.hstack([dy, 0.0]) < 0.0) - & (np.hstack([0.0, dy]) > 0.0) - & (np.greater(y, thres)) - )[0] - - # handle multiple peaks, respecting the minimum distance - if peaks.size > 1 and min_dist > 1: - highest = peaks[np.argsort(y[peaks])][::-1] - rem = np.ones(y.size, dtype=bool) - rem[peaks] = False - - for peak in highest: - if not rem[peak]: - sl = slice(max(0, peak - min_dist), peak + min_dist + 1) - rem[sl] = True - rem[peak] = False - - peaks = np.arange(y.size)[~rem] - - return peaks + #from peakutils + #from https://bitbucket.org/lucashnegri/peakutils/raw/f48d65a9b55f61fb65f368b75a2c53cbce132a0c/peakutils/peak.py + ''' + The MIT License (MIT) + + Copyright (c) 2014-2022 Lucas Hermann Negri + + Permission is hereby granted, free of charge, to any person obtaining a copy + of this software and associated documentation files (the "Software"), to deal + in the Software without restriction, including without limitation the rights + to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + copies of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be included in + all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN + THE SOFTWARE. + ''' + if isinstance(y, np.ndarray) and np.issubdtype(y.dtype, np.unsignedinteger): + raise ValueError("y must be signed") + + if not thres_abs: + thres = thres * (np.max(y) - np.min(y)) + np.min(y) + + min_dist = int(min_dist) + + # compute first order difference + dy = np.diff(y) + + # propagate left and right values successively to fill all plateau pixels (0-value) + zeros, = np.where(dy == 0) + + # check if the signal is totally flat + if len(zeros) == len(y) - 1: + return np.array([]) + + if len(zeros): + # compute first order difference of zero indexes + zeros_diff = np.diff(zeros) + # check when zeros are not chained together + zeros_diff_not_one, = np.add(np.where(zeros_diff != 1), 1) + # make an array of the chained zero indexes + zero_plateaus = np.split(zeros, zeros_diff_not_one) + + # fix if leftmost value in dy is zero + if zero_plateaus[0][0] == 0: + dy[zero_plateaus[0]] = dy[zero_plateaus[0][-1] + 1] + zero_plateaus.pop(0) + + # fix if rightmost value of dy is zero + if len(zero_plateaus) and zero_plateaus[-1][-1] == len(dy) - 1: + dy[zero_plateaus[-1]] = dy[zero_plateaus[-1][0] - 1] + zero_plateaus.pop(-1) + + # for each chain of zero indexes + for plateau in zero_plateaus: + lplat = plateau.shape[0] + left = dy[plateau[0] - 1] + right = dy[plateau[lplat-1] + 1] + middle = lplat // 2 + # always in sorted order + median = plateau[middle] + for p in plateau: + if p < median: + dy[p] = left + else: + dy[p] = right + + # find the peaks by using the first order difference + peaks = np.where( + (np.hstack([dy, 0.0]) < 0.0) + & (np.hstack([0.0, dy]) > 0.0) + & (np.greater(y, thres)) + )[0] + + # handle multiple peaks, respecting the minimum distance + if peaks.size > 1 and min_dist > 1: + highest = peaks[np.argsort(y[peaks])][::-1] + rem = np.ones(y.size, dtype=bool) + rem[peaks] = False + + for peak in highest: + if not rem[peak]: + sl = slice(max(0, peak - min_dist), peak + min_dist + 1) + rem[sl] = True + rem[peak] = False + + peaks = np.arange(y.size)[~rem] + + return peaks + +def compute_wavelengths(width, pixels, wavelengths, errors): + wavelengthData = [] + calibration_error = 1000.0 + message = None + + if len(pixels) == 3: + print("Calculating second order polynomial...") + coefficients = np.poly1d(np.polyfit(pixels, wavelengths, 2)) + print(coefficients) + C1 = coefficients[2] + C2 = coefficients[1] + C3 = coefficients[0] + print("Generating Wavelength Data!\n\n") + for pixel in range(width): + wavelength=((C1*pixel**2)+(C2*pixel)+C3) + wavelength = round(wavelength,6) #because seriously! + wavelengthData.append(wavelength) + print("Done! Note that calibration with only 3 wavelengths will not be accurate!") + if errors == 1: + message = 0 #return message zero(errors) + else: + message = 1 #return message only 3 wavelength cal secodn order poly (Inaccurate) + + if len(pixels) > 3: + print("Calculating third order polynomial...") + coefficients = np.poly1d(np.polyfit(pixels, wavelengths, 3)) + print(coefficients) + #note this pulls out extremely precise numbers. + #this causes slight differences in vals then when we compute manual, but hey ho, more precision + #that said, we waste that precision later, but tbh, we wouldn't get that kind of precision in + #the real world anyway! 1/10 of a nm is more than adequate! + C1 = coefficients[3] + C2 = coefficients[2] + C3 = coefficients[1] + C4 = coefficients[0] + ''' + print(C1) + print(C2) + print(C3) + print(C4) + ''' + print("Generating Wavelength Data!\n\n") + for pixel in range(width): + wavelength=((C1*pixel**3)+(C2*pixel**2)+(C3*pixel)+C4) + wavelength = round(wavelength,6) + wavelengthData.append(wavelength) + + #final job, we need to compare all the recorded wavelenths with predicted wavelengths + #and note the deviation! + #do something if it is too big! + predicted = [] + #iterate over the original pixelnumber array and predict results + for i in pixels: + px = i + y=((C1*px**3)+(C2*px**2)+(C3*px)+C4) + predicted.append(y) + + #calculate 2 squared of the result + #if this is close to 1 we are all good! + corr_matrix = np.corrcoef(wavelengths, predicted) + corr = corr_matrix[0,1] + R_sq = corr**2 + calibration_error = R_sq + + print("R-Squared={}".format(R_sq)) + + message = 2 #Multiwavelength cal, 3rd order poly + + return wavelengthData, message, calibration_error def readcal(width): - #read in the calibration points - #compute second or third order polynimial, and generate wavelength array! - #Les Wright 28 Sept 2022 - errors = 0 - message = 0 #variable to store returned message data - try: - print("Loading calibration data...") - file = open('caldata.txt', 'r') - except: - errors = 1 - - try: - #read both the pixel numbers and wavelengths into two arrays. - lines = file.readlines() - line0 = lines[0].strip() #strip newline - pixels = line0.split(',') #split on , - pixels = [int(i) for i in pixels] #convert list of strings to ints - line1 = lines[1].strip() - wavelengths = line1.split(',') - wavelengths = [float(i) for i in wavelengths]#convert list of strings to floats - except: - errors = 1 - - try: - if (len(pixels) != len(wavelengths)): - #The Calibration points are of unequal length! - errors = 1 - if (len(pixels) < 3): - #The Cal data contains less than 3 pixels! - errors = 1 - if (len(wavelengths) < 3): - #The Cal data contains less than 3 wavelengths! - errors = 1 - except: - errors = 1 - - if errors == 1: - print("Loading of Calibration data failed (missing caldata.txt or corrupted data!") - print("Loading placeholder data...") - print("You MUST perform a Calibration to use this software!\n\n") - pixels = [0,400,800] - wavelengths = [380,560,750] - - - #create an array for the data... - wavelengthData = [] - - if (len(pixels) == 3): - print("Calculating second order polynomial...") - coefficients = np.poly1d(np.polyfit(pixels, wavelengths, 2)) - print(coefficients) - C1 = coefficients[2] - C2 = coefficients[1] - C3 = coefficients[0] - print("Generating Wavelength Data!\n\n") - for pixel in range(width): - wavelength=((C1*pixel**2)+(C2*pixel)+C3) - wavelength = round(wavelength,6) #because seriously! - wavelengthData.append(wavelength) - print("Done! Note that calibration with only 3 wavelengths will not be accurate!") - if errors == 1: - message = 0 #return message zero(errors) - else: - message = 1 #return message only 3 wavelength cal secodn order poly (Inaccurate) - - if (len(pixels) > 3): - print("Calculating third order polynomial...") - coefficients = np.poly1d(np.polyfit(pixels, wavelengths, 3)) - print(coefficients) - #note this pulls out extremely precise numbers. - #this causes slight differences in vals then when we compute manual, but hey ho, more precision - #that said, we waste that precision later, but tbh, we wouldn't get that kind of precision in - #the real world anyway! 1/10 of a nm is more than adequate! - C1 = coefficients[3] - C2 = coefficients[2] - C3 = coefficients[1] - C4 = coefficients[0] - ''' - print(C1) - print(C2) - print(C3) - print(C4) - ''' - print("Generating Wavelength Data!\n\n") - for pixel in range(width): - wavelength=((C1*pixel**3)+(C2*pixel**2)+(C3*pixel)+C4) - wavelength = round(wavelength,6) - wavelengthData.append(wavelength) - - #final job, we need to compare all the recorded wavelenths with predicted wavelengths - #and note the deviation! - #do something if it is too big! - predicted = [] - #iterate over the original pixelnumber array and predict results - for i in pixels: - px = i - y=((C1*px**3)+(C2*px**2)+(C3*px)+C4) - predicted.append(y) - - #calculate 2 squared of the result - #if this is close to 1 we are all good! - corr_matrix = np.corrcoef(wavelengths, predicted) - corr = corr_matrix[0,1] - R_sq = corr**2 - - print("R-Squared="+str(R_sq)) - - message = 2 #Multiwavelength cal, 3rd order poly - - - if message == 0: - calmsg1 = "UNCALIBRATED!" - calmsg2 = "Defaults loaded" - calmsg3 = "Perform Calibration!" - if message == 1: - calmsg1 = "Calibrated!!" - calmsg2 = "Using 3 cal points" - calmsg3 = "2nd Order Polyfit" - if message == 2: - calmsg1 = "Calibrated!!!" - calmsg2 = "Using > 3 cal points" - calmsg3 = "3rd Order Polyfit" - - returndata = [] - returndata.append(wavelengthData) - returndata.append(calmsg1) - returndata.append(calmsg2) - returndata.append(calmsg3) - return returndata - - -def writecal(clickArray): - calcomplete = False - pxdata = [] - wldata = [] - print("Enter known wavelengths for observed pixels!") - for i in clickArray: - pixel = i[0] - wavelength = input("Enter wavelength for: "+str(pixel)+"px:") - pxdata.append(pixel) - wldata.append(wavelength) - #This try except serves two purposes - #first I want to write data to the caldata.txt file without quotes - #second it validates the data in as far as no strings were entered - try: - wldata = [float(x) for x in wldata] - except: - print("Only ints or decimals are allowed!") - print("Calibration aborted!") - - pxdata = ','.join(map(str, pxdata)) #convert array to string - wldata = ','.join(map(str, wldata)) #convert array to string - f = open('caldata.txt','w') - f.write(pxdata+'\r\n') - f.write(wldata+'\r\n') - print("Calibration Data Written!") - calcomplete = True - return calcomplete + #read in the calibration points + #compute second or third order polynimial, and generate wavelength array! + #Les Wright 28 Sept 2022 + errors = 0 + message = 0 #variable to store returned message data + calibrated_frame_width = None + try: + print("Loading calibration data...") + with open('caldata.txt', 'r') as f: + #read both the pixel numbers and wavelengths into two arrays. + lines = f.readlines() + calibrated_frame_width = lines[0].strip() + calibrated_frame_width = float(calibrated_frame_width) + ratio = width / calibrated_frame_width + line0 = lines[1].strip() #strip newline + pixels = line0.split(',') #split on , + pixels = [int(int(i) * ratio) for i in pixels] #convert list of strings to ints + line1 = lines[2].strip() + wavelengths = line1.split(',') + wavelengths = [float(i) for i in wavelengths]#convert list of strings to floats + except: + errors = 1 + + try: + if (len(pixels) != len(wavelengths)): + #The Calibration points are of unequal length! + errors = 1 + if (len(pixels) < 3): + #The Cal data contains less than 3 pixels! + errors = 1 + if (len(wavelengths) < 3): + #The Cal data contains less than 3 wavelengths! + errors = 1 + except: + errors = 1 + + if errors == 1: + print("Loading of Calibration data failed (missing caldata.txt or corrupted data!") + print("Loading placeholder data...") + print("You MUST perform a Calibration to use this software!\n\n") + pixels = [0,width/2,width] + wavelengths = [380,560,750] + + wavelengthData, _, _ = compute_wavelengths(width, pixels, wavelengths, errors) + + if message == 0: + calmsg1 = "UNCALIBRATED!" + calmsg2 = "Defaults loaded" + calmsg3 = "Perform Calibration!" + if message == 1: + calmsg1 = "Calibrated!!" + calmsg2 = "Using 3 cal points" + calmsg3 = "2nd Order Polyfit" + if message == 2: + calmsg1 = "Calibrated!!!" + calmsg2 = "Using > 3 cal points" + calmsg3 = "3rd Order Polyfit" + + return [wavelengthData, calmsg1, calmsg2, calmsg3] + + +def closest(lst, K): + # https://www.geeksforgeeks.org/python-find-closest-number-to-k-in-given-list/ + return lst[min(range(len(lst)), key = lambda i: abs(lst[i]-K))] + + + +def writecal(clickArray, frameWidth): + pxdata = [] + wldata = [] + + print("Enter known wavelengths for observed pixels!") + for _, _, pixel, wavelength_displayed in sorted(clickArray, key=lambda x: x[0]): + while True: + wavelength = input("Enter wavelength for {}px ({:.1f}nm): ".format(pixel, wavelength_displayed)) + try: + wavelength = float(wavelength) + except: + print("Only ints or decimals are allowed!") + continue + + pxdata.append(pixel) + wldata.append(wavelength) + break + + _, _, calibration_error = compute_wavelengths(frameWidth, pxdata, wldata, 0) + while True: + save = input("R-squared={}. Overwrite current calibration data and save new values? ".format(calibration_error)) + if save == "n": + return False + elif save != "y": + print("Invalid value '{}' entered. Please enter 'y' or 'n'.".format(save)) + continue + else: + break + + pxdata = ','.join(map(str, pxdata)) # convert array to string + wldata = ','.join(map(str, wldata)) # convert array to string + with open('caldata.txt','w') as f: + f.write(str(frameWidth)+'\r\n') + f.write(pxdata+'\r\n') + f.write(wldata+'\r\n') + + print("Calibration Data Written!") + return True def generateGraticule(wavelengthData): - low = wavelengthData[0] #get lowet number in list - high = wavelengthData[len(wavelengthData)-1] #get highest number - #round and int these numbers so we have our range of numbers to look at - #give a margin of 10 at each end for good measure - low = int(round(low))-10 - high = int(round(high))+10 - #print('...') - #print(low) - #print(high) - #print('...') - returndata = [] - #find positions of every whole 10nm - tens = [] - for i in range(low,high): - if (i%10==0): - #position contains pixelnumber and wavelength - position = min(enumerate(wavelengthData), key=lambda x: abs(i - x[1])) - #If the difference between the target and result is <9 show the line - #(otherwise depending on the scale we get dozens of number either end that are close to the target) - if abs(i-position[1]) <1: - #print(i) - #print(position) - tens.append(position[0]) - returndata.append(tens) - fifties = [] - for i in range(low,high): - if (i%50==0): - #position contains pixelnumber and wavelength - position = min(enumerate(wavelengthData), key=lambda x: abs(i - x[1])) - #If the difference between the target and result is <1 show the line - #(otherwise depending on the scale we get dozens of number either end that are close to the target) - if abs(i-position[1]) <1: - labelpos = position[0] - labeltxt = int(round(position[1])) - labeldata = [labelpos,labeltxt] - fifties.append(labeldata) - returndata.append(fifties) - return returndata + low = wavelengthData[0] + high = wavelengthData[-1] + + # find the closest 10th to the start and end + first_ten = int(round(np.ceil(low/10)*10)) # closest without going under + last_ten = int(round(np.floor(high/10)*10)) # closest without going over + + # find positions of every whole 10nm + tens = [] + fifties = [] + for tenth_wl in range(first_ten, last_ten+1, 10): + pixel, wl = min(enumerate(wavelengthData), key=lambda x: abs(tenth_wl - x[1])) + # If the difference between the target and result is <9 show the line + # (otherwise depending on the scale we get dozens of number either end that are close to the target) + if abs(tenth_wl - wl) < 1: + if tenth_wl % 50 == 0: + fifties.append([pixel, int(round(wl))]) + else: + tens.append(pixel) + + return [tens, fifties] background = 'iVBORw0KGgoAAAANSUhEUgAAAyAAAABQCAYAAADhuhE0AAAABmJLR0QA/wD/AP+gvaeTAAAACXBIWXMAAAsTAAALEwEAmpwYAAAAB3RJTUUH5goFFDgj33B8iQAAAB1pVFh0Q29tbWVudAAAAAAAQ3JlYXRlZCB3aXRoIEdJTVBkLmUHAAAgAElEQVR42u2deXxU1fn/P3e2rJM9kIQEEgKEEEIgrEkgiNSq0C+i2BZBlOKGWrTWal1a/epXbf1ZtOJSpWqxImoFXKkoEkzYZA2QhYSwJWQh62Qyk2Uyy/P7IwuZmXvvzISEBHner9d5zcydc892z73n+Zx7FgEAgWEYhmEYhmEY5hKg4CJgGIZhGIZhGIYFCMMwDMMwDMMwLEAYhmEYhmEYhmFYgDAMwzAMwzAMwwKEYRiGYRiGYRiGBQjDMAzDMAzDMCxAGIZhGIZhGIZhAcIwDMMwDMMwDMMChGEYhmEYhmEYFiAMwzAMwzAMwzAsQBiGYRiGYRiGYQHCMAzDMAzDMAwLEIZhGIZhGIZhGBYgDMMwDMMwDMOwAGEYhmEYhmEYhmEBwjAMwzAMwzAMCxCGYRiGYRiGYViAMAzDMAzDMAzDsABhGIZhGIZhGIYFCMMwDMMwDMMwDAsQhmEYhmEYhmEGEyougr5D6xOKxZOeRbRyFhR1sWjT+aDFoICpVYDgDSh8bFAHmaEMaoLFtwY6ysOZxmwcPr0ZhpZaLkCGYa5oYkMDcfuU0Zgb4YdYwYRQQz28jHoIBiOENhNI4w+rxh+GgHCcUAXg67pm/HX3YViJuPAYhmEuIwQAffrkJnwLQNmpbbo+FQ6fPb8rQWoFbFolLFpAN6QZBb4n8FHdRrxbsO6yKMRA33A8MmULrIcno82ogAWAVcRJHbcpCeEJRiA8H2uz07lWMgxzxTAvORYPTR+BdNt5+JwrhmCTeVhanB+mlrBIfDlkNH61dScLEYZhmCtXgHzXKTwUEiJEKXNMafe7Nd6Cfw5Zhwf3PjRoCzBj1I24qe0TNJWrXYoNCwCbeBva7b6BwLWSYZgrBrpB5iHp+EC1ST9cDfGJmFVSj6PVNVyoDMMwV64AUbopNhROwsPRX+W0GkzKn4yalsHVsEyIzsRyUxaMtUqXbzkc/+tqZx2PswC5xJwl5+raVQW1fC0YcTaQ2JOqw13L1cazNmMh5N94SD0sRQRL27A4TCzTo7i+gQuWYRhmENMPk9BJRNNQp9ZBj/+EHp+Cgyay9xe1PxKFYwqh1WgHVeHdE74JrbVKp9QLEqUipfzYXhmMVhEXAcPV5JIXoth3qaZDBO+yM/h24hguU4ZhmCtPgEDE1JZrTcit80KOhOK/ad8OmoK7KeVBGHPDJNtSKSEiJ7WYQYLAqpDx7OkG7kzom/tNAEihQOXIaVgTnIb0qij45qow4nQA/qCdiMqx08T7t3p8H3HoRywfP47LlWEYZhDTD6tgOVpvBPn3Au6MAus4d8beGZgaOQ0HqvYPeMFlBCyHWaT968pl0JTzKFJ8gNzKr3GuvgDNbY0I9BuKUG0MRg3NQJR2OgKsqdCfiYWuihcjGzDIzWMMI/JUIq46fXMPdhbcudjpWJ5dhaxt9s/4Ml0TVu86gtUAPluYiYWFOdKKj4D7I8OwLp+LlmEY5goSIATnYVgXfgtItvPtpdQgKSwRd8Qvw4qqxfA+oxRpyqmjnWkX8PjoJ3FT1Q0DXnB+5xOgl2j/vGbvwZ+zM5yO64yV0BkrcbJqn93xcSOuQXrso7CUzwROcaUctMKEYbiK9AvmkGg8eDYI/1i7z6XfGz/PQfkN0zAsb7+kMhxXe54LlWEY5soSIHLjV5ybbJO1HYerj+Jw9VE85fMCzo49Bv8iP8kw0gwZg6Lg2us1djnuybvHlnsUVmHpNhSWbut7keQXgdSUZ2FqX4Dq2iDU6jRQqghhQ9oRGVWO9pb1OLTvORBZex1H7Ph5iJm4HObA6ahtD4fBqkFjmwIqLyDI14JA3zaEowgtJ79G7lcvwmoxXVSeYmdci5jrl8MyagZq/YbCIKihF5QwKwBvFSFEbUEImhHQUAZb8QGc/XwdKg7uHnQ3XnjyGCT+773QX5WIBn8F6tRmKMiC4S3Aca0LgS0ImPDbGxB842Q0TQhCg9aCJlUbmoRWaElAgEWJYKMSQfkG6L/Ix9FX/wuy2volH3HXj0HM8wmoSKxFhVclbDBjiCUYw8sj0PJGI46sPix5bvozU+C/XIXKqFI0KevgS94IM4ci/HQ0zr9twN5Xj15U2lLmxWLa8hiETzfDFl4Lq8YAs6IRGqjgbQmCT1sgTEXhKPq6BZtezIXZZO3V066/xMh182Jx2/IYpE03Y2h4LTRKAxRoBCwqWNqD0NYciKKScHy9tQUvvpILU7u1z9OQPCoc/3tbIq4aoYe/pQHq9jqQRYEW9XBoVx7vkzjOjEjHzH8fQ2VjudvnvFzbjtVSF4QAL109t+4MwzCDmH7eB6Tnvh8da8QImCh7/sMTH8DfjjwLsaWJLABMw8zwr9DYnVM6zgBVob/Toig2AHmzs/Hr7KvcSvtrs79DZPY1oguteCe0YFmxX7ff1UobLFZBdGndj8Mn4GxtXr9dtESQ5KpbNZ0m0ZyZ7+N43q1o0Cu60yi2GNnYsc2A8R4U5X/oURomXvV7CCMfR255mPgiZiKrMA8LMWOk/lPsXrscNqvZs/gWPwxh/mPItYaJr9rsGF/PY+k9JOIxcrXomsTCbQSnJXqEAJEb4DQk1woVpkChVmHW9+9hz6zhMAtmJz8KssKmuFlSeKS//TucXzYWp70NkF702f53rEmL6I9qsPvOf7kvROgqOK/b1rkOqpAHQSkgc9sN2HnVMdi682GB47JE08tTcXTCIbTp2rqDHnH1cAz5TIuSgALZRbrHV0zHt6l5MNa0eFRXfv77iUh6XEBdWK7Ldfi6XKB5GKo+HYlXl++G1excRm+TeHVxq+o4+E12MVHkgd9PxMOPCwjvkX4lIL2UngUwW4bh0y0jsfz+3TBbXF9j2iFdfYSbAbVKge9fnoVZgXsgmM1OfsmmgOJh24A1XNGBWpwLM0guPUhQQdFo4RaeYRhmkNKPk9AdV74it7TOv09scDjHHnWD80ubTUEfiCoqAjCuOAMqhdqtFI+snNkdhqNKOxH+lX06/ElUwQkAVozdcEkVpL2NqsTPZx7HgV23Qa9XiPrpeXLRST+c1n+AKRlPuxWfxkuLWcuP4QhWI7c0TPrSi/yu0KuxU7UEE/63CiHDEtyKz1sbjJlv5uPI1X9DrilMPi7HdQ0EdwoMrhdu6/Yo9OKK9Kgzfj6YXLYN2ZmxMAvic6Ok7pDgUdGYfP5d7LlrJE57t7odJwCc9WrEruVKTDz/GIJHR3hw/4rln6BQKzDl3AJkz8mDTbA5FKD96nb7og9gZOlo+Ib7dojnW8eifZsBRQEFsnkmAIXD9uHqkjho/NRu1hUNHjg2C0Grj+B8WK5T1ZCqDgKAZnUFgpfsxCtVExCTECLpT7jI3hs5//5aDfYem4VnVh9BcFiu6+rV47taqMCS63ei6sgEJIwK6d1DpPPTz0eNsvWTkRmcDYEkOgoGeMO/5nazbIGS/+BaMZFhGIa5JAJEuIgzBUmDT5AI+k8HHgHCrKKx03kVfpv6uMt400Zcg/YSH9EwlF6Ed/IeszfEI9qcDJuuT+PO8Xg2wYhVs/6D1Nh56O+1cXqW1LUzD2H/rrHiMYoZ2QS0tws4cu5ppEx7UDYeL+9AjLmxHDtPJotfbnJhnXX6O1IXCt+78hAakygbn09AKOL+XopdSJK2IOX0gTs2kqs1Enpdt50DnX5oEw5EaDy+Z4JHx8C/4HkcGtImo5YcC905vNyw8/ArWIyQhMhedih0kHFgEQ5EHofrdVE7/i/UFmLKlimInBaJxnW10Ct0Li9d1+8zAQW4cUu6y9T5BHphRfkYVCTvlEyRq/UGCIAx9AgezfPFiMRQj86FjA52h4BAL+SVj8HoHuknTzLQmeFQvyPI+8YXiaND3bukIgL+0NrpiLAecLfPaEDIGBEl2y/QNCyGW3eGYZgrT4BImQCujbflCUsdWkb7FtAc4vxavcXcjMNJu+za1J7t0i8VK13Guyz2j9LtdOpZVOnP2h0zDi20W8PLsSNdX+wHxc5fIv3sFjzoZ8GTyY14dNYR/Gbmv5Aav6DPS7sr3gM7UzwzlTs9WSwCaqwvQRsYLelx0s15yC8LkF/QDC5ESOfxcoMakff9CLW3n3hsCiWS/5aH481a9yw7ciM9rgqGXIgqj6+I/cm7EkJk7gsSjU7l44WY/X/BOU2bGxeyZ80XV1blah2G/bgEKh+1G+E5rI/aGebOlAKRgpSvbbun7EbwjgA0KOtd6lRHSjP3IyDKXzapK/ImoTIg323tSTI5aFWX4/EfI+Ht5psXT55yoldfAHblTYKvQ/oFqfrsYi1vtVCOHzdHws9X7d4ldiiMBM0uNzIxsMpkVWKMbOHuCwjm1p1hGObKEyAkYxVKE+oTgqdaH5FsVQiAIVwveu4fC1cBKnIa9AUAvociMG5IimzcceXpkjuWfNf6upP/T8r+BEEgJ/tAzG5tb1agNi8Q9TtTYNu1HONOfYG7wy14KOMErp6wCp6+IXHH98Qp5zFmwiqo/SLh5TsU48evxPRpFdLigYCqajUmXbVRNLy0+Wvw4wnpXsXkaB3SVGsQdTQd6i/84f1FCOILF2J2wHb4a8THiufrApD2W/G5J2n3v4H9zZHixlan0eSrJMyu34Ox7yyD9q4YaH4ZgMhHZmDGtlcwo6UUSldGkuCGYBHgUR2WP1lO+UgsurDucRwLanOQ1j1CJCBznwLD534OTcBqRM/+CrP3ekFBUq+lBOQF1SDt/dt6cU8LMvc1yebfKlhxwrfYbZ3Y83ib0Iqrn5eeO7ZoTRrKYn6UjD1cl4yWNWnYnh6Ftf5qvBPijV0L46HcPhsqm79o6o0B+Xj0wzTRXLrKPUmIK6ma85c1aQgVSX/XOY26ZKx5LQ3pM6PgH6JGSJQ3Fi6Jx/bds2Ejf9GTAjT5+PC1NPcf1WIJFCQuygC/FYnQ+uHq6mLpLHl54Y/HjoNhGIYZvPTDJPTvOnWN2ER0BQRMsvOvVqi6l+G9s3oJvE979zjvwtRNghI2AF9kfoFFOQtF4z6Ydgahe2OdJqJbAByZ/S1+k32d6Hmz4q7Ho2f+6zSN1gpAGW3G0nIvkEgxvZR5GOacSaIT0bvmQ9p6hGkV8WcDEDqmBad8/4ztR152q4zHgWSnHafOyse2ncmi52bOOYicHydLThoPCLTBciYCLc213ecolV4Ydl0Tymo1ojNvM6O+Qs466bc60YmZsMzbjvMtKqdzQ3ytaHslCi36mm7/Ko0PIl7Vo7xNLZnOaF8zFC9eh7L9WZLxRk1Kw7CXN+HAnCjpwjxLzpPYu1xAL4bO0RnAqUbYzxrW2gipPxSh8s3/oHL3IZibWxE5YwJibr0OjQsnIj+wYxUsr0B/BNatRY2qBU6zbDvDnvVxLXbe8jenZMzc8FvsukUJqQnxQyw+0Id9AJO+VSIfV4nUrh75IhtmFo1D9VMnULbtFHyH+GH8c5Ox55dHYBXaITY7uOdlTCvOQMWfa3Hmu1Joo/wx9cUkFP5iFwTBJjqxe9zJdGwYvccpmSovJW5pGgadpkx0DQLtV5n4x4Icycs1PjMaS7Zb0K467xSnrzUET0a1QScyCX4DiU8wVwK41oNqo/FS4kTTMKgd0t/19Pvuq0zcKJP+zIxobN9kgcp23qnIrdYQRGW2oabOOf12k9CtIlXWAtgUWvxQlYo3v6rE7qOVaG4zY8a4SNw6JwYL4xsRuHJgNtoouHkmxuXtklx7IWfKbMz+Nptbd4ZhmEFMP7wBEZuMeuE4IQ+EQhCKQTiJdlsJcmu24Ld7b4X3abVDV5tgr5S8CH8teUEy5r83PSeZmuRTc3rML7Fn6fA/iHYAEoCq+BxR8QEAT+7JgP/0Csl+bsfRDWIT1glAzQlf+B75G+6cvcOjEnbsmASAiCgzcvbPlDx3767ZiI6WXoGqyahASvpTdsdSZj6EsvMaUck6I65UVnwAQPnxHMRWrRHtQG8wKTFhkf38mgk3P4zyFrVkD6wKBN+3lsqKDwCozN0rLz66wuzzHl1B8vsIkwD/tPuQPfcOlGz6Fs3n69BuaEbptr3YdfvT3eIDAFIeuRU1KpNkDUho0YqKDwDYtfR1jGkJgdRbyBqVASmPXOthwVwYhpW5ZSx2jfsMJRsLYNK3QVdSj52//g4ZhydL5P0Cad/Pxvaxu1H06QmY9CbUHa/HNwtykHJwpl3MPVPQPKxWNKw5D6WgXlMmmvLI0hmy4gMA8nPKcXJNrOh91a5swOLHJnh0tT3l3odSoBBJPwCUl86QFR8AkLO7HGveixV9MCjRgMfun+Be/XdYvcPkNQJpz/lj7u+zsWlHCc43NMPQ0o5tB0tx+0u7Bkx8bL4xE+MKdknes4Yxibhu++BbdpthGIbpdwHiqqkWWzHHPfbN+FF2F/T1Be/CnGgUD7VcgzsmPiR6Xlx5utMQawIAgfDeiSck42u3tOLBfTFozfwaKj+baE7FzFHxlXQEVGZfhVszN16U7IsZvRMmk17Sn9ncjPjh2bJqSek/x+4c/8hFkurHWPiiW+nL/2a1s1rqRDVmnn18E2+UrT7T6AxObPu077VCvwjxCwWmIgGqXz+Fqv3u7W/h9YtUGXEvYOi35bLRR2zTyWbYa/5wFwUj/oI0oWU4chZ8IS5gP6qTVXZxbSOxY95O0XNPr66WTK1B0yB6zohF0nNDDr1odKucP1+dL1kVkuep3LrSvdWx8zrTL9aJ8bKb6V+9RloMzJulcl3/HaYMkaDCr19TYX9B1aBqrDYtysTCkhzxjgkCzEOjkHayFq0WXn6XYRhmsKPqv6DlVuuBiBCRH1ndMLEB1++91mWsXwSvx41YKdqoL/Z5AO/AfphTZuw8WE95ixoAmtQ6FBza7yKXhBdy/gfDQxOxbMoa+J/KhL5cYxe/WK6lxr3r9t2EqLAkVNYVuLQZxDREdb1rAVNftxnAz5ytp86Aapti7fzXtY2FlLrK930TuP5N8U0Revw2KqUzXq20N4Tr/McArZAcXG/e/WHfVlN3V9i9SHUz9awee79wf8PJ2vggAM0i4XUUTN2Xh2TPr/v8KHBDkuQ9VzvKXRFlz5CvNSiWsLZr9lQ51FL7t5iR26JRbD4tem7J1rOIlCi9FoVB9Bzb2DrJkh/2Zj6Wvym/R0fHp1FyXQPN8Gq3r3RvqlBcZ/rFnpYvv5mPNW867AMiuu+FseNT5NINH1Lt8SP7rHUqvsjeO2gaKQHAnlsyMeN4jqQHS9gQXGNQoqC2jlt1hmGYywDFwEUtuPG741jltCoknkiEod3gMtSnD/yhe0leR9s14NBwjAi2t7puG/GIpDDYpVzndm7K6o/j+exr8GS5F7Ji58M86z8ISSuFNrrd5ZKgPf83mwRck/SKx2ZhlyYoq8xymdaq8h9kA240eNkd1vX8LTayTm5uslTie/ynb7dfmlZn8ZJ+SSYA53K+7H+t3Cf12z5Q23++9yiERm+bbIFW7ZIfBlO1u0jiHuvIsN671Q1lBie5W73xrHSaTzbAWdleOFL7mbSBaNKboCSV6GUxC+2i5zR76WQFPSA/z9pxBJJjrq0avcv7UHAp2WR6gLx0op0V5E71EiA//58AjULv+hHs8DD6z96B22DQER+1CiW3pmPGiRz726pHoVvCI/CLFl9knz3HLTrDMMyVLUCEXvzn3N3dFm/Ca2n/wLD9MahpqXEr5hZzM/KSdjtImM52q13Ao8kv2fkfcW6GqFGhDrbhndxnepX7g2f/izU7f42n9sbiL+Ve2Bw8GnVTnkZQ5o8IjTO53LJCqZ/sls0sRpOh3OW5Bhd+9EalgwBRyasnuRF1blhkDSb7aqgzq2Sriu5cyaXXx30QUPlX2z0TIEqbRFgdNchQUSt7vrGyQTY9OmWLGwLKuSCq9pZKi4gmk4RS7aByZ5WLem0T6YKQplmlc1nyLlattZNbjstqWxQNHmlYT6uNQqVzWe1ITjhIjZTr/E9ha3D9EHEYhvnVznIMBsL9fVH66wmIL9kjqe3bo0Ygs1aBb0+d5dacYRjmMqIfhmBJLUop3cySimDTAlatFQ1DmlDoV4KP6jbinYJ/Aac8T8Gfjj+ATapckMV5l+nU8gsrYc2NvwHWU97ixl/yEbTlNPdJiZzXncTmg88CeBYAMD/1MUSceR4GncLJ2BIANFYE9qvNTGRz57LI/++qe1lwYUH27MG0Ce5nsL+WAO3NJoYeBtpQdPIiwxCbni1znW0kE6YguSiDdIF3fDdWNUmeYTVZZS+YvlQvG6O1cxUsqSrmbrWQKiHy4CnWkR6LyyvTV6P43N1eR3J9AJHfgmDxOPKiMw0YaBKGhuDgnCHwP3VYUlwZY8di2rFaHK+t55acYRjmMqMf3oDIb04mIBkCxkFAAgSMgoA4KCzDodINg1fZMEQeHIO52fM7xEcvKag5ioappeKrTp32xuLxdwMAlkY/JGnTrq94tt8Kfcvhv6JxzEuS0q1Fr+i1DgjQRruMPyBguGy3cKC//YDyYK3F2fARZKwxx0+poVkS1mCw2uJsb/cwPIJjRg/iW0paybTWN3oUUpBVIaLmLjjtsHDZ87XRoXDeneaCC7L69Op+Jmvv90KxmW0exebqrYKfJVi2qknLL7i1IaK7e1b2dhSfzRIsWtWl0umWuO+tkO70W69vHdA7KHN0NI5l+sO/okiyYtQkTMHIH8tYfDAMw1ymqPo3+IHbseqNpufxJP4p2kYvC/wDPsZaxJVNRzuch2poxhmRXfhFv6ZvR+Hr+Bn+aFdSXfErVAS0y5eqlA0yPOpq1DcUy8YdGX0V6k9LBx6kNaHaToCYUNUovqtyeMEE1Fbk9WnZBKlMqGpX2xtbPdIXM2sBzh8/OEirqqO66n2/eFCbAuf9xUzrju+RM8dDd7JC+jpnjIVOcllsILDNB1Vu5ediC0e4qBKUi93XFIwmdZXowgzZE8JxLq+232UmXcRVtpiCoVJXiYY7bUI48vsp/bJqagBZPCUB64fXQ1lbJzm/rChpJiZu3guT1cotOMMwzGVKP09CH7jW7MOCd2AeaxRNReDhUVg8/h5Yz3g7pVQAcCTkM5fhv5XYhHvT1kCj8ulV+oYEjpS0A/xDrC5LVay3lAAMDV3kMu7Q8BulLTwBCA84a+c/zKdYPAEARs+4t8+vXbjxhGzm1TOX9k9VdXsGsDt13tXkGDfK4ZT8G5OwBfJzhcJumADx1QI6wz/pjvDoi3uYLuqpIZcCVXGYqF8CkHnvYH5T1sHZ4jDJlxl3X6r0X8ws+j7k0Wsm4cPIcij1dfa3UI/07UicjcRPd7H4YBiGYQEiZ2z0lQHTO7aGfCg+DrxVwH2Na0RHDSl8CG8eecxl2MbjWsTuXYW/hjTh/zL3YO7Y29xOl0rphfnR6yUHq/mHGXqd54qTmdBotNJxq3xwumy27Bhyq9F+Q0Rj5SbJy1rjuxwqde9EWMr//A7JLzqPNzce+czZDu7xe78iDqPn3twndUTZM3CHlwRKL6+LrP8XJ0JMXx+WVUvVPx8me371NSEy6RNg2lLmwirtq0k3Qq9Lz6UBv8loV016nh++vAYan9695L3pdylY25Ask1Cl6C1B6Njd3F3+u8noNMqw6/uty2vg08v0/+7uFDTkJnsmOvptOWrXvHZzGv7qnQeFqVk0naRS4r2YDFz9H97hnGEYhgWI2yJkYLrVnjnwMBRhVlEZZC3XiJqGrZNOodZY6XYcbTUqUE4aZhW9j79EmvFCegkeytyI68bfi3HRmQgLGA6FoESAbzhSR/4CK2Z+gFVxOlT+GCNp8lkDD7hVumLDQCor1LhquvTuyemZP+DcOY3IEmEdnwH+NhzdYz//5dievyNmaLuzXU3AyWofTL7jBHy14W6Vl19QBDLu/CfGP9eIo2NfQV5jsJOfYxtXI9rXLGkUWWwCWu/dgJgps2XjChudhBnbymT9+MvUzdhrru2jewC9ugeOvrQeQyxqSauw2K8ZMz/4veh/M/99H4r9dJIdAUOsWhx96Vs38iD0WkD1xVPE1fyKH/5+DKHtMaK2c53PSTx0YjICwn3dii84wg+r/pmBNxrHI/2Vo2gOlh5aqCR/0RWpAWDmtbFu5/Gtvx+DrT3G6RkAAF4+J1FwYjLC3Ex/xBA//POVDDQWjscrDx9FsFeeZ0pvgN5+fLEiE/eb9kKwWESHXNo0vviz70Tc8RXvcM4wDPNToR/mgAgiRsvAdKu1mJtxPGk3RmVnivZWirW7Xxpf7VWOAaC5SgVb1SioMQpJWISx6LFvWAtgPd3h6uRkmkA4WPaqW/GJ9bULAPbvnIjMqZXQtz6HktMbARDi4m6Av/Yp5OyK6bjqYl2uBEyMP4icfPtx5xZzK6KFt3EOq0TTs68sGkMXVGKqbzbqij5C+fHtMOgq4BsQjqCh8QiNnYyghGvREpCCPGMEdpMA6KTlr6W9FcOL16E89i7JmcHlJjV8/5iF2fV7cf6L11F+4Ae0NxsQEpeA4T9bCK85N+Ng+Fj8KMjXvXCLCXqlj2gXuvbtf2HMgytR+t13MDXpe1krem/ZmfRGJGwuQc2voiWtxz1Lw5EZ/3848/inOH/wBIamxiP+rzdhZ1prZ83racZfSFvCZgt2ejTZ+OKtU+rlOXId8+ZWC1rfjgatEt8Dojp6H5ZXDoUieyqOfVSH/O3lqK8wIDDcF5HxQRg5ORTjrw1CWEoLTBF5IGE3jOjc+E/uwWkKh8VHL9oh8Pi/tPBeOQY7vyuFQW+S78BoteDjt6OxrDP9jmUUFb0PFZVDsTN7Kj7aUIft28tRUWlAeJgv4kcGYfLEUFw7JwgpY1sQEZAHwbb7wkPH3erZt5fZYxY05Dh1bPRMp8LUgudaDuG5iXDYgBHOmzJaANg6PgUDGIZhmEFKny9sSvgW9nsP22+JLWDSJc1g0pCJWN9wGFaL4Ng+werwXTGiHfNL3Rt2sxok1/ZJb1os0/Tm0mIAAA1tSURBVG5aAcRkFGPt7rGycY8TiVsq/p7HSeyS9NgWOjLKAuPJOBj05aJVZfrSUuw7E2O/lbRS9nJLH+u5JfVfRFZLUygx5c1zONAeKR6G45bWcvGkS4uQjD0nsDt6tGh5OIXXs1QFiV5pOitx1S2AMMXzHgIfL4yrfBfHgppd1Kqex7pqgXiNSG4MxfGoNbC0mmVu5Kuk4xNcLDpAwT3OuXCuEoDVRV/ECJIu/kMySzrfVzodlTH77C6hVBWRq46O/z0kEee/T2RAM3q3nX9Xt4QCQKpIeIIAHCmdjqGd6Zeq5i4fJCKXWxgncYl2dPq1Ol8q4eZL2wjRQg8eZFaHh7fUcRYgDMMwg5pLMARrYCmoOYKmqWV2iksqpSdjszzKoTuTZd0dVi0AGBLXhs8Lf+F26Yr1sU+ZddT1ZmYi2z6rVIQhyj9IiI+OEHI3JiF5hE48UMFFIXi47wjZrMh/NBmJfgbxl2py53tQ/Uyfrb9wDqGf9gPpfSCWVhPOTXscMe0+LmqV2Ds95wKKNgehYsYGefFx0Wnv/bJKUsVPLqJ7LykXUbpklzsQ9fyPenmvAkDOepPYaCG5Tcml00JARlIuWnTJTksGkNz55OaDyZ3L1FebmVzsbeLuEmgCRJfo7te9ghiGYZjBLkCEQZXJtYYXJIVBd4OvILx5/DG3w/SNNru10Rm5WUIRY43Y0ZKBGt1Jt0pXao+Eb3dNxpSZxR4pGI2GMHH4Mzi6X37oV7vJgKJNMciMP9ARiNRapO4YAG74aW2qx5nfjcBMFLpnifZixN+hl59Hkqmpj21woU/vB13JORiTnsDkWh8PEum8Td6kuiFoTvoYDcVVvcjPwN3T7tjUbYZ2vBtThKgDmRDo4uNzxXvPH4KqKcmtW8wd295oaMf4mCLkHcjsvrVIKv9ye/F4YoBf/EJt/dd0CBJCQ6oT4mI2Y2EYhmF+CgLkYvYf7h8+zF8L6lyS1zFV3SlLrUZxzVG3w7yvXIOscb+BIvMg/ONMooaCOzs4B0WZEZz5Od47EYKy6sMet9WOC9gQWbFt11hMnfkR/Pxt9vaFSCOeEN+MkUG34eDuZ9yK09zejJz3pyHJdDemx1VAoSB5a1FUFBDGBLdgtiILsdvmuzAsddh1XxImZT+KSZo6z9Wei4tAVitqb5qNcSZDr97UQL5W9Z0IOVmOQ0NXIP2fZzCyzddFQdsfH2EKwsx1FhyJeBG6kvO9tAQv5h52X8DI2ZquMDWb8fq0HJy/OwmRFdOhIIVsNRDdY54AbcsYtGbNxifzYyXjslkJT8yuhdowzqVgcjcPLc1mzJmWg4fuTkJ1xXSAFK71rcMekwSgxToGWUdnY/6qWM8e1YPBgHf1allwcUEH+i0OwzAM47EN20ftx1aID9gfmDkgXTyf/hZ+tuceySHUW6e9gDX7n+x1+MNDEzFn1HKM8E6Dt2kEbLowtDZq0NaiQGuzAEEDKP1tUPpboApvQKvvCRQ1fobvj67x+BKMA0nOI6np0fpqtdGYlPJ/aG2dh+raINTWq6HUEMKGmhEVVY72tg9xaO+zIOr9mvoRI6YgfvpKCGEZqLPGQG/WwGBWosUiwNuH4OdLCPY1IVhVD29LGVrOZKEkZx0aq0/3Kr7YGdciZt4KmOOno9ZnCAyCGnqFEmYB8FETgtUWhFIzAvTlsBXtx9nP/4WKg65Xz1Go1JjyxNNQLVqCqqgo1GnUMKoEkEroMQi/x8B0yTkgZyA5gL0Xc0DE71oBE367EME3TYY+OQANWiuaVG0wCG3wJwEBFiVCjAoEFRig/zwfR1/9L8hq8/BGvgqiA+s9ngNyIf8KADY35oBIPT0OeWhYjpwSgVkr4xGZIUCIqYNVo4dFaYBNaIEXecOL/OBrCoaiPhitZd44mdWC79eVoOq0+zvWq9QKrHx6CuYsUUEbVQVBXQcIRqgEEp0HMsmDPKROicA9K+ORkSEgNqYOGpUeShgg2FpAZm+QzQ+mlmDUNwSjrNwbWTtbsO7DEpwudZ1+2uF8WbvnjSy6tM/m7jkgUhPmHB9yFhf/d+ZLaOIGnmEY5ooRIIMVP40WP/g2wtyocGrrEGbFL3R+MFtNl0VepASIFUA1d/8xDMMwDMMwgxjFlZLRof5RoGbBTnl1UTXu4GUjPhzVI8sNhmEYhmEYhgXIIOTF8etAZkHEgCe8Xfqnyy4/Yjs/MwzDMAzDMAwLkAHEW+WDjOFX45PZ2RiZM13UeLem1mFv6feXbR75DQjDMAzDMAxzOaH6qWXomWn/D0v3P9IxJ8ICWMoAa1nHXEtHo10A8I7p8csuj45LfDIMwzAMwzDM5cJP8g2I44KkYntmCACa009hU8G7l2X+xPYXYDHCMAzDMAzDsAAZYKR2RlaMbMOdxzJ/Uvnj4VgMwzAMwzAMC5ABMsrF9r/reltgS9bhN/pU1BgrfxJ5ZBiGYRiGYZjLhZ/cHBDq8dk9V8KboAizoDG6Elus/8bfDzwNuswHLAkOeWUY5sqGqONJIAjCJTmPYRiGYS7WZmfHjh27K9J1cSnjPH36NBERpaSkdB+LjIzsTsvIkSO7j0+YMIGIiE6dOtUv+XDnPHfDTklJoVdffZWKi4upra2Namtrad++fXT//feTUqm087t06VIqKiqitrY2KioqoltuucWjcDyJix07duzYDS6nYP3FMAxzacnKygIAzJ07t/vYvHnzur8vWLCg+3tmZsdctR07dsiGKQjCgL/FWLt2LY4fP4758+dDq9Vi0qRJKCwsxOuvv45//OMfdnlav349srKyEBkZiR07dmD9+vVIT093Oxx342IYhmH4DQg7duzYXXZvQO688046evQoGQwGKiwspPvuu8/u/+nTp9MPP/xAOp2OmpubaevWrfTzn/9cMrwlS5YQEdHXX3/dfWzz5s1UW1tLlZWVtGPHju7jn3zyCRERLVmyxC699913H5WWlpLVapXMx4oVK6igoIBaW1spPz+fli1b5uSvi9tvv52Ki4uptbWVDh8+TGlpaXb/98ST8g0ICCAiIqPR2H3syy+/JCKi0aNHEwAaM2YMERF9/vnnHoXTGz/s2LFjx25QOC4EduzYsQCR+n/lypVERPT+++9TQEAAPfvss0REdPfdd3f7OXnyJBERLViwgLy9vWnWrFn01VdfSYYZERFBRERNTU2kVCpJrVaTXq+n9957j9566y0ym80UEhJCAKiyspKIiCIiIuzS+95771FAQIBkPm6//XYiItqyZQtFRkZSZGQkbd26VVKArFu3jgIDA+mWW24hIqKCgoI+GaZ23XXXERFRUVFR97HS0lIiIvLy8iIA5O3tTUREZ86c8Sic3vhhx44dO3YsQNixY8duUAuQgoICIiKKj48nABQYGEhERPn5+d1+6uvrqa2tjSZPnkwajcateLvCTUtLo7lz53YLmC4jetmyZTR69GgiIiosLHRKb5cgkcrHsWPHiIgoISGh+1hiYqKkAImKiiIApFKpiIjIYrFctACJi4ujkpISslgsNH/+/O7jra2tduEJgkBERK2trR6F46kfduzYsWPHAoQdO3bsBr0A6TKWHTGZTN1+7rnnHjIYDEREZDabaf/+/TRnzhzZeF977TUiInryySfp5ZdfJqPRSN7e3t1vQzZu3EgrVqwgIqI33njDZXodjxuNRiIiO0Hk5eUlKUDkwuqNALnrrruoqamJzGYz3XbbbXb/Ob4B8fHxkXwDIheOJ37YsWPHjh0LEHbs2LG7LATIqVOn7N4QSDmNRkOTJ0+mJ554goiIysvLZf3feOONRESUlZVFRUVFtHHjxu7/Pv74YzIYDPTRRx8REdGiRYs8FiBdb0BGjRrl1huQvhIgERERtGXLFiIi2rdvn91KX45zQLreziQkJDjNAXEnHHf8sGPHjh07FiDs2LFjd1kJkAceeICIiDZs2EChoaHk7+9P1113HW3durXbz6ZNm2jq1KmkVqtp6tSpRERUXFwsG29QUBBZLBYym81ERLR06dLu/xYvXtw9DMpqtVJoaKjHAqRrDsjmzZtp6NChFBkZSd98802vBEhNTQ0REcXFxcnm6Ve/+hXV1dWRXq+nVatWkSAIov4yMzOJiOitt96i4OBgeuutt8hqtVJGRobb4bgbFzt27NixYwHCjh07doNSgMit9rR48WLat28fNTY2UmNjI23ZsoWuueaa7v+vv/56ysrKoubmZtLpdJSVlUWpqaku4z548CAREbW3t1NQUFD3ca1WS21tbURElJub65ZgEjt+xx13UGFhIRmNRsrNzaU777yze5iYJwJk5cqVVF1d7VKsuaKn31tvvZWKi4vJZDJRcXGx3T4g7oTjSVzs2LFjx25wuZ4bajMMwzA/YZKSkpCfn4+ioiIkJiZygTAMwzADAm9EyDAM8xNl48aNSE1NhUajwZgxY/Daa68BAF544QUuHIZhGIYFCMMwDNO3bNiwAW+//TYMBgP2798PQRCwcOFCfPDBB1w4DMMwzIDBQ7AYhmEYhmEYhrlk8BsQhmEYhmEYhmFYgDAMwzAMwzAMwwKEYRiGYRiGYRiGBQjDMAzDMAzDMIOf/w/zpJ7quaBv/gAAAABJRU5ErkJggg=='