Dose comparison / Normalization

[Marcraven]

I am running scans at different orientation angles of the same sample, using Lorentz and structure_factor corrections.

I would like the results from the different scans to be comparable to each other, as if every tetrahedron "sees" the same photon intensity.

Is that the case by default? Does it depend on the deposition method?

If not. Any idea to perform dose normalization afterwards?

import os
from xrd_simulator.motion import RigidBodyMotion
from xrd_simulator.polycrystal import Polycrystal
from xrd_simulator.beam import Beam
from xrd_simulator.detector import Detector
import tifffile
import numpy as np
import time
import copy
import utilities
import json



dir_this_file=os.path.dirname(os.path.abspath(__file__))
beams_directory = os.path.join('artifacts','beams')
detectors_directory = os.path.join('artifacts','detectors')
samples_directory = os.path.join('artifacts','samples')

n_angles=180
angle_step=1

beam_size = 50
width = 300
thickness = 300
grain_size = 1
method = 'centroid' # centroid centroid_with_scintillator project

scan_width = np.sqrt(width**2+thickness**2)+beam_size*3
n_steps=int(scan_width/beam_size)*2
threads = min(n_steps,40)

beam_file = f'23keV_{beam_size}x{beam_size}x2000'
detector_file = '227mm_distance_431mm_448mm_size'
sample_file = f'{beam_size}.0x{width}.0x{thickness}.0_{grain_size}.0um'
run_name = utilities.run_name()
tiffs_directory = os.path.join(dir_this_file,'artifacts','tiffs',run_name,'raw_data')

info_dict={
    'run_name':run_name,
    'beam':beam_file,
    'sample':sample_file,
    'detector':detector_file,
    'method':method,
    'n_angles':n_angles,
    'angle_step[º]':angle_step,
    'n_steps':n_steps,
    'scan_width[um]':scan_width,
    'threads':threads
}
if not os.path.exists(tiffs_directory):
    os.makedirs(tiffs_directory)
with open(os.path.join(tiffs_directory,'info.json'), "w") as file:
    json.dump(info_dict, file)

print(f'Starting run {run_name}')

sample = Polycrystal.load(os.path.join(dir_this_file,samples_directory ,f'{sample_file}.pc'))
beam = Beam.load(os.path.join(dir_this_file,beams_directory ,f'{beam_file}.beam'))
detector = Detector.load(os.path.join(dir_this_file,detectors_directory ,f'{detector_file}.det'))

scan_start = RigidBodyMotion(rotation_axis=np.array([0, 0, 1]),
                         rotation_angle=np.radians(0.0001),
                         translation=np.array([0, -scan_width*0.5, 0]))

scan_motion = RigidBodyMotion(rotation_axis=np.array([0, 0, 1]),
                         rotation_angle=np.radians(angle_step),
                         translation=np.array([0, scan_width, 0]))



for angle in range(48,n_angles):
    start_time = time.time() 
    copied_sample = copy.deepcopy(sample)
    print(f'Sample {angle*angle_step} took {time.time()-start_time} seconds to copy')
    start_time = time.time()
    if angle != 0:
        copied_sample.transform(
            RigidBodyMotion(rotation_axis=np.array([0, 0, 1]),
                            rotation_angle=np.radians(angle*angle_step),
                            translation=np.array([0, 0, 0])),time=1)
    print(f'Sample {angle*angle_step} took {time.time()-start_time} seconds to rotate')
    start_time = time.time()
    copied_sample.transform(scan_start,time=1)
    print(f'Sample {angle*angle_step} took {time.time()-start_time} seconds to translate')
    start_time = time.time() 

    copied_sample.diffract(beam=beam,
                    detector=detector,
                    rigid_body_motion=scan_motion,
                    number_of_processes=threads,
                    number_of_frames=n_steps)
                    # proximity = False,
                    # BB_intersection = True)

    print(f'Angle {angle*angle_step} took {time.time()-start_time} seconds to diffract')

    start_time = time.time()
    diffraction_pattern = np.float32(detector.render(frames_to_render='all',
                                    lorentz=True,
                                    polarization=False,
                                    structure_factor=True,
                                    method=method,#"centroid""centroid_with_scintillator",#"project",
                                    number_of_processes=threads,
                                    verbose=False))
    print(f'Angle {angle*angle_step} took {time.time()-start_time} seconds to render')
    file_dir = os.path.join(tiffs_directory,f'angle_{angle}.tif')
    start_time = time.time()
    tifffile.imwrite(file_dir,                     
                     diffraction_pattern,
                     metadata={'angle':angle*angle_step,
                               'angles':n_angles,
                               'steps':n_steps,
                               'beam':beam_file,
                               'sample':sample_file,
                               'detector':detector_file},
                     imagej=True,
                     maxworkers=1)
    print(f'Angle {angle*angle_step} took {time.time()-start_time} seconds to write')
    detector.frames = [] # We empty this variable in case it was used from a previous diffract event.
    del diffraction_pattern, copied_sample

    
        

[AxelHenningsson]

Hi,

That is a nice little simulation example you got coded up there - very happy to see the api in use. 🥳 🚀

In general, the answer is yes - the intensity is comparable between frames with respect to rendering method (project, centroid or scintillator). It is, of course, not comparable in cases when the flags for polarization, Lorentz, or structure factors are toggled on and off between different frames. But as long as the rendering method stays fixed, the intention of xrd_simulator is to keep the intensity comparable between frames as the sample transforms by a sequence of motions.

The rendered intensity of a single diffraction event (originating from a single scattering unit) is always proportional to two variables, namely scattering_unit.volume and intensity_scaling_factor. The first of these is simply the volume of the clipping between the beam and the scattering tetrahedron. The second of these variables is a scaling factor that is nominally unity (1.0) and changes whenever the Lorentz, polarization, or structure factor flags are turned on. To be very precise, the intensity_scaling_factor is defined by the following function:

    def _get_intensity_factor(
            self,
            scattering_unit,
            lorentz,
            polarization,
            structure_factor):
        intensity_factor = 1.0
        if lorentz:
            intensity_factor *= scattering_unit.lorentz_factor
        if polarization:
            intensity_factor *= scattering_unit.polarization_factor
        if structure_factor:
            if scattering_unit.phase.structure_factors is not None:
                intensity_factor *= (scattering_unit.real_structure_factor **
                                    2 + scattering_unit.imaginary_structure_factor**2)
            else:
                raise ValueError("Structure factors have not been set, .cif file is required at sample instantiation.")

        return intensity_factor

This means that the summed intensity (over the detector plane) of single diffraction spot (if it originates from a single scattering unit) is comparable between frames. (It should be noted that diffraction peaks as seen on the detector surface typically are the result of many scattering units diffracting to the same position. Imagine for instance a crystal grain comprised of a 1000 tetrahedrons all with a close to similar orientation. The observed diffraction peak will feature overlap between different parts of the scattering grain volume. Still, the results are very much comparable in terms of intensity between frames. Just a side note to clarify how things work. 🙂 )

Hope that helps,
Cheers
Axel

[Marcraven]

Thanks Axel!

The code is very well written and easy to follow. I have to tidy up a bit my vectorized version of it, but eventually I will send a pull request, although it might also be interesting to keep the original since the readibility is quite higher...

Any reason why i might see inifinte values on some pixels? Have you tried saving diffracted images in int16 vs float32 vs float64?

Cheers

[AxelHenningsson]

Thanks,

Looking forward to the pull request. 👍 Maintaining duplicates codes will likely be a hassle. A pro tip when writing compressed (unreadable) python code is to leave the equivalent unwrapped loops as comments next to the code to explain for the devs what is happening. Another thought; It will be good to have one or two end to end benchmark scripts to run with your edits so we can see how much speed we are gaining for a typical use case. This would be nice to have also going forward in terms of optimization.

The np.inf is a feature, not a bug 😎 . When the Lorentz factors diverge the default rendering is to set np.inf to make sure the user does not trust values in these regions. (this corresponds to having the diffracting planes constantly in their Bragg condition throughout the motion). In reality you will of course get a finite intensity which is something like: structure factor * lorentz factor * polarization factor * grain volume * motion duration * beam intensity. However, for most practical applications such peaks will be saturated on the frame so it is anyway useless data.

For more info on design choices check out the published companion paper:

xrd_simulator: 3D X-ray diffraction simulation software supporting 3D polycrystalline microstructure morphology descriptions Henningsson, A. & Hall, S. A. (2023). J. Appl. Cryst. 56, 282-292