From f14d8882e6f8eb1b393220907e11879d5c4e9da5 Mon Sep 17 00:00:00 2001
From: Daniel Toyra <dtoyra@star.sr.bham.ac.uk>
Date: Wed, 12 Aug 2015 19:00:33 +0100
Subject: [PATCH] Added function that reads the binary metroPro files
---
pykat/optics/maps.py | 224 ++++++++++++++++++++++++++++++++++++-------
1 file changed, 191 insertions(+), 33 deletions(-)
diff --git a/pykat/optics/maps.py b/pykat/optics/maps.py
index 4e7de57..1c427e6 100644
--- a/pykat/optics/maps.py
+++ b/pykat/optics/maps.py
@@ -399,6 +399,24 @@ class surfacemap(object):
return fig
+ def rescaleData(self,scaling):
+ '''
+ Rescaling data. Assumes that the current map is scaled correctly, i.e.,
+ surfacemap.scaling*surfacemap.data would give the data in meters.
+
+ Input: scaling
+ scaling - float number that representing the new scaling. E.g.,
+ scaling = 1.0e-9 means that surfacemap.data is converted to
+ being given in nano meters.
+ '''
+ if isinstance(scaling, float):
+ self.data[self.notNan] = self.data[self.notNan]*self.scaling/scaling
+ self.scaling = scaling
+ else:
+ print('Error: Scaling factor must be of type float, {:s} found'.
+ format(type(scaling)))
+
+
def rms(self, w=None):
if w is None:
return math.sqrt((self.data[self.notNan]**2).sum())/self.notNan.sum()
@@ -1821,44 +1839,186 @@ def read_map(filename, mapFormat='finesse', scaling=1.0e-9):
# map with surfacemap.plot().
data[isNan] = 0
- return surfacemap(name, maptype, size, center, step, scaling, data, notNan)
+ return surfacemap(name, maptype, (0,0), center, step, scaling, data, notNan)
+
+ if mapFormat == 'metroPro':
+ # Remove '.dat' for output name
+ name = filename.split('.')[0]
+ # Reading file
+ data,hData,x1,y1 = readMetroProData(filename)
+ # xy-stepsize in meters
+ stepSize = (hData['cameraRes'], hData['cameraRes'])
+ # Temporary center
+ center = ((len(data[0,:])-1)/2.0,(len(data[:,0])-1)/2.0)
+ # Keeping track of nan's
+ isNan = np.isnan(data)
+ notNan = isNan==False
+ data[isNan] = 0
+ # Change this:
+ mType = 'phase'
+ fType = 'both'
+ maptype = ' '.join([mType, fType])
+ # Unnecessary size, remove from class later
+ size = data.shape[::-1]
+ smap = surfacemap(name, maptype, size, center, stepSize, 1.0, data, notNan)
+ smap.rescaleData(scaling)
+ smap.recenter()
- # TODO: Add options for reading virgo maps, and .xyz zygo
- # maps (need .xys file for this). Binary ligo-maps?
- # The intensity data is not used to anything here. Remove
- # or add to pykat?
-
+ return smap
+def readMetroProData(filename):
+ '''
+ Reading the metroPro binary data files.
+
+ Translated from 'LoadMetroProData.m' by Hiro Yamamoto.
+ '''
+ f = BinaryReader(filename)
+ # Read header
+ hData = readHeaderMP(f)
+ if hData['format'] < 0:
+ print('Error: Format unknown to readMetroProData()\nfilename: {:s}'.format(filename))
+ return 0
+ # Read phase map data
+ # Skipping header and intensity data
+ f.seek(hData['size']+hData['intNBytes'])
+ # Reading data
+ dat = f.read('int32',size=hData['Nx']*hData['Ny'])
+ # Marking unmeasured data as NaN
+ dat[dat >= hData['invalid']] = np.nan
+ # Scale data to meters
+ dat = dat*hData['convFactor']
+ # Reshaping into Nx * Ny matrix
+ dat = dat.reshape(hData['Ny'], hData['Nx'])
+ # Flipping up/down, i.e., change direction of y-axis.
+ dat = dat[::-1,:]
+ # Auxiliary data to return
+ dxy = hData['cameraRes']
+
+ print(dxy)
+ x1 = dxy * np.arange( -(len(dat[0,:])-1)/2, (len(dat[0,:])-1)/2 + 1)
+ y1 = dxy * np.arange( -(len(dat[:,0])-1)/2, (len(dat[:,1])-1)/2 + 1)
-## def rnm(n,m,rho):
-## '''
-## Based on 'FT_Rnm.m'.
-## Calculates radial part of the Zernike polynomial for given n and m, and rho.
-## n - Radial coordinate
-## m - Azimuthal coordinate
-## rho - Matrix with normalised radial coordinates, i.e., rho[i,j] <= 1.
-## '''
+ return dat, hData, x1, y1
+
+
+def readHeaderMP(f):
+ '''
+ Reads header of the metroPro format files.
-## m = abs(m)
-## Rnm = 0
-## # If n<=25 we use original formula, otherwise recurrence formula as factorials lead
-## # to very large numbers.
-## if n<=25:
-## # Radial term
-## S = int((n-m)/2)
-## for k in range(S+1):
-## a = ((-1)**k)*math.factorial(n-k)/\
-## (math.factorial(k)*math.factorial((n+m)/2.0-k)*math.factorial((n-m)/2.0-k))
-## p = a*rho**(n-2*k)
-## Rnm = Rnm+p
-## else:
-## # Use recurrence formula
-## pass
+ Translated from the readerHeader() function within the
+ 'LoadMetroProData.m' function written by Hiro Yamamoto.
+ '''
+ hData = {}
+ hData['magicNum'] = f.read('uint32')
+ hData['format'] = f.read('int16')
+ hData['size'] = f.read('int32')
+ # Check if the magic string and format are known ones.
+ if not (hData['format'] >= 1 and hData['format']<=3 and
+ hData['magicNum']-hData['format'] == int('881B036E',16)):
+ hData['format'] = -1
+ # Read necessary data
+ hData['invalid'] = int('7FFFFFF8',16)
+ f.seek(60)
+ # Intensity data, which we will skip over.
+ hData['intNBytes'] = f.read('int32')
+ # Top-left coordinates, which are useless.
+ hData['X0'] = f.read('int16')
+ hData['Y0'] = f.read('int16')
+ # Number of data points alng x and y
+ hData['Nx'] = f.read('int16')
+ hData['Ny'] = f.read('int16')
+ # Total data, 4*Nx*Ny
+ hData['phaNBytes'] = f.read('int32')
+ f.seek(218)
+ # Scale factor determined by phase resolution tag
+ phaseResTag = f.read('int16')
+ if phaseResTag == 0:
+ phaseResVal = 4096
+ elif phaseResTag == 1:
+ phaseResVal = 32768
+ elif phaseResTag == 2:
+ phaseResVal = 131072
+ else:
+ phaseResVal = 0
+ f.seek(164)
+ intf_scale_factor = f.read('float')
+ hData['waveLength'] = f.read('float')
+ f.seek(176)
+ obliquity_factor = f.read('float')
+ # Eq. in p12-6 in MetroPro Reference Guide
+ hData['convFactor'] = intf_scale_factor*obliquity_factor*\
+ hData['waveLength']/phaseResVal
+ # Bin size of each measurement
+ f.seek(184)
+ hData['cameraRes'] = f.read('float')
-## return Rnm
-
+ return hData
+
+
+import struct
+
+class BinaryReaderEOFException(Exception):
+ def __init__(self):
+ pass
+ def __str__(self):
+ return 'Not enough bytes in file to satisfy read request'
+
+class BinaryReader:
+ # Map well-known type names into struct format characters.
+ typeNames = {
+ 'int8' :'b',
+ 'uint8' :'B',
+ 'int16' :'h',
+ 'uint16' :'H',
+ 'int32' :'i',
+ 'uint32' :'I',
+ 'int64' :'q',
+ 'uint64' :'Q',
+ 'float' :'f',
+ 'double' :'d',
+ 'char' :'s'}
+
+ def __init__(self, fileName):
+ self.file = open(fileName, 'rb')
+
+ def read(self, typeName, size=None):
+ typeFormat = BinaryReader.typeNames[typeName.lower()]
+ typeFormat = '>'+typeFormat
+ typeSize = struct.calcsize(typeFormat)
+ if size is None:
+ value = self.file.read(typeSize)
+ if typeSize != len(value):
+ raise BinaryReaderEOFException
+ unpacked = struct.unpack(typeFormat, value)[0]
+ else:
+ value = self.file.read(size*typeSize)
+ if size*typeSize != len(value):
+ raise BinaryReaderEOFException
+ unpacked = np.zeros(size)
+ for k in range(size):
+ i = k*typeSize
+ unpacked[k] = struct.unpack(typeFormat,value[i:i+typeSize])[0]
+ return unpacked
+
+ def seek(self, offset, refPos=0):
+ '''
+ offset in bytes and refPos gives reference position, where 0
+ means origin of the file, 1 uses current position and 2 uses
+ the end of the file.
+ '''
+ self.file.seek(offset, refPos)
+
+
+ def __del__(self):
+ self.file.close()
+
+# TODO: Add options for reading virgo maps, and .xyz zygo
+ # maps (need .xys file for this). Binary ligo-maps?
+ # The intensity data is not used to anything here. Remove
+ # or add to pykat?
+
@@ -1870,5 +2030,3 @@ def read_map(filename, mapFormat='finesse', scaling=1.0e-9):
# [map3,x_tilt,y_tilt,offset2]=FT_remove_piston_from_mirror_map(map2,w, display_style);
# map3=FT_invert_mirror_map(map3, invert);
-# Understand the internal coordinate system of the
-# maps/matrices.
--
GitLab