# Fileset [Ptychography.py](https://mdr.nims.go.jp/filesets/e52826ff-5356-4fb8-af6d-32cd29460593/download) ## Creator [Nakazawa Katsuaki](https://orcid.org/0000-0002-6056-5615), [Mitsuishi Kazutaka](https://orcid.org/0000-0002-9361-4057) ## Rights MIT License ## Other metadata [Ptychographical reconstruction code by Wigner Distribution Deconvolution （WDD）](https://mdr.nims.go.jp/datasets/0b2d3761-5a37-4b85-b988-b487431778b6) ## Fulltext '''*********************************************************************Ptychography code for WDD reconstruction and aberration measurement.*********************************************************************The codes were originally developed byKatsuaki Nakazawa(ICYS fellow at NIMS, Japan)e-mail:NAKAZAWA.Katsuaki@nims.go.jpand debugged and modified byKazutaka Mitsuishi (NIMS, Japan)e-mail:MITSUISHI.Kazutaka@nims.go.jpbased on the following papers and books:Pennycook, T. J. et al. Efficient phase contrast imaging in STEM using a pixelated detector.Part 1: experimental demonstration at atomic resolution. Ultramicroscopy 151, 160-167,doi:10.1016/j.ultramic.2014.09.013 (2015).Yang, H. et al. Simultaneous atomic-resolution electron ptychography and Z-contrast imaging oflight and heavy elements in complex nanostructures. Nat Commun 7, 12532,doi:10.1038/ncomms12532 (2016).Rodenburg, J. & Maiden, A. in Springer handbook of microscopy(eds P. W. Hawkes & John C. H. Spence) Ch. 17, (Springer Nature, 2019).Rodenburg, J. M. Advances in Imaging and Electron Physics 87-184 (2008).**********************************************************************Disclaimer*THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDINGBUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE ANDNONINFRINGEMENT. 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.*********************************************************************'''from multiprocessing import Poolfrom multiprocessing import Processimport osimport globimport mathimport stem4d as stem4dimport cv2import numpy as npimport hyperspy.api as hsfrom PIL import Imagefrom scipy import ndimagefrom skimage.restoration import unwrap_phaseimport matplotlib.pyplot as pltimport cupy as cpimport sysimport timefrom datetime import datetimeimport gcfrom multiprocessing import Poolimport importlibimportlib.reload(stem4d)class Ptychograph(stem4d.Image4d): def __init__(self, foldername): super().__init__(foldername) if foldername[-1] != '/': foldername += '/' self.foldername = foldername self.foldername2 = foldername self.fft_image4d ="" self.aberr_image ="" self.modes="" self.aperture_zero ="" self.D0 = 100 self.tifnames='' # This value is needed to be adjusted depending on noise level. # See for example C.M. O'Leary et. al. Ultramicroscopy 221 (2021)113189. self.WDD_epciron = 0.01 def calc_abs_phase(self,image,wrapped=False): abs_image = np.abs(image) if wrapped==True: phase_image=np.angle(image) else: phase_image = unwrap_phase(np.angle(image)) # return abs_image, phase_image def get_strong_g_FFT_image(self, num): FFT_abs = Image.open(self.foldername+'fft_abs.tif') FFT_abs = np.array(FFT_abs) intensity = np.sort(np.unique(FFT_abs))[::-1] # for i in range(num): intensity = np.sort(np.unique(FFT_abs))[::-1] ys, xs = np.where(FFT_abs == intensity[i]) y = ys[0] x = xs[0] FFT_abs[y-3:y+4, x-3:x+4] = 0 G = self.get_FFT_image(y, x, saveimg=True) print("G image at=(y,x):(",y,x,")") plt.imshow(np.abs(G)) plt.show() return def get_FFT_image(self, y, x, saveimg=True): name = self.foldername + "FFTs/line_" + str(y)+ ".npy" lineimage = np.load(name) lineimage = lineimage.transpose(1,2,0) image = lineimage[x, :, :] if saveimg == True: abs_image, phase_image = self.calc_abs_phase(image) Image.fromarray(abs_image).save(self.foldername + 'FFT_abs_image_y=' + str(y) + '_x=' + str(x) + '.tif') Image.fromarray(phase_image).save(self.foldername + 'FFT_phase_image_y=' + str(y) + '_x=' + str(x) + '.tif') return image def FFT4D_CPU_Pool(self, cpu_num=os.cpu_count()-5,saveFFT4DData = False):# now = datetime.now().time() print("*******FFT 4D starts*****************:",now) os.makedirs(self.foldername+'FFTs', exist_ok=True) abs_image = np.zeros((self.y, self.x)) image4d = np.zeros((self.y, self.x, self.ky, self.kx//self.bin)) tifnames = glob.glob(self.foldername + "Lines/*.npy") abs_image = np.zeros((self.y, self.x)) now = datetime.now().time() print("*******Reading line data*****************:",now) for name in tifnames: y = int(name[name.rfind('line_')+5:name.rfind('.npy')]) lineimage = np.load(name) image4d[y, :, :, :] = self.bin_lineimage(lineimage) del(lineimage,tifnames) now = datetime.now().time() print("*******FFT starts*****************:",now) fft_image4d = np.fft.fftshift( np.fft.fft( np.fft.fftshift( np.fft.fft(image4d, axis=1), axes=1), axis=0), axes=0) del(image4d) now = datetime.now().time() print("*******FFT finished. Start saving*****************:",now) input = [(y,fft_image4d[y,:,:,:]) for y in range(self.y)] now = datetime.now().time() print("******* input data ready start Pool *****************:",now) with Pool(cpu_num) as p: p.map(self.Pool_save_fft_images_y, input) p.close() p.terminate() del(input) gc.collect() now = datetime.now().time() print("finish saving fft_ line by line data:",now) now = datetime.now().time() print("Start making fft_abs image with Pool:",now) with Pool(cpu_num) as p: result = p.map(self.Pool_make_fft_abs_image, range(self.y)) p.close() p.terminate() abs_image = np.array(result).reshape((self.y,self.x)) plt.imshow(np.abs(abs_image)) plt.show() now = datetime.now().time() print("finish making fft_abs image with Pool:",now) im = Image.fromarray(abs_image.astype(np.float32)) im.save(self.foldername+'fft_abs.tif') def Pool_save_fft_images_y(self, input): y,fft_images_y =input fft_images_y = fft_images_y.transpose(2,0,1) np.save(self.foldername+'FFTs/line_'+str(y),fft_images_y) def Pool_make_fft_abs_image(self, y): fft_images_y = np.load(self.foldername+'FFTs/line_'+str(y)+".npy") fft_images_y = fft_images_y.transpose(1,2,0) abs_images = np.abs(fft_images_y) abs_image_x = abs_images.sum(2).sum(1) return abs_image_x def measure_center_and_radius(self, y=None, x=None): if y == None: y = self.y//2 if x == None: x = self.x//2 image = self.get_FFT_image(y, x) thresh = np.zeros((self.ky, self.kx//self.bin)) thresh[np.where(image >= np.mean(image))] = 1 Image.fromarray(thresh).save(self.foldername + 'thresh_image.tif') cen_ky, cen_kx = ndimage.measurements.center_of_mass(thresh) radius = math.sqrt(np.sum(thresh)/np.pi) * self.pixelsize_k plt.imshow(thresh) plt.show() self.cen_kx = cen_kx self.cen_ky = cen_ky self.radius = radius print("cen_ky=",cen_ky," cen_kx=",cen_kx,"radius=",radius) return cen_ky, cen_kx, radius def calc_distanceratio_and_rotation(self, diff_ky, diff_kx, diff_y, diff_x, cen_ky=None, cen_kx=None, cen_y=None, cen_x=None, num= 5): if cen_ky == None: cen_ky = self.cen_ky if cen_kx == None: cen_kx = self.cen_kx if cen_y == None: cen_y = self.y//2 if cen_x == None: cen_x = self.x//2 diff1 = math.sqrt((cen_y-diff_y)**2 + (diff_x-cen_x)**2) # diff2 = math.sqrt((cen_ky-diff_ky)**2 + (diff_kx-cen_kx)**2)# distance_ratio = diff2/diff1 angle1 = math.atan2(cen_y - diff_y, diff_x-cen_x) angle2 = math.atan2(cen_ky - diff_ky, diff_kx-cen_kx) rotation = angle2 - angle1 self.distance_ratio = distance_ratio self.rotation = rotation print("Acc=",self.accvol," CL=",self.cl," CLA=",self.clap1," mag=",self.mag," x=",self.x," bin=",self.bin) print("distance_ratio=",distance_ratio, " rotation[rad]=",rotation) print("rotaion[deg]=",rotation*180/np.pi) G = self.get_FFT_image(self.y//2,self.x//2)# print("G image at=(y,x):(",self.y//2,self.x//2,")") plt.imshow(np.abs(G)) plt.show() aperture_circle = self.make_circularfilter_circle(self.cen_ky,self.cen_kx,self.radius) plt.imshow(np.abs(G)+ np.max(np.abs(G))*aperture_circle) plt.show() FFT_abs = hs.load(self.foldername+'fft_abs.tif').data FFT_abs[self.y//2-3:self.y//2+4, self.x//2-3:self.x//2+4] = 0 # intensity = np.sort(np.unique(FFT_abs))[::-1] for i in range(num): intensity = np.sort(np.unique(FFT_abs))[::-1] ys, xs = np.where(FFT_abs == intensity[i]) y = ys[0] x = xs[0] FFT_abs[y-3:y+4, x-3:x+4] = 0 G = self.get_FFT_image(y, x) print("G image at=(y,x):(",y,x,")") plt.imshow(np.abs(G)) plt.show() offset_ky, offset_kx = self.calc_offsets(y, x) print("y,x, offset y,x=",y,x,offset_ky,offset_kx) aperture0 = self.make_circularfilter(self.cen_ky, self.cen_kx, 0, self.radius) aperture_plus_k = self. make_circularfilter(self.cen_ky + offset_ky, self.cen_kx + offset_kx, 0, self.radius) aperture_minus_k = self. make_circularfilter(self.cen_ky - offset_ky, self.cen_kx - offset_kx, 0, self.radius) A_plus = aperture0*aperture_plus_k*np.abs((aperture_minus_k-1)) aperture0_circle = self.make_circularfilter_circle(self.cen_ky, self.cen_kx, self.radius) aperture_plus_k_circle = self. make_circularfilter_circle(self.cen_ky + offset_ky, self.cen_kx + offset_kx, self.radius) aperture_minus_k_circle = self. make_circularfilter_circle(self.cen_ky - offset_ky, self.cen_kx - offset_kx, self.radius) A_plus_circle = aperture0_circle+aperture_plus_k_circle+aperture_minus_k_circle plt.imshow( A_plus_circle*np.max(np.abs(G))+np.abs(G)) plt.show() return distance_ratio, rotation def draw_distanceratio_and_rotation(self, num= 10): print("Acc=",self.accvol," CL=",self.cl," CLA=",self.clap1," mag=",self.mag," x=",self.x," bin=",self.bin) print("distance_ratio=",self.distance_ratio, " rotation[rad]=",self.rotation) G = self.get_FFT_image(self.y//2,self.x//2) print("G image at=(y,x):(",self.y//2,self.x//2,")") plt.imshow(np.abs(G)) plt.show() aperture_circle = self.make_circularfilter_circle(self.cen_ky,self.cen_kx,self.radius) plt.imshow(np.abs(G)+ np.max(np.abs(G))*aperture_circle) plt.show() FFT_abs = hs.load(self.foldername+'fft_abs.tif').data FFT_abs[self.y//2-3:self.y//2+4, self.x//2-3:self.x//2+4] = 0 intensity = np.sort(np.unique(FFT_abs))[::-1] for i in range(num): intensity = np.sort(np.unique(FFT_abs))[::-1] ys, xs = np.where(FFT_abs == intensity[i]) y = ys[0] x = xs[0] FFT_abs[y-3:y+4, x-3:x+4] = 0 G = self.get_FFT_image(y, x) print("G image at=(y,x):(",y,x,")") plt.imshow(np.abs(G)) plt.show() offset_ky, offset_kx = self.calc_offsets(y, x) print("y,x, offset y,x=",y,x,offset_ky,offset_kx) aperture0 = self.make_circularfilter(self.cen_ky, self.cen_kx, 0, self.radius) aperture_plus_k = self. make_circularfilter(self.cen_ky + offset_ky, self.cen_kx + offset_kx, 0, self.radius) aperture_minus_k = self. make_circularfilter(self.cen_ky - offset_ky, self.cen_kx - offset_kx, 0, self.radius) A_plus = aperture0*aperture_plus_k*np.abs((aperture_minus_k-1)) aperture0_circle = self.make_circularfilter_circle(self.cen_ky, self.cen_kx, self.radius) aperture_plus_k_circle = self. make_circularfilter_circle(self.cen_ky + offset_ky, self.cen_kx + offset_kx, self.radius) aperture_minus_k_circle = self. make_circularfilter_circle(self.cen_ky - offset_ky, self.cen_kx - offset_kx, self.radius) A_plus_circle = aperture0_circle+aperture_plus_k_circle+aperture_minus_k_circle plt.imshow( A_plus_circle*np.max(np.abs(G))+np.abs(G)) plt.show() return def show_aberration_figure(self, aberrations): aberr_image = self.make_aberration_image(0,0,aberrations) print("aberr image") plt.imshow(np.angle(aberr_image)) plt.show() def calc_offsets(self, y ,x): cen_y = self.y//2 cen_x = self.x//2 distance = math.sqrt((cen_y - y)**2 + (x - cen_x)**2)*self.distance_ratio theta = math.atan2(cen_y - y, x - cen_x) + self.rotation offset_ky = -1 * distance * math.sin(theta) offset_kx = distance * math.cos(theta) return offset_ky, offset_kx def CTF(self,filename=""): #rotation must be radian #radius must be wrriten in mrad # ky and kx must be written in integer num_of_segment = 2 seg_i_angle = [0,15] seg_o_angle =[0,30] seg_azu_angle = np.array([[0, 90], [180, 270]]) k_CoM = np.array([7.5,15+7.5]) startTime = time.time() CTF = np.zeros((self.y, self.x)).astype(np.complex) for qy in range(self.y): for qx in range(self.x): offset_ky, offset_kx = self.calc_offsets(qy, qx) aperture0 = self.make_circularfilter(self.cen_ky, self.cen_kx, 0, self.radius) aperture_plus_k = self. make_circularfilter(self.cen_ky + offset_ky, self.cen_kx + offset_kx, 0, self.radius) aperture_minus_k = self. make_circularfilter(self.cen_ky - offset_ky, self.cen_kx - offset_kx, 0, self.radius) A_plus = aperture0*aperture_plus_k*np.abs((aperture_minus_k-1)) A_minus = aperture0*aperture_minus_k*np.abs((aperture_plus_k-1)) for segment_index in range(num_of_segment): D_k = self.make_circularfilter(0,0,seg_i_angle[segment_index],seg_o_angle[segment_index],seg_azu_angle[segment_index]) CTF[qy, qx] += np.sum(D_k*A_minus * k_CoM[segment_index]-D_k*A_plus* k_CoM[segment_index]) Image.fromarray(abs_phi).save(self.foldername + 'CTF_'+ filename + '.tif') elapsed_time = time.time()- startTime print("elapsed time:{0}".format(elapsed_time)) def get_phi_from_fft4D(self, qy): phi_qy = cp.zeros(self.x).astype(cp.complex) name = self.foldername + "FFTs/line_" + str(qy)+ ".npy" lineimage = np.load(name) lineimage = lineimage.transpose(1,2,0) for qx in range(self.x): G = lineimage[qx, :, :] #[Qx, ky, kx]ã�®é �ã�ªã�®ã�§ã��Gã�¯[ky,kx] offset_ky, offset_kx = self.calc_offsets(qy, qx) aperture0 = self.make_circularfilter(self.cen_ky, self.cen_kx, 0, self.radius) aperture_plus_k = self. make_circularfilter(self.cen_ky + offset_ky, self.cen_kx + offset_kx, 0, self.radius) aperture_minus_k = self. make_circularfilter(self.cen_ky - offset_ky, self.cen_kx - offset_kx, 0, self.radius) A_plus = aperture0*aperture_plus_k*np.abs((aperture_minus_k-1)) A_minus = aperture0*aperture_minus_k*np.abs((aperture_plus_k-1)) phi_qy[qx] = np.sum(G*A_minus-G*A_plus) if qy == self.y//2 and qx == self.x//2: phi_qy[qx] = np.sum(G*aperture0) return phi_qy def WDD_multi(self, aberrations=np.zeros(12),filename="",cpu_num=os.cpu_count()-5): #rotation must be radian #radius must be wrriten in mrad # ky and kx must be written in integer phi_k = cp.zeros((self.y, self.x)).astype(cp.complex) self.aberr_image = self.make_aberration_image(0,0,aberrations) self.aperture_zero = self.make_circularfilter(self.cen_ky, self.cen_kx, 0, self.radius) now = datetime.now().time() print("*******WDD starts*****************:",now) with Pool(cpu_num) as p: result = p.map(self.get_phi_for_WDD, range(self.y)) p.close() # add this. p.terminate() # add this. phi_k = cp.array(result).reshape((self.y,self.x)) now = datetime.now().time() print("*******Pool calculation finished. FFT by GPU starts*****************:",now) phi = cp.conj(phi_k) phi = cp.fft.ifft2(cp.fft.ifftshift(phi)) phi = cp.flip(phi) abs_phi, phase_phi = self.calc_abs_phase(cp.asnumpy(phi)) Image.fromarray(abs_phi).save(self.foldername2 + 'Phi_abs_WDD_multi'+ filename + '.tif') Image.fromarray(phase_phi).save(self.foldername2 +'Phi_phase_WDD_multi'+ filename + '.tif') plt.figure(figsize=(6,6)) plt.imshow( phase_phi ) plt.show() def get_phi_for_WDD(self, qy): phi_k_qy = np.zeros(self.x).astype(np.complex) name = self.foldername + "FFTs/line_" + str(qy)+ ".npy" lineimage = np.load(name) lineimage = lineimage.transpose(1,2,0) for qx in range(self.x): G = lineimage[qx, :, :] offset_ky, offset_kx = self.calc_offsets(qy, qx) H = np.fft.ifft2(np.fft.ifftshift(G)) aperture_plus_k = self. make_circularfilter(self.cen_ky + offset_ky, self.cen_kx + offset_kx, 0, self.radius) probe_zero = self.aberr_image * self.aperture_zero probe_plus_k = np.roll(self.aberr_image, (round(offset_ky), round(offset_kx)), axis=(0, 1))#offsetå��ã��aberr_imageç�»å��ã��ã�¹ã�¯ã�­ã�¼ã�«ã��ã��ã�� probe_plus_k *= aperture_plus_k# qy,qxã�ªã��ã�»ã��ã��ã��ã��ä½�ç½®ã�®å��å·®é�¢æ�°ã�¨çµ�ã��é�¢æ�°ã�� A = probe_zero * np.conj(probe_plus_k)# Ï�ï¼�K)ï½� Ï�ï¼�ï¼�K+Q) if np.sum(A) != 0: kai_A = np.fft.ifft2(np.fft.ifftshift(A))# kai_phi = np.conj(kai_A)*H/(np.abs(kai_A)**2+self.WDD_epciron) D = np.fft.fftshift(np.fft.fft2(kai_phi))# phi_k_qy[qx] = D[self.ky//2, self.kx//2//self.bin]# if qy == self.y//2 and qx == self.x//2: self.D0 = np.sqrt(D[self.ky//2, self.kx//2//self.bin]) return phi_k_qy def calc_aberrations(self, y, x): aberration = np.array([(x**2 + y**2)/2, (x**2 - y**2)/2, x*y, (x**3 - 3*x*y**2)/3, (3*x**2*y - y**3)/3, (x**3 + x*y**2)/3, (y**3 + x**2*y)/3, (x**4 + y**4 + 2*x**2*y**2)/4, (x**4 + y**4 - 6*x**2*y**2)/4, x**3*y - x*y**3, (x**4 - y**4)/4, (x**3*y + x*y**3)/2]) aberration *= 2*math.pi/(self.wavelength * 1e+9) return aberration def measure_aberration_by_specific_Q_iterative(self, num, Qs, cnt_thresh=5): aberr_mat = np.zeros((0,12)) b_array = np.zeros(0) y_array = np.zeros((0),dtype=np.int) x_array = np.zeros((0),dtype=np.int) aberr_array = np.zeros((12)) beta = np.zeros(12) identity_part = np.zeros((0,num)) offset_ky_kx_list = list() Q_K_list = list() for i in range(num): y = Qs[i][0] x = Qs[i][1] y_array = np.hstack((y_array, y)) x_array = np.hstack((x_array, x)) G = self.get_FFT_image(y, x) offset_ky, offset_kx = self.calc_offsets(y, x) offset_ky_kx_list.append([y,x,offset_ky,offset_kx]) overlap = self.make_overlap_region_with_aberration(offset_ky, offset_kx, np.zeros(12)) aperture_minus_k = self. make_circularfilter(self.cen_ky - offset_ky, self.cen_kx - offset_kx, 0, self.radius) A_plus = overlap*np.abs((aperture_minus_k-1)) kys, kxs = np.where(np.abs(A_plus) != 0) Q_K_list.append([y,x,kys,kxs]) aberr_mat = np.vstack((aberr_mat, self.make_aberr_matrix(offset_ky, offset_kx, kys, kxs))) identity_part_tmp = np.zeros((kys.size, num)) identity_part_tmp[:,i] = 1 identity_part = np.vstack((identity_part, identity_part_tmp)) b_array = np.hstack((b_array, self.make_b_array(G, kys, kxs))) aberr_mat = np.hstack((aberr_mat, identity_part)) aberr_mat_inv = self.make_inv_svd(aberr_mat) print("aberr_mat.shape=",aberr_mat.shape) for c in range(1,13): beta[:c] = 0.5 cnt = 0 while cnt <= cnt_thresh: aberr_array += aberr_mat_inv.dot(b_array.T)[:12]*beta b_array = np.zeros(0) for Q_K in Q_K_list: y, x, kys,kxs= Q_K G = self.get_FFT_image(y, x, saveimg=False) offset_ky, offset_kx = self.calc_offsets(y, x) overlap = self.make_overlap_region_with_aberration(offset_ky, offset_kx, aberr_array) aperture_minus_k = self. make_circularfilter(self.cen_ky - offset_ky, self.cen_kx - offset_kx, 0, self.radius) A_plus = overlap*np.abs((aperture_minus_k-1)) G_diff = G*np.conj(A_plus) b_array = np.hstack((b_array, self.make_b_array(G_diff, kys, kxs))) overlap = self.make_disk_with_aberration(0,0, aberr_array) abs_overlap, phase_overlap = self.calc_abs_phase(overlap,wrapped=True) cnt += 1 overlap = self.make_disk_with_aberration(0,0, aberr_array) abs_G_diff, phase_G_diff = self.calc_abs_phase(overlap,wrapped=True) Image.fromarray(abs_G_diff).save(self.foldername + 'G_diff_abs.tif') Image.fromarray(phase_G_diff).save(self.foldername + 'G_diff_phase.tif') return aberr_array def make_overlap_region_with_aberration(self, dky, dkx, aberr_array): disk_zero = self.make_disk_with_aberration(0, 0, aberr_array) disk_plus = self.make_disk_with_aberration(dky, dkx, aberr_array) disk_minus = self.make_disk_with_aberration(-dky, -dkx, aberr_array) overlap = disk_zero * np.conj(disk_plus) * np.abs(disk_minus-1) overlap -= disk_minus * np.conj(disk_zero) * np.abs(disk_plus-1) return overlap def make_disk_with_aberration(self, dky, dkx, aberr_array): aberr = self.make_aberration_image(dky, dkx, aberr_array, saveimg=False) aperture = self.make_circularfilter(self.cen_ky + dky, self.cen_kx + dkx, 0, self.radius) return aperture * aberr def make_aberration_image(self, dky, dkx, aberr_array, saveimg=False): aberr_array = np.asarray(aberr_array,dtype="float32") aberr_phase = np.zeros((self.ky, self.kx//self.bin)) for ky in range(self.ky): for kx in range(self.kx//self.bin): aberr_phase[ky, kx] = self.calc_aberrations((ky-dky-self.cen_ky )* self.pixelsize_k * 0.001, (kx-dkx-self.cen_kx )* self.pixelsize_k * 0.001).dot(aberr_array) aberr = np.exp(1j * aberr_phase) if saveimg == True: abs_aberr, phase_aberr = self.calc_abs_phase(aberr) Image.fromarray(abs_aberr).save(self.foldername + 'aberr_abs.tif') Image.fromarray(phase_aberr).save(self.foldername + 'aberr_phase.tif') return aberr def make_inv_svd(self, matrix): y, x = matrix.shape s, v, d = np.linalg.svd(matrix) if y > x: v = np.hstack((np.diag(1/v), np.zeros((x, y-x)))) elif y < x: v = np.vstack((np.diag(1/v), np.zeros((x-y, y)))) elif y == x: v = np.diag(1/v) matrix_inv = d.T.dot(v).dot(s.T) return matrix_inv def make_b_array(self, G, kys, kxs): _, phase = self.calc_abs_phase(G) b_array = np.zeros((kys.size)) ind = 0 for ky, kx in zip(kys, kxs): b_array[ind] = phase[ky, kx] ind += 1 return b_array def make_aberr_matrix(self, dky, dkx, kys, kxs): aberr_mat_sizem = kys.size aberr_mat_sizen = 12# å��å·®ä¿�æ�°ã�¯12å��ã�� aberr_mat = np.zeros((aberr_mat_sizem, aberr_mat_sizen)) ind = 0 for ky, kx in zip(kys, kxs): aberr_mat[ind, :] = self.calc_aberrations((ky-self.cen_ky) * self.pixelsize_k * 0.001, (kx-self.cen_kx) * self.pixelsize_k * 0.001) aberr_mat[ind, :] -= self.calc_aberrations((ky-dky-self.cen_ky) * self.pixelsize_k * 0.001, (kx-dkx-self.cen_kx) * self.pixelsize_k * 0.001) ind += 1 return aberr_mat