From f33a2ec958ebb04840d3b655f50de72a2035cab4 Mon Sep 17 00:00:00 2001 From: Cees Bassa Date: Tue, 27 Feb 2018 22:53:13 +0100 Subject: [PATCH] New python-based data acquisition script --- python/acquire.py | 176 +++++++++++++++++++++++++++++++++++++++++++++ python/capture.py | 75 ------------------- python/compress.py | 135 ---------------------------------- 3 files changed, 176 insertions(+), 210 deletions(-) create mode 100644 python/acquire.py delete mode 100644 python/capture.py delete mode 100644 python/compress.py diff --git a/python/acquire.py b/python/acquire.py new file mode 100644 index 0000000..e48176d --- /dev/null +++ b/python/acquire.py @@ -0,0 +1,176 @@ +#!/usr/bin/env python +import sys +import numpy as np +import cv2 +import time +import ctypes +import multiprocessing +from astropy.time import Time +from astropy.io import fits + +# Capture images +def capture(buf,z1,t1,z2,t2,device,nx,ny,nz): + # Array flag + first=True + + # For ever loop + while True: + # Get frames + for i in range(nz): + # Get frame + ret,frame=device.read() + + # Skip lost frames + if ret==True: + # Store time + t=float(time.time()) + + # Convert image to grayscale + gray=np.asarray(cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)).astype(np.uint8) + + # Store buffer + z=gray + + # Store results + if first: + z1[i]=z + t1[i]=t + else: + z2[i]=z + t2[i]=t + + # Assign buffer ready + if first: + buf.value=1 + else: + buf.value=2 + + # Swap flag + first=not first + +def compress(buf,z1,t1,z2,t2,nx,ny,nz): + # Flag to keep track of processed buffer + process_buf=1 + + # Start processing + while True: + # Wait for buffer to become available + while buf.value!=process_buf: + time.sleep(1.0) + + # Process first buffer + if buf.value==1: + t=t1 + z=z1.astype('float32') + process_buf=2 + elif buf.value==2: + t=t2 + z=z2.astype('float32') + process_buf=1 + + # Compute statistics + zmax=np.max(z,axis=0) + znum=np.argmax(z,axis=0) + zs1=np.sum(z,axis=0)-zmax + zs2=np.sum(z*z,axis=0)-zmax*zmax + zavg=zs1/float(nz-1) + zstd=np.sqrt((zs2-zs1*zavg)/float(nz-2)) + + # Convert to float and flip + zmax=np.flipud(zmax.astype('float32')) + znum=np.flipud(znum.astype('float32')) + zavg=np.flipud(zavg.astype('float32')) + zstd=np.flipud(zstd.astype('float32')) + + # Format time + nfd="%s.%03d"%(time.strftime("%Y-%m-%dT%T", time.gmtime(t[0])),int((t[0]-np.floor(t[0]))*1000)) + t0=Time(nfd,format='isot') + dt=t-t[0] + + # Generate fits + fname="%s.fits"%nfd + + # Format header + hdr=fits.Header() + hdr['DATE-OBS']="%s"%nfd + hdr['MJD-OBS']=t0.mjd + hdr['EXPTIME']=dt[-1]-dt[0] + hdr['NFRAMES']=nz + hdr['CRPIX1']=float(nx)/2.0 + hdr['CRPIX2']=float(ny)/2.0 + hdr['CRVAL1']=0.0 + hdr['CRVAL2']=0.0 + hdr['CD1_1']=1.0/3600.0 + hdr['CD1_2']=0.0 + hdr['CD2_1']=0.0 + hdr['CD2_2']=1.0/3600.0 + hdr['CTYPE1']="RA---TAN" + hdr['CTYPE2']="DEC--TAN" + hdr['CUNIT1']="deg" + hdr['CUNIT2']="deg" + hdr['CRRES1']=0.0 + hdr['CRRES2']=0.0 + hdr['EQUINIX']=2000.0 + hdr['RADECSYS']="ICRS" + hdr['COSPAR']=4171 + hdr['OBSERVER']="Cees Bassa" + for i in range(nz): + hdr['DT%04d'%i]=dt[i] + for i in range(10): + hdr['DUMY%03d'%i]=0.0 + + # Write fits file + hdu=fits.PrimaryHDU(data=np.array([zavg,zstd,zmax,znum]),header=hdr) + hdu.writeto(fname) + print("Compressed",fname) + + # Main function +if __name__ == '__main__': + # Settings + devid=int(sys.argv[1]) + nx=720 + ny=576 + nz=250 + + # Initialize device + device=cv2.VideoCapture(devid) + + # Set properties + device.set(3,nx) + device.set(4,ny) + + # Initialize arrays + z1base=multiprocessing.Array(ctypes.c_uint8,nx*ny*nz) + z1=np.ctypeslib.as_array(z1base.get_obj()).reshape(nz,ny,nx) + t1base=multiprocessing.Array(ctypes.c_double,nz) + t1=np.ctypeslib.as_array(t1base.get_obj()) + z2base=multiprocessing.Array(ctypes.c_uint8,nx*ny*nz) + z2=np.ctypeslib.as_array(z2base.get_obj()).reshape(nz,ny,nx) + t2base=multiprocessing.Array(ctypes.c_double,nz) + t2=np.ctypeslib.as_array(t2base.get_obj()) + buf=multiprocessing.Value('i',0) + + pcapture=multiprocessing.Process(target=capture,args=(buf,z1,t1,z2,t2,device,nx,ny,nz)) + pcompress=multiprocessing.Process(target=compress,args=(buf,z1,t1,z2,t2,nx,ny,nz)) + + pcapture.start() + pcompress.start() + pcapture.join() + pcompress.join() + +# print("Starting acquisition") +# for i in range(5): + # Capture +# tstart=time.time() +# t,z=capture(device,nx,ny,nz) +# print("Captured %d %dx%d frames in %.3f s"%(nz,nx,ny,time.time()-tstart)) + + # Compress +# tstart=time.time() +# compress(t,z) +# print("Compressed %d %dx%d frames in %.3f s"%(nz,nx,ny,time.time()-tstart)) + + # Release device + device.release() + + diff --git a/python/capture.py b/python/capture.py deleted file mode 100644 index e641549..0000000 --- a/python/capture.py +++ /dev/null @@ -1,75 +0,0 @@ -#!/usr/bin/env python - -import numpy as np -import cv2 -import time -import subprocess -import os - -# Make directory -if not os.path.exists("/dev/shm/video0"): - os.mkdir("/dev/shm/video0") - -# Get sunset and sunrise times (sattools allnight program) -tsolar=subprocess.check_output(['allnight']).replace("\n","").split(" ") - -# Convert to time structs -tset=time.strptime(tsolar[0]+" UTC", "%Y-%m-%dT%H:%M:%S %Z") -trise=time.strptime(tsolar[1]+" UTC", "%Y-%m-%dT%H:%M:%S %Z") -tnow=time.gmtime() -dtset=time.mktime(tset)-time.mktime(tnow) -dtrise=time.mktime(trise)-time.mktime(tnow) - -# Wait for sunset -#if dtset>0: -# print time.strftime("%FT%T",time.gmtime())+" waiting for sunset (%ds)"%dtset -# time.sleep(dtset) - -# Settings -device=cv2.VideoCapture(0) -device.set(3,720) -device.set(4,576) - -# Set counter -iframe=1 - -# Start capture -print time.strftime("%FT%T",time.gmtime())+" start capture" -while dtrise>0: - # Get frame - ret,frame=device.read() - - # Skip lost frames - if ret==True: - # Get time - t=float(time.time()) - - # Format time - nfd="%s.%03d"%(time.strftime("%Y-%m-%dT%T", time.gmtime(t)),int((t-np.floor(t))*1000)) - - # Get Size - ny,nx=frame.shape[0],frame.shape[1] - - # Convert image to grayscale - gray=np.asarray(cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)).astype(np.int8) - - # Open output file - f=open("/dev/shm/video0/img%06d.pgm"%iframe,"w") - f.write("P5\n# %s\n%d %d\n255\n"%(nfd,nx,ny)) - f.write(gray) - f.close() - - # Log - if iframe%250==1: - print "%06d: %s %dx%d"%(iframe,nfd,nx,ny) - - # Increment - iframe+=1 - - # Refresh time - tnow=time.gmtime() - dtrise=time.mktime(trise)-time.mktime(tnow) - - -# End capture -print time.strftime("%FT%T",time.gmtime())+" end capture" diff --git a/python/compress.py b/python/compress.py deleted file mode 100644 index 1f8be92..0000000 --- a/python/compress.py +++ /dev/null @@ -1,135 +0,0 @@ -#!/usr/bin/env python -import numpy as np -from astropy.time import Time -from astropy.io import fits -import time -import glob -import os - -def read_pgm(fname): - # Open file - f=open(fname,"r") - - # Read lines - line1=f.readline() - line2=f.readline() - line3=f.readline() - line4=f.readline() - - # Read parameters - nfd=line2.split(" ")[1] - nx,ny=int(line3.split(" ")[0]),int(line3.split(" ")[1]) - - # Read image - z=np.fromfile(f,dtype='uint8',count=nx*ny).astype('float32') - - # Close file - f.close() - - return z,nfd - -# Create fits file -def create_fits_file(nx,ny,nz,path,iframe): - # Allocate - z=np.empty(nz*ny*nx,dtype='float32').reshape(nz,nx*ny) - mjd=np.empty(nz,dtype='float64') - dt=np.empty(nz,dtype='float32') - - # Loop over frames - tstart=time.time() - for i in xrange(nz): - z[i],nfd=read_pgm(path+"/img%06d.pgm"%(iframe+i)) - - # Format time - t=Time(nfd,format='isot') - if i==0: - t0=t - mjd[i]=t.mjd - dt[i]=86400.0*(mjd[i]-mjd[0]) - print "Files read in %.3f s"%(time.time()-tstart) - - # Compute statistics - tstart=time.time() - zmax=np.max(z,axis=0) - znum=np.argmax(z,axis=0) - z1=np.sum(z,axis=0)-zmax - z2=np.sum(z*z,axis=0)-zmax*zmax - zavg=z1/float(nz-1) - zstd=np.sqrt((z2-z1*zavg)/float(nz-2)) - print "Statistics in %.3f s"%(time.time()-tstart) - - # Reshape, reformat and flip - tstart=time.time() - zmax=np.flipud(zmax.astype('float32').reshape(ny,nx)) - znum=np.flipud(znum.astype('float32').reshape(ny,nx)) - zavg=np.flipud(zavg.astype('float32').reshape(ny,nx)) - zstd=np.flipud(zstd.astype('float32').reshape(ny,nx)) - z=np.array([zavg,zstd,zmax,znum]) - print "Reshape in %.3f s"%(time.time()-tstart) - - # Filename - fname="%s.fits"%t0 - - # Format header - hdr=fits.Header() - hdr['DATE-OBS']="%s"%t0 - hdr['MJD-OBS']=t0.mjd - hdr['EXPTIME']=dt[-1]-dt[0] - hdr['NFRAMES']=nz - hdr['CRPIX1']=float(nx)/2.0 - hdr['CRPIX2']=float(ny)/2.0 - hdr['CRVAL1']=0.0 - hdr['CRVAL2']=0.0 - hdr['CD1_1']=1.0 - hdr['CD1_2']=0.0 - hdr['CD2_1']=0.0 - hdr['CD2_2']=1.0 - hdr['CTYPE1']="RA---TAN" - hdr['CTYPE2']="DEC--TAN" - hdr['CUNIT1']="deg" - hdr['CUNIT2']="deg" - hdr['CRRES1']=0.0 - hdr['CRRES2']=0.0 - hdr['EQUINIX']=2000.0 - hdr['RADECSYS']="ICRS" - hdr['COSPAR']=4171 - hdr['OBSERVER']="Cees Bassa" - for i in xrange(nz): - hdr['DT%04d'%i]=dt[i] - for i in xrange(10): - hdr['DUMY%03d'%i]=0.0 - - # Write fits file - hdu=fits.PrimaryHDU(data=z,header=hdr) - hdu.writeto(fname,clobber=True) - - return fname - -# Main function -if __name__ == '__main__': - # Settings - nx=720 - ny=576 - nz=250 - path="/dev/shm/video0" - - # Start forever loop - while True: - # Find files - files=sorted(glob.glob(path+"/img??????.pgm")) - - # Enough files - if len(files)>nz: - # Get first frame - iframe=int(files[0].replace(path+"/img","").replace(".pgm","")) - - tstart=time.time() - fname=create_fits_file(nx,ny,nz,path,iframe) - tend=time.time() - print "Created %s in %.2f"%(fname,tend-tstart) - - # Remove files - for i in xrange(nz): - os.remove(path+"/img%06d.pgm"%(iframe+i)) - else: - time.sleep(1)