ModuleDefects

21class Defect(ABC):
22    """
23    Abstract class representing a defect on the optical surface.
24    """
25
26    @abstractmethod
27    def RMS(self):
28        pass
29
30    @abstractmethod
31    def PV(self):
32        pass

34class MeasuredMap(Defect):
35    """
36    Class describing a defect map measured experimentally (for instance using a Hartmann wavefront sensor). Note that the provided image should cover the entire support. 
37    """
38    def __init__(self, Support, Map):
39        self.deformation = Map
40        self. Support = Support
41        rect = Support._CircumRect()
42        X = np.linspace(-rect[0], rect[0], num=self.deformation.shape[0])
43        Y = np.linspace(-rect[1], rect[1], num=self.deformation.shape[1])
44        self.DerivX, self.DerivY = np.gradient(self.deformation, rect/self.deformation.shape)
45        DerivInterpX = scipy.interpolate.RegularGridInterpolator((X, Y), np.transpose(self.DerivX), method="linear")
46        DerivInterpY = scipy.interpolate.RegularGridInterpolator((X, Y), np.transpose(self.DerivY), method="linear")
47        SurfInterp = scipy.interpolate.RegularGridInterpolator((X, Y), np.transpose(deformation), method="linear")
48
49        self.rms = np.std(deformation)
50        self.deformation = deformation
51        self.DerivInterp = lambda x: np.array([DerivInterpX(x[:2]), DerivInterpY(x[:2])])
52        self.SurfInterp = lambda x: SurfInterp(x[:2])
53    def get_normal(self, Point):
54        Grad = self.DerivInterp(Point)
55        dX, dY =  Grad.flatten()
56        norm = np.linalg.norm([dX, dY, 1])
57        dX /= norm
58        dY /= norm
59        return np.array([dX, dY, np.sqrt(1 - dX**2 - dY**2)])
60
61    def get_offset(self, Point):
62        return self.SurfInterp(Point)
63
64    def RMS(self):
65        return self.rms
66
67    def PV(self):
68        pass

 70class Fourrier(Defect):
 71    """
 72    Class describing surface defects quantified by their power spectral density
 73    """
 74    def __init__(self, Support, RMS, slope=-2, smallest=0.01, biggest=10):  # The sizes are the wavelength in mm
 75        rect = Support._CircumRect()
 76        k_max = np.pi / smallest
 77        k_min = np.pi / biggest
 78        ResX = int(2 ** math.ceil(math.log2(k_max * rect[0] / 2)))
 79        ResY = int(2 ** math.ceil(math.log2(k_max * rect[1] / 2)))
 80
 81        # TODO Ensure that we don't get astigmatic pixels because of log/ceil
 82
 83        kX = np.linspace(0, k_max, num=ResX)
 84        kY = np.linspace(0, k_max, num=ResY)
 85        kXX, kYY = np.meshgrid(kX, kY)
 86        D = np.sqrt(kXX**2 + kYY**2)  # Matrix of k norms
 87        mask = (D > k_min) * (D < k_max)
 88
 89        if slope == -2:
 90            integ = np.log(k_max / k_min)
 91        else:
 92            a = slope + 2
 93            integ = (k_max**a - k_min**a) / a
 94        b = np.log(2 / np.pi * RMS**2 / integ)
 95        FFT = np.zeros_like(D)
 96        FFT[mask] = D[mask] ** slope * np.exp(b)
 97        FFT = np.sqrt(FFT)
 98
 99        def tiled(x):
100            return np.block([[np.transpose(np.rot90(x)), x], [np.rot90(x, k=2), np.transpose(np.rot90(x, k=-1))]][::-1])
101
102        phases = np.random.uniform(0, 2 * np.pi, size=FFT.shape)
103        phase = tiled(phases)
104        FFT_tiled = tiled(FFT)
105        FFT_tiled = FFT_tiled * np.exp(1j * phase)
106        kX2 = np.linspace(-k_max, k_max, num=ResX * 2)
107        kY2 = np.linspace(-k_max, k_max, num=ResY * 2)
108        kXX2, kYY2 = np.meshgrid(kX2, kY2)
109        DerivFFTX_tiled = FFT_tiled * 1j * kXX2
110        DerivFFTY_tiled = FFT_tiled * 1j * kYY2
111        deformation = np.fft.irfft2(np.fft.fftshift(FFT_tiled), FFT_tiled.shape) * kX2.shape[0] * (k_max - k_min)
112        coeff = RMS / np.std(deformation)
113        deformation = deformation * coeff
114        DerivX = (
115            np.fft.irfft2(np.fft.fftshift(DerivFFTX_tiled), DerivFFTX_tiled.shape)
116            * kX2.shape[0]
117            * coeff
118            * (k_max - k_min)
119        )
120        DerivY = (
121            np.fft.irfft2(np.fft.fftshift(DerivFFTY_tiled), DerivFFTY_tiled.shape)
122            * kY2.shape[0]
123            * coeff
124            * (k_max - k_min)
125        )
126
127        X = np.linspace(-rect[0], rect[0], num=kX2.shape[0])  # TODO why not divided by 2
128        Y = np.linspace(-rect[1], rect[1], num=kY2.shape[0])
129        self.DerivX = DerivX
130        self.DerivY = DerivY
131
132        DerivInterpX = scipy.interpolate.RegularGridInterpolator((X, Y), np.transpose(self.DerivX), method="linear")
133        DerivInterpY = scipy.interpolate.RegularGridInterpolator((X, Y), np.transpose(self.DerivY), method="linear")
134        SurfInterp = scipy.interpolate.RegularGridInterpolator((X, Y), np.transpose(deformation), method="linear")
135
136        self.rms = np.std(deformation)
137        self.deformation = deformation
138        self.DerivInterp = lambda x: np.array([DerivInterpX(x[:2]), DerivInterpY(x[:2])])
139        self.SurfInterp = lambda x: SurfInterp(x[:2])
140
141    def get_normal(self, Point):
142        Grad = self.DerivInterp(Point)
143        dX, dY =  Grad.flatten()
144        norm = np.linalg.norm([dX, dY, 1])
145        dX /= norm
146        dY /= norm
147        return np.array([dX, dY, np.sqrt(1 - dX**2 - dY**2)])
148
149    def get_offset(self, Point):
150        return self.SurfInterp(Point)
151
152    def RMS(self):
153        return self.rms
154
155    def PV(self):
156        pass

159class Zernike(Defect):
160    def __init__(self, Support, coefficients):
161        self.coefficients = coefficients
162        self.max_order = np.max([k[0] for k in coefficients])
163        self.support = Support
164        self.R = Support._CircumCirc()
165
166    def get_normal(self, Point):
167        x, y, z = Point / self.R
168        val, derivX, derivY = zernike_gradient([x], [y], self.max_order)
169        dX = 0
170        dY = 0
171        for k in self.coefficients:
172            dX += self.coefficients[k] * derivX[k][0][1][0]
173            dY += self.coefficients[k] * derivY[k][0][1][0]
174        dX /= self.R
175        dY /= self.R
176        return np.array([-dX, -dY, 1])
177
178    def get_offset(self, Point):
179        x, y, z = Point / self.R
180        val, derivX, derivY = zernike_gradient([x], [y], self.max_order)
181        Z = 0
182        for k in self.coefficients:
183            Z += self.coefficients[k] * val[k][0][1][0]
184        return Z
185
186    def RMS(self):
187        return np.sqrt(np.sum([i**2 for i in self.coefficients.values()]))
188
189    def PV(self):
190        pass

193def normal_add(N1, N2):
194    """
195    Simple function that takes in two normal vectors of a deformation and calculates the total normal vector if the two deformations were individually applied.
196    """
197    normal1 = N1 / np.linalg.norm(N1)
198    normal2 = N2 / np.linalg.norm(N2)
199    grad1X = -normal1[0] / normal1[2]
200    grad1Y = -normal1[1] / normal1[2]
201    grad2X = -normal2[0] / normal2[2]
202    grad2Y = -normal2[1] / normal2[2]
203    gradX = grad1X + grad2X
204    gradY = grad1Y + grad2Y
205    return np.array([-gradX, -gradY, 1])