diff --git a/tests/end_to_end_tests/fast_s3dxrd.py b/tests/end_to_end_tests/fast_s3dxrd.py index c1ba0fa..1ad499b 100644 --- a/tests/end_to_end_tests/fast_s3dxrd.py +++ b/tests/end_to_end_tests/fast_s3dxrd.py @@ -40,7 +40,7 @@ def orientation_density_function(x, q): return 1. / (np.pi**2) # uniform ODF. os.path.join( os.path.dirname(__file__), 'saves'), - 'fast_polycrystal_from_odf') + 'fast_polycrystal_from_odf.pc') if 1: def strain_tensor(x): return np.array( diff --git a/tests/end_to_end_tests/laue_powder.py b/tests/end_to_end_tests/laue_powder.py deleted file mode 100644 index a7a057e..0000000 --- a/tests/end_to_end_tests/laue_powder.py +++ /dev/null @@ -1,82 +0,0 @@ -import numpy as np -from xfab import tools -import matplotlib.pyplot as plt -from scipy.signal import convolve -from xrd_simulator import laue, utils -from xrd_simulator.motion import _RodriguezRotator - -"""Simple simulation of 50 random quartz grains in powder diffraction style only using laue.py -and no spatial functions, i.e not considering grain shapes and the like. This is a check to -see that we have our basic crystal equations under control. -""" - -np.random.seed(5) -U = np.eye(3, 3) -strain_tensor = np.zeros((6,)) -unit_cell = [4.926, 4.926, 5.4189, 90., 90., 120.] -B = utils._epsilon_to_b(strain_tensor, unit_cell) -wavelength = 0.285227 -D = 142938.28756189224 # microns -detector = np.zeros((1024, 1024)) -pixsize = 75 # microns - -x = np.array([1., 0., 0.]) -omega = np.linspace(0., np.pi / 2., 3) -ks = np.array([np.array([[np.cos(om), -np.sin(om), 0], - [np.sin(om), np.cos(om), 0], [0, 0, 1]]).dot(x) for om in omega]) -ks = 2 * np.pi * ks / wavelength - -hklrange = 3 -for ii in range(50): # sample of 10 crystals - print('Crystal no ', ii, 'of total ', 50) - phi1, PHI, phi2 = np.random.rand(3,) * 2 * np.pi - U = tools.euler_to_u(phi1, PHI, phi2) - for hmiller in range(-hklrange, hklrange + 1): - for kmiller in range(-hklrange, hklrange + 1): - for lmiller in range(-hklrange, hklrange + 1): - G_hkl = np.array([hmiller, kmiller, lmiller]) - for i in range(len(ks) - 1): - - G = laue.get_G(U, B, G_hkl) - theta = laue.get_bragg_angle(G, wavelength) - - rotation_axis = np.array([0, 0, 1]) - rotator = _RodriguezRotator(rotation_axis) - rotation_angle = omega[i + 1] - omega[i] - rho_0, rho_1, rho_2 = laue.get_tangens_half_angle_equation( - ks[i], theta, G, rotation_axis) - t1, t2 = laue.find_solutions_to_tangens_half_angle_equation( - rho_0, rho_1, rho_2, rotation_angle) - - for j, s in enumerate([t1, t2]): - if s is not None: - - wavevector = rotator(ks[i], s * rotation_angle) - kprime = G + wavevector - - ang = rotation_angle * s - sin = np.sin(-(omega[i] + ang)) - cos = np.cos(-(omega[i] + ang)) - R = np.array( - [[cos, -sin, 0], [sin, cos, 0], [0, 0, 1]]) - kprime = R.dot(kprime) - khat = kprime / np.linalg.norm(kprime) - sd = D / khat[0] - - yd = khat[1] * sd - zd = khat[2] * sd - - col = (-(yd / pixsize) + - detector.shape[1] // 2).astype(int) - row = (-(zd / pixsize) + - detector.shape[0] // 2).astype(int) - - if col > 0 and col < detector.shape[1] and row > 0 and row < detector.shape[0]: - detector[col, row] += 1 - - -kernel = np.ones((4, 4)) -detector = convolve(detector, kernel, mode='full', method='auto') -plt.imshow(detector, cmap='gray') -plt.title("Hits: " + str(np.sum(detector) / np.sum(kernel))) -plt.show() diff --git a/tests/end_to_end_tests/powder.py b/tests/end_to_end_tests/powder.py index 3de4212..3ffd0d1 100644 --- a/tests/end_to_end_tests/powder.py +++ b/tests/end_to_end_tests/powder.py @@ -47,8 +47,17 @@ translation = np.array([0, 0, 0]) motion = RigidBodyMotion(rotation_axis, rotation_angle, translation) -polycrystal.diffract(beam, detector, motion) -diffraction_pattern = detector.render(frames_to_render=0, method='project') +polycrystal.diffract(beam, + detector, + motion, + proximity=True, + BB_intersection=True) + +diffraction_pattern = detector.render(frames_to_render=0, + method='project', + lorentz=False, + polarization=False, + structure_factor=False) plt.imshow(diffraction_pattern > 0, cmap='gray') plt.title("Hits: " + str(len(detector.frames[0]))) diff --git a/tests/end_to_end_tests/s3dxrd_odf_sample.py b/tests/end_to_end_tests/s3dxrd_odf_sample.py index e8b46c5..7581015 100644 --- a/tests/end_to_end_tests/s3dxrd_odf_sample.py +++ b/tests/end_to_end_tests/s3dxrd_odf_sample.py @@ -12,16 +12,16 @@ parameters = { "detector_distance": 191023.9164, - "detector_center_pixel_z": 256.2345, - "detector_center_pixel_y": 255.1129, - "pixel_side_length_z": 181.4234, - "pixel_side_length_y": 180.2343, - "number_of_detector_pixels_z": 512, - "number_of_detector_pixels_y": 512, + "detector_center_pixel_z": 1024.938, + "detector_center_pixel_y": 1020.4516, + "pixel_side_length_z": 45.35585, + "pixel_side_length_y": 45.058575, + "number_of_detector_pixels_z": 2048, + "number_of_detector_pixels_y": 2048, "wavelength": 0.285227, - "beam_side_length_z": 512 * 200., - "beam_side_length_y": 512 * 200., - "rotation_step": np.radians(10.0), + "beam_side_length_z": 102400.0, + "beam_side_length_y": 102400.0, + "rotation_step": np.radians(2.0), "rotation_axis": np.array([0., 0., 1.0]) } @@ -29,55 +29,53 @@ unit_cell = [4.926, 4.926, 5.4189, 90., 90., 120.] sgname = 'P3221' # Quartz +path_to_cif_file = None def orientation_density_function(x, q): return 1. / (np.pi**2) # uniform ODF. -number_of_crystals = 500 -sample_bounding_cylinder_height = 256 * 180 / 128. -sample_bounding_cylinder_radius = 256 * 180 / 128. +number_of_crystals = 10000 +sample_bounding_cylinder_height = parameters['pixel_side_length_z'] +sample_bounding_cylinder_radius = parameters['pixel_side_length_y'] maximum_sampling_bin_seperation = np.radians(10.0) -# Linear strain gradient along rotation axis. - +# Linear strain gradient along rotation axis. def strain_tensor(x): return np.array( [[0, 0, 0], [0, 0, 0], [0, 0, 0.02 * x[2] / sample_bounding_cylinder_height]]) -# Make the beam much smaller than the sample -# vertices = beam.vertices.copy() -# vertices[:,1:] = 0.180*vertices[:,1:]/np.max(vertices[:,1:]) -# beam.set_beam_vertices(vertices) - - polycrystal = templates.polycrystal_from_odf(orientation_density_function, number_of_crystals, sample_bounding_cylinder_height, sample_bounding_cylinder_radius, unit_cell, sgname, + path_to_cif_file, maximum_sampling_bin_seperation, strain_tensor) path = os.path.join( os.path.join( os.path.dirname(__file__), 'saves'), - 'polycrystal_from_odf') + 'polycrystal_from_odf.pc') polycrystal.save(path, save_mesh_as_xdmf=True) polycrystal = Polycrystal.load(path) -# Full field diffraction. +# Full-field diffraction. polycrystal.diffract( beam, detector, motion, min_bragg_angle=0, max_bragg_angle=None, - verbose=True) + verbose=True, + proximity=True, + BB_intersection=True) + diffraction_pattern = detector.render( frames_to_render=0, lorentz=False, polarization=False, structure_factor=False, - method="centroid", + method="centroid_with_scintillator", verbose=True) pr.disable() @@ -86,6 +84,5 @@ def strain_tensor(x): return np.array( ps.print_stats(15) print("") -diffraction_pattern[diffraction_pattern > 0] = 1 -plt.imshow(diffraction_pattern, cmap='gray') +plt.imshow(diffraction_pattern, cmap='jet') plt.show() diff --git a/tests/end_to_end_tests/saves/fast_polycrystal_from_odf.h5 b/tests/end_to_end_tests/saves/fast_polycrystal_from_odf.h5 index a1c0e21..00b8d86 100644 Binary files a/tests/end_to_end_tests/saves/fast_polycrystal_from_odf.h5 and b/tests/end_to_end_tests/saves/fast_polycrystal_from_odf.h5 differ diff --git a/tests/end_to_end_tests/saves/fast_polycrystal_from_odf.pc b/tests/end_to_end_tests/saves/fast_polycrystal_from_odf.pc index 41dd1e9..67ee7f2 100644 Binary files a/tests/end_to_end_tests/saves/fast_polycrystal_from_odf.pc and b/tests/end_to_end_tests/saves/fast_polycrystal_from_odf.pc differ diff --git a/tests/end_to_end_tests/saves/fast_polycrystal_from_odf.xdmf b/tests/end_to_end_tests/saves/fast_polycrystal_from_odf.xdmf index 5764c8f..c8cbc05 100644 --- a/tests/end_to_end_tests/saves/fast_polycrystal_from_odf.xdmf +++ b/tests/end_to_end_tests/saves/fast_polycrystal_from_odf.xdmf @@ -1 +1 @@ -fast_polycrystal_from_odf.h5:/data0fast_polycrystal_from_odf.h5:/data1fast_polycrystal_from_odf.h5:/data2fast_polycrystal_from_odf.h5:/data3fast_polycrystal_from_odf.h5:/data4fast_polycrystal_from_odf.h5:/data5fast_polycrystal_from_odf.h5:/data6fast_polycrystal_from_odf.h5:/data7fast_polycrystal_from_odf.h5:/data8fast_polycrystal_from_odf.h5:/data9fast_polycrystal_from_odf.h5:/data10fast_polycrystal_from_odf.h5:/data11 +fast_polycrystal_from_odf.h5:/data0fast_polycrystal_from_odf.h5:/data1fast_polycrystal_from_odf.h5:/data2fast_polycrystal_from_odf.h5:/data3fast_polycrystal_from_odf.h5:/data4fast_polycrystal_from_odf.h5:/data5fast_polycrystal_from_odf.h5:/data6fast_polycrystal_from_odf.h5:/data7fast_polycrystal_from_odf.h5:/data8fast_polycrystal_from_odf.h5:/data9fast_polycrystal_from_odf.h5:/data10fast_polycrystal_from_odf.h5:/data11fast_polycrystal_from_odf.h5:/data12 \ No newline at end of file diff --git a/tests/end_to_end_tests/saves/polycrystal_from_odf.h5 b/tests/end_to_end_tests/saves/polycrystal_from_odf.h5 index e226dd5..3367a18 100644 Binary files a/tests/end_to_end_tests/saves/polycrystal_from_odf.h5 and b/tests/end_to_end_tests/saves/polycrystal_from_odf.h5 differ diff --git a/tests/end_to_end_tests/saves/polycrystal_from_odf.pc b/tests/end_to_end_tests/saves/polycrystal_from_odf.pc new file mode 100644 index 0000000..dfd4737 Binary files /dev/null and b/tests/end_to_end_tests/saves/polycrystal_from_odf.pc differ diff --git a/tests/end_to_end_tests/saves/polycrystal_from_odf.xdmf b/tests/end_to_end_tests/saves/polycrystal_from_odf.xdmf index aeb7d56..f8dbd21 100644 --- a/tests/end_to_end_tests/saves/polycrystal_from_odf.xdmf +++ b/tests/end_to_end_tests/saves/polycrystal_from_odf.xdmf @@ -1 +1 @@ -polycrystal_from_odf.h5:/data0polycrystal_from_odf.h5:/data1polycrystal_from_odf.h5:/data2polycrystal_from_odf.h5:/data3polycrystal_from_odf.h5:/data4polycrystal_from_odf.h5:/data5polycrystal_from_odf.h5:/data6polycrystal_from_odf.h5:/data7polycrystal_from_odf.h5:/data8polycrystal_from_odf.h5:/data9polycrystal_from_odf.h5:/data10polycrystal_from_odf.h5:/data11 \ No newline at end of file +polycrystal_from_odf.h5:/data0polycrystal_from_odf.h5:/data1polycrystal_from_odf.h5:/data2polycrystal_from_odf.h5:/data3polycrystal_from_odf.h5:/data4polycrystal_from_odf.h5:/data5polycrystal_from_odf.h5:/data6polycrystal_from_odf.h5:/data7polycrystal_from_odf.h5:/data8polycrystal_from_odf.h5:/data9polycrystal_from_odf.h5:/data10polycrystal_from_odf.h5:/data11polycrystal_from_odf.h5:/data12 \ No newline at end of file diff --git a/tests/end_to_end_tests/voroni_powder.py b/tests/end_to_end_tests/voroni_powder.py deleted file mode 100644 index 8cd99a5..0000000 --- a/tests/end_to_end_tests/voroni_powder.py +++ /dev/null @@ -1,155 +0,0 @@ -import matplotlib.pyplot as plt -import meshio -import numpy as np -from xrd_simulator.polycrystal import Polycrystal -from xrd_simulator.mesh import TetraMesh -from xrd_simulator.phase import Phase -from xrd_simulator.detector import Detector -from xrd_simulator.beam import Beam -from xrd_simulator.motion import RigidBodyMotion -from xfab import tools -import os -import cProfile -import pstats - -np.random.seed(23) - -centroids = [] -for i in range(2, 66): - path = os.path.join( - os.path.dirname(__file__), - '../../extras/voroni_sample/grain_meshes') - path = os.path.normpath(path) - grainmeshfile = os.path.join(path, 'grain' + str(i).zfill(4) + '.xdmf') - mesh = meshio.read(grainmeshfile) - centroids.append(list((1 / 256.) * np.mean(mesh.points, axis=0))) -centroids = np.array(centroids) -sample_diameter = 1.0 -coord, enod = [], [] -k = 0 -for c in centroids: - coord.append([c[0], c[1], c[2]]) - coord.append([c[0] + 0.01, c[1], c[2]]) - coord.append([c[0], c[1] + 0.01, c[2]]) - coord.append([c[0], c[1], c[2] + 0.01]) - enod.append([k, k + 1, k + 2, k + 3]) - k += 3 -coord = np.array(coord) -enod = np.array(enod) -mesh = TetraMesh.generate_mesh_from_vertices(coord, enod) -print(mesh.coord) - -print("") -print('nelm:', mesh.number_of_elements) -print("") - -pixel_size = 5 * sample_diameter / 256. -detector_size = pixel_size * 1024 -detector_distance = 50 * sample_diameter -det_corner_0 = np.array( - [detector_distance, -detector_size / 2., -detector_size / 2.]) -det_corner_1 = np.array( - [detector_distance, detector_size / 2., -detector_size / 2.]) -det_corner_2 = np.array( - [detector_distance, -detector_size / 2., detector_size / 2.]) - -detector = Detector( - pixel_size, - pixel_size, - det_corner_0, - det_corner_1, - det_corner_2) - -# data = os.path.join( os.path.join(os.path.dirname(__file__), 'data' ), 'Fe_mp-150_conventional_standard.cif' ) -unit_cell = [3.64570000, 3.64570000, 3.64570000, 90.0, 90.0, 90.0] -sgname = 'Fm-3m' # Iron -phases = [Phase(unit_cell, sgname)] - -grain_avg_rot = np.max([np.radians(1.0), np.random.rand() * 2 * np.pi]) -euler_angles = grain_avg_rot + \ - np.random.normal(loc=0.0, scale=np.radians(20), size=(mesh.number_of_elements, 3)) -orientation = np.array([tools.euler_to_u(ea[0], ea[1], ea[2]) - for ea in euler_angles]) -element_phase_map = np.zeros((mesh.number_of_elements,)).astype(int) -polycrystal = Polycrystal(mesh, orientation, strain=np.zeros( - (3, 3)), phases=phases, element_phase_map=element_phase_map) - -w = detector_size # full field beam -beam_vertices = np.array([ - [-detector_distance, -w, -w], - [-detector_distance, w, -w], - [-detector_distance, w, w], - [-detector_distance, -w, w], - [detector_distance, -w, -w], - [detector_distance, w, -w], - [detector_distance, w, w], - [detector_distance, -w, w]]) -wavelength = 0.115227 -xray_propagation_direction = np.array([1, 0, 0]) * 2 * np.pi / wavelength -polarization_vector = np.array([0, 1, 0]) -beam = Beam( - beam_vertices, - xray_propagation_direction, - wavelength, - polarization_vector) - -rotation_angle = 90.0 * np.pi / 180. -rotation_axis = np.array([0, 0, 1]) -translation = np.array([0, 0, 0]) -motion = RigidBodyMotion(rotation_axis, rotation_angle, translation) - -print("Diffraction computations:") -pr = cProfile.Profile() -pr.enable() -polycrystal.diffract(beam, detector, motion) -pr.disable() -pr.dump_stats('tmp_profile_dump') -ps = pstats.Stats('tmp_profile_dump').strip_dirs().sort_stats('cumtime') -ps.print_stats(10) -print("") - -print("Detector centroid rendering:") -pr = cProfile.Profile() -pr.enable() -diffraction_pattern1 = detector.render( - frames_to_render=0, - lorentz=False, - polarization=False, - structure_factor=False, - method="centroid") -pr.disable() -pr.dump_stats('tmp_profile_dump') -ps = pstats.Stats('tmp_profile_dump').strip_dirs().sort_stats('cumtime') -ps.print_stats(10) -print("") - -print("Detector project rendering:") -pr = cProfile.Profile() -pr.enable() -diffraction_pattern2 = detector.render( - frames_to_render=0, - lorentz=False, - polarization=False, - structure_factor=False, - method='project') -pr.disable() -pr.dump_stats('tmp_profile_dump') -ps = pstats.Stats('tmp_profile_dump').strip_dirs().sort_stats('cumtime') -ps.print_stats(10) - -# diffraction_pattern[ diffraction_pattern<=0 ] = 1 -# diffraction_pattern = np.log(diffraction_pattern) - -# from scipy.signal import convolve - -# kernel = np.ones((4,4)) -# diffraction_pattern1 = convolve(diffraction_pattern1, kernel, mode='full', method='auto') - -fig, ax = plt.subplots(1, 2) -ax[0].imshow(diffraction_pattern1, cmap='gray') -ax[1].imshow(diffraction_pattern2, cmap='gray') -ax[0].set_title("Fast delta peak rendering") -ax[1].set_title("Full projection rendering") -ax[0].set_xlabel("Hits: " + str(len(detector.frames[0]))) -ax[1].set_xlabel("Hits: " + str(len(detector.frames[0]))) -plt.show()