diff --git a/python/coupledwave.py b/python/coupledwave.py deleted file mode 100644 index 3002667..0000000 --- a/python/coupledwave.py +++ /dev/null @@ -1,324 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Thu Sep 6 22:14:36 2018 - -@author: david -""" - -import numpy as np -import matplotlib.pyplot as plt - -#evaluate a plane wave for all points in R = [Rx, Ry, Rz] -# E0 [Ex, Ey, Ez] is the electric field vector -# k is the k vector -# n is the refractive index of the material -# R is a list of coordinates in x, y, and z -def planewave2field(E0, k, n, R): - knr = n * (k[0] * R[0] + k[1] * R[1] + k[2] + R[2]) - exp_2pknr = np.exp(1j * 2 * np.pi * knr) - Ex = E0[0] * exp_2pknr - Ey = E0[1] * exp_2pknr - Ez = E0[2] * exp_2pknr - - return [Ex, Ey, Ez] - -def orthogonalize(E, k): - s = np.cross(E, k) - return np.cross(k, s) - -class layersample: - - #generate a layered sample based on a list of layer refractive indices - # layer_ri is a list of complex refractive indices for each layer - # z is the list of boundary positions along z - def __init__(self, layer_ri, z): - self.n = np.array(layer_ri).astype(np.complex128) - self.z = np.array(z) - - #calculate the index of the field component associated with a layer - # l is the layer index [0, L) - # c is the component (x=0, y=1, z=2) - # d is the direction (0 = transmission, 1 = reflection) - def idx(self, l, c, d): - i = l * 6 + d * 3 + c - 3 - print(i) - return i - - #create a matrix for a single plane wave specified by k and E - def solve(self, k, E): - - #calculate the wavenumber - kn = np.linalg.norm(k) - - #s is the plane wave direction scaled by the refractive index - s = k / kn * self.n[0] - - #allocate space for the matrix - L = len(self.n) - M = np.zeros((6*(L-1), 6*(L-1)), dtype=np.complex128) - - #allocate space for the RHS vector - b = np.zeros(6*(L-1), dtype=np.complex128) - - #initialize a counter for the equation number - ei = 0 - - #calculate the sz component for each layer - self.sz = np.zeros(L, dtype=np.complex128) - for l in range(L): - self.sz[l] = np.sqrt(self.n[l]**2 - s[0]**2 - s[1]**2) - - #set constraints based on Gauss' law - for l in range(0, L): - #sz = np.sqrt(self.n[l]**2 - s[0]**2 - s[1]**2) - - #set the upward components for each layer - # note that layer L-1 does not have a downward component - # David I, Equation 7 - if l != L-1: - M[ei, self.idx(l, 0, 1)] = s[0] - M[ei, self.idx(l, 1, 1)] = s[1] - M[ei, self.idx(l, 2, 1)] = -self.sz[l] - ei = ei+1 - - #set the downward components for each layer - # note that layer 0 does not have a downward component - # Davis I, Equation 6 - if l != 0: - M[ei, self.idx(l, 0, 0)] = s[0] - M[ei, self.idx(l, 1, 0)] = s[1] - M[ei, self.idx(l, 2, 0)] = self.sz[l] - ei = ei+1 - - #enforce a continuous field across boundaries - for l in range(1, L): - sz0 = self.sz[l-1] - sz1 = self.sz[l] - A = np.exp(1j * kn * sz0 * (self.z[l] - self.z[l-1])) - if l < L - 1: - B = np.exp(-1j * kn * sz1 * (self.z[l] - self.z[l+1])) - - - #if this is the second layer, use the simplified equations that account for the incident field - if l == 1: - M[ei, self.idx(0, 0, 1)] = 1 - M[ei, self.idx(1, 0, 0)] = -1 - if L > 2: - M[ei, self.idx(1, 0, 1)] = -B - - b[ei] = -A * E[0] - ei = ei + 1 - - M[ei, self.idx(0, 1, 1)] = 1 - M[ei, self.idx(1, 1, 0)] = -1 - if L > 2: - M[ei, self.idx(1, 1, 1)] = -B - b[ei] = -A * E[1] - ei = ei + 1 - - M[ei, self.idx(0, 2, 1)] = s[1] - M[ei, self.idx(0, 1, 1)] = sz0 - M[ei, self.idx(1, 2, 0)] = -s[1] - M[ei, self.idx(1, 1, 0)] = sz1 - if L > 2: - M[ei, self.idx(1, 2, 1)] = -B*s[1] - M[ei, self.idx(1, 1, 1)] = -B*sz1 - b[ei] = A * sz0 * E[1] - A * s[1]*E[2] - ei = ei + 1 - - M[ei, self.idx(0, 0, 1)] = -sz0 - M[ei, self.idx(0, 2, 1)] = -s[0] - M[ei, self.idx(1, 0, 0)] = -sz1 - M[ei, self.idx(1, 2, 0)] = s[0] - if L > 2: - M[ei, self.idx(1, 0, 1)] = B*sz1 - M[ei, self.idx(1, 2, 1)] = B*s[0] - b[ei] = A * s[0] * E[2] - A * sz0 * E[0] - ei = ei + 1 - - #if this is the last layer, use the simplified equations that exclude reflections from the last layer - elif l == L-1: - M[ei, self.idx(l-1, 0, 0)] = A - M[ei, self.idx(l-1, 0, 1)] = 1 - M[ei, self.idx(l, 0, 0)] = -1 - ei = ei + 1 - - M[ei, self.idx(l-1, 1, 0)] = A - M[ei, self.idx(l-1, 1, 1)] = 1 - M[ei, self.idx(l, 1, 0)] = -1 - ei = ei + 1 - - M[ei, self.idx(l-1, 2, 0)] = A*s[1] - M[ei, self.idx(l-1, 1, 0)] = -A*sz0 - M[ei, self.idx(l-1, 2, 1)] = s[1] - M[ei, self.idx(l-1, 1, 1)] = sz0 - M[ei, self.idx(l, 2, 0)] = -s[1] - M[ei, self.idx(l, 1, 0)] = sz1 - ei = ei + 1 - - M[ei, self.idx(l-1, 0, 0)] = A*sz0 - M[ei, self.idx(l-1, 2, 0)] = -A*s[0] - M[ei, self.idx(l-1, 0, 1)] = -sz0 - M[ei, self.idx(l-1, 2, 1)] = -s[0] - M[ei, self.idx(l, 0, 0)] = -sz1 - M[ei, self.idx(l, 2, 0)] = s[0] - ei = ei + 1 - #otherwise use the full set of boundary conditions - else: - M[ei, self.idx(l-1, 0, 0)] = A - M[ei, self.idx(l-1, 0, 1)] = 1 - M[ei, self.idx(l, 0, 0)] = -1 - M[ei, self.idx(l, 0, 1)] = -B - ei = ei + 1 - - M[ei, self.idx(l-1, 1, 0)] = A - M[ei, self.idx(l-1, 1, 1)] = 1 - M[ei, self.idx(l, 1, 0)] = -1 - M[ei, self.idx(l, 1, 1)] = -B - ei = ei + 1 - - M[ei, self.idx(l-1, 2, 0)] = A*s[1] - M[ei, self.idx(l-1, 1, 0)] = -A*sz0 - M[ei, self.idx(l-1, 2, 1)] = s[1] - M[ei, self.idx(l-1, 1, 1)] = sz0 - M[ei, self.idx(l, 2, 0)] = -s[1] - M[ei, self.idx(l, 1, 0)] = sz1 - M[ei, self.idx(l, 2, 1)] = -B*s[1] - M[ei, self.idx(l, 1, 1)] = -B*sz1 - ei = ei + 1 - - M[ei, self.idx(l-1, 0, 0)] = A*sz0 - M[ei, self.idx(l-1, 2, 0)] = -A*s[0] - M[ei, self.idx(l-1, 0, 1)] = -sz0 - M[ei, self.idx(l-1, 2, 1)] = -s[0] - M[ei, self.idx(l, 0, 0)] = -sz1 - M[ei, self.idx(l, 2, 0)] = s[0] - M[ei, self.idx(l, 0, 1)] = B*sz1 - M[ei, self.idx(l, 2, 1)] = B*s[0] - ei = ei + 1 - - #store the matrix and RHS vector (for debugging) - self.M = M - self.b = b - - #evaluate the linear system - P = np.linalg.solve(M, b) - - #save the results (also for debugging) - self.P = P - - #store the coefficients for each layer - self.Ptrans = [] - self.Prefl = [] - for l in range(L): - if l == 0: - self.Ptrans.append((E[0], E[1], E[2])) - else: - px = P[self.idx(l, 0, 0)] - py = P[self.idx(l, 1, 0)] - pz = P[self.idx(l, 2, 0)] - self.Ptrans.append((px, py, pz)) - - if l == L-1: - self.Prefl.append((0, 0, 0)) - else: - px = P[self.idx(l, 0, 1)] - py = P[self.idx(l, 1, 1)] - pz = P[self.idx(l, 2, 1)] - self.Prefl.append((px, py, pz)) - - #this converts the coefficients into a NumPy array for convenience later on - self.Ptrans = np.array(self.Ptrans) - self.Prefl = np.array(self.Prefl) - - #store values required for evaluation - #store k - self.k = kn - - #store sx and sy - self.s = np.array([s[0], s[1]]) - - self.solved = True - - #evaluate a solved homogeneous substrate - def evaluate(self, X, Y, Z): - - if not self.solved: - print("ERROR: the layered substrate hasn't been solved") - return - - #this code is a bit cumbersome and could probably be optimized - # Basically, it vectorizes everything by creating an image - # such that the value at each pixel is the corresponding layer - # that the pixel resides in. That index is then used to calculate - # the field within the layer - - #allocate space for layer indices - LI = np.zeros(Z.shape, dtype=np.int) - - #find the layer index for each sample point - L = len(self.z) - LI[Z < self.z[0]] = 0 - for l in range(L-1): - idx = np.logical_and(Z > self.z[l], Z <= self.z[l+1]) - LI[idx] = l - LI[Z > self.z[-1]] = L - 1 - - #calculate the appropriate phase shift for the wave transmitted through the layer - Ph_t = np.exp(1j * self.k * self.sz[LI] * (Z - self.z[LI])) - - #calculate the appropriate phase shift for the wave reflected off of the layer boundary - LIp = LI + 1 - LIp[LIp >= L] = 0 - Ph_r = np.exp(-1j * self.k * self.sz[LI] * (Z - self.z[LIp])) - Ph_r[LI >= L-1] = 0 - - #calculate the phase shift based on the X and Y positions - Ph_xy = np.exp(1j * self.k * (self.s[0] * X + self.s[1] * Y)) - - #apply the phase shifts - Et = self.Ptrans[LI] * Ph_t[:, :, None] - Er = self.Prefl[LI] * Ph_r[:, :, None] - - #add everything together coherently - E = (Et + Er) * Ph_xy[:, :, None] - - #return the electric field - return E - - -#This sample code produces a field similar to Figure 2 in Davis et al., 2010 -#set the material properties -depths = [-50, -15, 15, 40] -n = [1.0, 1.4+1j*0.05, 1.4, 1.0] - -#create a layered sample -layers = layersample(n, depths) - -#set the input light parameters -d = np.array([0.2, 0, 1]) -d = d / np.linalg.norm(d) -l = 5.0 -k = 2 * np.pi / l * d -E0 = [0, 1, 0] -E0 = orthogonalize(E0, d) - -#solve for the substrate field -layers.solve(k, E0) - -#set the simulation domain -N = 512 -M = 1024 -D = [-50, 80] - -x = np.linspace(D[0], D[1], N) -z = np.linspace(D[0], D[1], M) -[X, Z] = np.meshgrid(x, z) -Y = np.zeros(X.shape) -E = layers.evaluate(X, Y, Z) - -Er = np.real(E) -I = Er[..., 0] ** 2 + Er[..., 1] **2 + Er[..., 2] ** 2 -plt.figure() -plt.set_cmap("afmhot") -plt.imshow(Er[..., 1]) diff --git a/python/muve-align.py b/python/muve-align.py new file mode 100644 index 0000000..61e435a --- /dev/null +++ b/python/muve-align.py @@ -0,0 +1,161 @@ +# -*- coding: utf-8 -*- +""" +Created on Fri May 4 17:15:07 2018 + +@author: david +""" + +import imagestack +import numpy +import matplotlib.pyplot as plt +import cv2 +import progressbar +import scipy.ndimage +import sys + + +def press(event): + global i + global ax + global fig + if event.key == 'z': + i = i - 1 + if i < 0: + i = 0 + if event.key == 'x': + i = i + 1 + if i == S.shape[0]: + i = S.shape[0] + if event.key == "up": + for n in range(i, S.shape[0]): + warps[n][1, 2] = warps[n][1, 2] + 1 + if event.key == "down": + for n in range(i, S.shape[0]): + warps[n][1, 2] = warps[n][1, 2] - 1 + if event.key == "left": + for n in range(i, S.shape[0]): + warps[n][0, 2] = warps[n][0, 2] + 1 + if event.key == "right": + for n in range(i, S.shape[0]): + warps[n][0, 2] = warps[n][0, 2] - 1 + aligned = cv2.warpAffine(S[i, :, :, :], warps[i], (S.shape[2], S.shape[1]), flags=cv2.INTER_LINEAR + cv2.WARP_INVERSE_MAP); + ax.cla() + ax.imshow(aligned.astype(numpy.uint8)) + ax.set_title('Slice ' + str(i)) + fig.canvas.draw() + sys.stdout.flush + +#apply the alignment adjustments to I based on the list of warp matrices in t +def apply_alignment(I, t): + A = numpy.zeros(I.shape, dtype=numpy.uint8) + for i in range(0, I.shape[0]): + A[i, :, :, :] = cv2.warpAffine(I[i, :, :, :], t[i], (I.shape[2], I.shape[1]), flags=cv2.INTER_LINEAR + cv2.WARP_INVERSE_MAP); + + return A + +#get the transform that aligns image B to image A +def align(A, B, max_power=5): + # Define the motion model + warp_mode = cv2.MOTION_TRANSLATION + + # Specify the number of iterations. + number_of_iterations = 5000; + + # Specify the threshold of the increment + # in the correlation coefficient between two iterations + termination_eps = 1e-10; + + # Define termination criteria + criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, number_of_iterations, termination_eps) + + #attempt to fit the image, increasing the blur as necessary + i = 0 + while i <= max_power: + #blur the image used for fitting + im1_blur = scipy.ndimage.filters.gaussian_filter(A, (2 ** i, 2 ** i), mode='constant') + im2_blur = scipy.ndimage.filters.gaussian_filter(B, (2 ** i, 2 ** i), mode='constant') + + # Run the ECC algorithm. The results are stored in warp_matrix. + # Define 2x3 matrix and initialize the matrix to identity + warp_matrix = numpy.eye(2, 3, dtype=numpy.float32) + try: + (cc, warp_matrix) = cv2.findTransformECC (im1_blur,im2_blur,warp_matrix, warp_mode, criteria) + except: + #print("Error aligning at p = " + str(i)) + i = i + 1 + else: + #print("Successful alignment at p = " + str(i)) + break + #enforce the fact that the x-axis is already aligned + #warp_matrix[0, 2] = 0 + return warp_matrix + + #if i > 0: + # (cc, warp_matrix) = cv2.findTransformECC (A,B,warp_matrix, warp_mode, criteria) + #warp_matrix[0, 2] = 0 + # return warp_matrix + +if len(sys.argv) < 3: + print("usage: muve-align input_mask output_directory") + print("ex: muve-align *.bmp ./aligned") + return + +fmask = sys.argv[1] +out_dir = sys.argv[2] + +if not os.path.isdir(out_dir): + os.mkdir(out_dir) + +#read the image stack for alignment +print("Loading image stack...") +S = imagestack.load(fmask, dtype=numpy.float32) +#convert to grayscale +G = imagestack.rgb2gray(S) + +# Define the motion model +warp_mode = cv2.MOTION_TRANSLATION + +# Define the motion model +warp_mode = cv2.MOTION_TRANSLATION + +# Specify the number of iterations. +number_of_iterations = 5000 + +# Specify the threshold of the increment +# in the correlation coefficient between two iterations +termination_eps = 1e-10 + +# Define termination criteria +criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, number_of_iterations, termination_eps) + +print("\n\nCalculating alignment transformations...") +bar = progressbar.ProgressBar(max_value=G.shape[0]) +#for each image in the stack, calculate a transformation matrix to align it to the previous image +warps = [numpy.eye(2, 3, dtype=numpy.float32)] +for i in range(1, G.shape[0]): + warps.append(align(G[i-1, :, :], G[i, :, :])) + bar.update(i+1) + + +print("\n\nApplying alignment transformation...") +# Define 2x3 matrix and initialize the matrix to identity +for i in range(1, G.shape[0]): + warps[i][0, 2] = warps[i][0, 2] + warps[i-1][0, 2] + warps[i][1, 2] = warps[i][1, 2] + warps[i-1][1, 2] + + + +i = 0 +aligned = cv2.warpAffine(S[i, :, :, :], warps[i], (S.shape[2], S.shape[1]), flags=cv2.INTER_LINEAR + cv2.WARP_INVERSE_MAP) + +fig, ax = plt.subplots() + +fig.canvas.mpl_connect('key_press_event', press) + +ax.imshow(aligned.astype(numpy.uint8)) +ax.set_title('Slice ' + str(i)) +ax.set_xlabel("Press 'z' and 'x' to change slices and arrow keys to align") +plt.show() + +I = apply_alignment(S, warps) +imagestack.save(I, out_dir + "/aligned_", ".bmp") \ No newline at end of file -- libgit2 0.21.4