Skip to content

Latest commit

 

History

History
478 lines (370 loc) · 22.6 KB

File metadata and controls

478 lines (370 loc) · 22.6 KB

User Manual: Diffraction Simulation Tools

Version: 11.0
Contact: hsharma@anl.gov


MIDAS provides two complementary simulation tools for High Energy Diffraction Microscopy (HEDM):

Tool Modality Output Use Case
ForwardSimulationCompressed Far-Field (FF) Zarr/ZIP detector images + spot CSV Simulating FF diffraction patterns with intensity, Gaussian blur, energy spread
simulateNF Near-Field (NF) Zarr/ZIP hit-count images + binary files Simulating NF diffraction spot coverage for reconstruction validation

1. Far-Field Simulation: ForwardSimulationCompressed

ForwardSimulationCompressed is a high-performance command-line tool for simulating far-field HEDM experiments. It takes a description of a crystalline sample (grain orientations, positions, and strains) and experimental geometry as input, and produces simulated detector images as output.

The simulation is highly parallelized for speed and includes advanced features to accurately model physical effects such as:

  • Detector distortions and tilts.
  • Wavelength/energy spread (Δλ/λ).
  • Strain within crystal lattices.
  • Relative peak intensities from powder diffraction data.
  • Sub-pixel spot positioning for high-fidelity images.

The output is a set of compressed image files in a Zarr/ZIP format, suitable for analysis or direct comparison with experimental data.

graph TD
    subgraph "Inputs"
        PF[Parameter File]
        GD["Grain Data<br>(Grains.csv, .vtk, etc.)"]
        HKL[hkls.csv]
        POS["positions.csv<br>(Optional)"]
    end

    subgraph "Simulation"
        FS[ForwardSimulationCompressed]
        PF --> FS
        GD --> FS
        HKL --> FS
        POS --> FS
    end

    subgraph "Outputs"
        ZIP["Zarr Zip Files<br>(scanNr_X.zip)"]
        SM["SpotMatrixGen.csv<br>(Optional)"]
        FS --> ZIP
        FS --> SM
    end
Loading

1.1. How to Run

The program is executed from the command line with two arguments: a parameter file and the number of CPU cores to use.

Command Syntax:

./ForwardSimulationCompressed <ParameterFile> <nCPUs>

Arguments:

  • <ParameterFile>: The path to a text file containing all simulation parameters. Details are in Section 1.2.
  • <nCPUs>: The number of CPU cores (threads) to use for the simulation. For best performance, this should be set to the number of physical cores on your machine.

Example:

./ForwardSimulationCompressed experiment_params.txt 16

This command runs the simulation using the settings in experiment_params.txt and parallelizes the computation across 16 CPU cores.

1.2. The Parameter File

The parameter file is a plain text file that controls every aspect of the simulation. Each line consists of a keyword followed by one or more values, separated by spaces. Lines starting with % or # can be used for comments.

1.2.1. Required Files and Their Formats

Before running, you must have several files in your working directory.

A. hkls.csv (Required) This file lists the crystallographic planes (HKLs) to be simulated for your material's crystal structure. It is generated by the GetHKLList utility (which is run automatically by the program). The format is a space-separated text file with a one-line header.

  • Format:
    h k l Intensity D-spacing RingNr
    1 1 1 1.00 2.1105 1
    2 0 0 0.49 1.8277 2
    ...
    
  • The simulation uses the h, k, l, and RingNr columns.

B. positions.csv (Optional, for Multi-Scan) If you are simulating a multi-scan experiment (controlled by the nScans parameter), this file must be present. It specifies the sample translation offset for each scan.

  • Format: A single column of numbers, with one value per scan. The units should match the units of the grain positions in your input file (typically microns).
    0.0
    50.0
    100.0
    ...
    

C. Input Grain Data File (Required) This file contains the core information about your sample's microstructure. The program is flexible and can read several formats, specified with the InFileName parameter.

  • Supported Formats:
    • Grains.csv: A standard format from reconstruction software.
    • EBSD File: A text file with columns for X Y Z Eul1 Eul2 Eul3. The Euler angles should be in degrees.
    • VTK File: A .vtk file from a finite element simulation (e.g., ABAQUS). The code is designed to parse specific fields for position, orientation, and strain.
    • Binary Format (.bin): A custom, high-speed binary format containing the orientation and strain data.

1.2.2. Parameter Keywords

Below is a complete list of keywords that can be used in the parameter file.

File I/O:

  • InFileName <filename>: Path to the input file containing grain data.
  • OutFileName <prefix>: The base name for the output ZIP files. The final output will be <prefix>_scanNr_X.zip.
  • IntensitiesFile <filename> (Optional): Path to a file containing relative peak intensities. If not provided, all peaks will have equal intensity. The format must match mpxrd_hkls.txt (header + 6 columns: h k l SpotInt D-spacing RingNr).

Experimental Geometry:

  • Lsd <value>: Sample-to-detector distance (in microns).
  • Wavelength <value>: Wavelength of the X-ray beam (in Angstroms).
  • OmegaStart <value>: Starting omega angle for the simulation (in degrees).
  • OmegaEnd <value>: Ending omega angle for the simulation (in degrees).
  • OmegaStep <value>: The omega step size per frame (in degrees).

Detector Parameters:

  • NrPixelsY <value>: Horizontal detector dimension in pixels (e.g., 2048).
  • NrPixelsZ <value>: Vertical detector dimension in pixels (e.g., 2048).
  • NrPixels <value>: (deprecated) Square detector shorthand — sets both NrPixelsY and NrPixelsZ to the same value.
  • px <value>: The physical size of a single pixel (in microns).
  • BC <y_center> <z_center>: The pixel coordinates of the beam center on the detector.

Detector Distortion Correction (Optional): These parameters model the geometric distortions of a real detector.

  • tx <value>, ty <value>, tz <value>: Detector tilt angles (in degrees).
  • RhoD <value>, p0 <value>, p1 <value>, p2 <value>: Coefficients for the radial distortion model.

Material and Simulation Parameters:

  • LatticeConstant <a> <b> <c> <alpha> <beta> <gamma>: The lattice parameters of the crystal structure. For cubic, all lengths are equal and angles are 90.
  • PeakIntensity <value>: The maximum intensity value for a single diffraction spot before normalization.
  • GaussWidth <value>: The standard deviation of the Gaussian used to blur the spots, in units of pixels. Controls the spot size.
  • EResolution <spread> <num_samples> (Optional): Simulates the wavelength spread.
    • <spread>: The fractional spread Δλ/λ (e.g., 0.001 for 0.1%).
    • <num_samples>: The number of discrete wavelengths to sample across the spread (an odd integer like 5 or 7 is recommended).
  • nScans <value>: The number of scans to simulate. If greater than 1, a positions.csv file is required.

Output Control:

  • WriteSpots <0_or_1>: If set to 1, a SpotMatrixGen.csv file will be generated, logging the detailed information for every simulated diffraction spot. 0 disables this.
  • WriteImage <0_or_1>: If set to 1 (default), the simulated detector images are written to a Zarr/ZIP file. If set to 0, image generation is skipped, which can save significant time and disk space if only SpotMatrixGen.csv is needed.
  • SimulationBatches <N> (Optional): Controls memory usage for image generation.
    • 0 (default): Auto-detect based on system RAM. If the full image array fits in 75% of system memory, uses a single pass (fastest). Otherwise, automatically selects enough batches to fit.
    • 1: Force single-pass mode — allocates the full image array in memory. Fastest, but requires ~22 GB for a 2048×2048 detector with 1440 frames.
    • N > 1: Process frames in N batches using a two-pass approach. Each batch allocates only (nFrames/N) × NrPixels² floats. For example, SimulationBatches 72 reduces peak memory from ~22 GB to ~320 MB.

1.2.3. Example Parameter File

# ==================================
# Example Parameters for Ti-7Al
# ==================================

# --- File Paths ---
InFileName        meso_grain_file.vtk
OutFileName       Ti7Al_simulation_results
IntensitiesFile   ti_alpha_intensities.txt

# --- Experimental Geometry ---
Lsd               800000.0   # 800 mm in microns
Wavelength        0.1653     # Angstroms
OmegaStart        0.0
OmegaEnd          180.0
OmegaStep         0.5

# --- Detector ---
NrPixelsY         2048
NrPixelsZ         2048
px                200.0      # 200 micron pixels
BC                1024.5 1023.0

# --- Distortion (example values) ---
tx                -0.12
ty                0.05
tz                0.01
RhoD              180000.0
p0                -0.0001
p1                0.00005
p2                0.0002

# --- Material & Simulation ---
LatticeConstant   2.95 2.95 4.68 90 90 120 # Hexagonal Ti-alpha
PeakIntensity     2000.0
GaussWidth        1.5
EResolution       0.001 5 # 0.1% spread, sampled with 5 points

# --- Multi-Scan ---
nScans            1

# --- Output Control ---
WriteSpots        1

1.3. Output Format

The primary output of the simulation is one or more ZIP files, one for each scan. These ZIP files are structured in the Zarr format, a modern standard for chunked, compressed n-dimensional arrays. This format is highly efficient and easily read by many scientific data analysis packages, especially in Python.

  • File Naming: OutFileName_scanNr_X.zip (e.g., Ti7Al_simulation_results_scanNr_0.zip)
  • Internal Structure: Inside the ZIP file, the data is organized as exchange/data/. The detector images are stored as individual, compressed files, one for each omega step (e.g., 0.0.0, 1.0.0, etc.).
  • Data Type: The image data is stored as 16-bit unsigned integers (uint16), normalized so that the brightest pixel in the entire scan corresponds to an intensity of 15000.

To read this data, you can use Python libraries such as zarr. The file structure is a standard Zarr v2 group stored within a Zip file.

Python Example: Reading and Plotting Data

import zarr
import matplotlib.pyplot as plt
import numpy as np

# Path to the output zip file
zip_path = 'Ti7Al_simulation_results_scanNr_0.zip'

# Open the Zip file as a Zarr store
store = zarr.ZipStore(zip_path, mode='r')

# Open the root group
root = zarr.group(store=store)

# Access the data array
# Structure is exchange/data
data = root['exchange']['data']

print("Data shape:", data.shape)
print("Data type:", data.dtype)

# Display the first frame (Omega step 0)
plt.imshow(data[0, :, :], cmap='viridis', vmax=100)
plt.title("Simulated Diffraction Pattern - Frame 0")
plt.colorbar()
plt.show()

# Close the store
store.close()

1.4. Technical Implementation Details

1.4.1. Physics & Rendering Engine

  • Analytic Integration: Instead of simple point sampling, the simulation uses an analytic Gaussian integration method. It models each diffraction spot as a 2D Gaussian distribution and calculates the exact integral of this distribution over the area of each pixel. This ensures high-fidelity simulations even for small spot sizes or coarse pixel grids, eliminating aliasing artifacts.
  • Energy Resolution: The code simulates the effect of a polychromatic beam (energy bandwidth) by discretizing the wavelength spread into multiple sub-samples (num_samples). Each reflection is calculated for these varying wavelengths, and their intensities are averaged, effectively broadening the spots radially.
  • Wedge Correction: The simulation accounts for the "wedge" effect (sample rotation during exposure) by calculating the intersection of the diffraction cone with the detector plane at the start and end angles of the frame.

1.4.2. Data Compression (Zarr/Zip)

  • Custom Writer: The program implements a lightweight, dependency-free Zarr writer in C using the libzip and c-blosc libraries.
  • Storage Layout: It creates a standard Zarr v2 hierarchy (.zgroup, .zattrs, .zarray) inside a ZIP container.
  • Compression: Image data is compressed using the Zstd algorithm with Bitshuffle (via Blosc), providing a high compression ratio for sparse diffraction images while maintaining fast read speeds.

1.4.3. Parallelization and Memory Management

  • OpenMP: The program is parallelized at two levels:
    1. Distortion Map Generation: The pixel-level distortion/tilt correction map (NrPixelsZ × NrPixelsY) is computed in parallel using collapse(2) to distribute the double loop across threads. All per-pixel variables are stack-local, so there are no data races.
    2. Main Simulation Loop: The core simulation loop is parallelized over the input grains/voxels. Each thread calculates the diffraction spots for a subset of grains, applies geometric corrections, and either atomically adds intensity to the shared image array (single-pass mode) or buffers lightweight spot records (batched mode).
  • Per-Thread Spot Buffering: When WriteSpots 1 is enabled, each thread buffers its spot data in a private, dynamically-growing array instead of writing to the output file inside a critical section. After the parallel region completes, all buffers are flushed sequentially. This eliminates lock contention that would otherwise serialize all threads.
  • Frame Batching (SimulationBatches): For large simulations (many frames and/or large detectors), the full image array can exceed available RAM. When SimulationBatches > 1, the simulation loop stores only lightweight spot records (~20 bytes each: frame, position, intensity) instead of the full image. After the simulation completes, frames are rendered and written in batches:
    1. Pass 1: Each batch of frames is rendered from the spot buffer to find the global maximum pixel intensity.
    2. Pass 2: Each batch is re-rendered, normalized, compressed, and written to the ZIP file.
graph TD
    A["OMP Simulation Loop<br/>(parallel over grains)"] --> B["Compute spot positions<br/>+ SpotMatrixGen.csv"]
    B --> C{"SimulationBatches"}

    C -->|"= 1<br/>(single pass)"| D["Accumulate Gaussians<br/>into float ImageArr<br/>(omp atomic)"]
    D --> E["Scan ImageArr for maxInt"]
    E --> F["Normalize → uint16 → blosc → ZIP"]

    C -->|"> 1<br/>(batched)"| G["Buffer ImageSpots<br/>(~20 bytes each)"]
    G --> H["Pass 1: Render each batch<br/>→ find global maxInt"]
    H --> I["Pass 2: Re-render each batch<br/>→ normalize → compress → ZIP"]

    style D fill:#2d6a4f,color:#fff
    style G fill:#1d3557,color:#fff
Loading
  • Float Precision: The image accumulation array uses float (4 bytes) instead of double (8 bytes), halving memory requirements. Since the final output is uint16, float precision (7 significant digits) is more than adequate.
  • Timing: Wall-clock time is measured using omp_get_wtime() to accurately report parallel speedup. Separate timers report distortion map, simulation, and batch rendering times.

1.5. See Also


2. Near-Field Simulation: simulateNF

simulateNF is the near-field counterpart to ForwardSimulationCompressed. While FF simulation produces continuous-intensity images of diffraction spots, NF simulation produces hit-count images — each pixel records how many diffraction spots illuminated it. This is the data format expected by the NF reconstruction pipeline.

Tip

If you are simulating near-field data, use simulateNF instead of ForwardSimulationCompressed. The NF simulator correctly handles multi-layer detectors, grain geometry, and the NF-specific spot projection model.

2.1. How to Run

./simulateNF <ParameterFile> <InputMicFile> <OutputPrefix> [nCPUs]

Arguments:

  • <ParameterFile>: Path to a text file containing simulation parameters (same format as FF, with NF-specific additions).
  • <InputMicFile>: Path to a .mic file describing the grain microstructure.
  • <OutputPrefix>: Base name for output files. Produces <prefix>, <prefix>.bin, SpotsInfo.bin, SimulatedSpots.csv, and optionally <prefix>.zip.
  • [nCPUs] (Optional): Number of CPU threads for OpenMP parallelization. Default: 1.

Example:

./simulateNF nf_params.txt microstructure.mic NF_simulation_output 16

2.2. Parameter File

The NF parameter file shares many keywords with the FF parameter file, plus some NF-specific ones.

Geometry (per layer):

  • nDistances <N>: Number of detector layers (distances). Must appear before Lsd and BC lines.
  • Lsd <value>: Sample-to-detector distance (microns). Specify one line per layer.
  • BC <y> <z>: Beam center (pixels). Specify one line per layer.

Omega Scanning:

  • OmegaStart <value>: Starting omega angle (degrees).
  • OmegaStep <value>: Omega step per frame (degrees).
  • StartNr <value>: First frame number.
  • EndNr <value>: Last frame number (inclusive). Total frames = EndNr − StartNr + 1.
  • OmegaRange <min> <max>: Omega range to include (degrees). Can appear multiple times.
  • BoxSize <y1> <y2> <z1> <z2>: Bounding box for each omega range. One per OmegaRange.

Material:

  • LatticeParameter <a> <b> <c> <α> <β> <γ>: Lattice constants.
  • SpaceGroup <number>: Space group number for HKL generation.
  • Wavelength <value>: X-ray wavelength (Angstroms).
  • MaxRingRad <value>: Maximum ring radius on the detector (microns).
  • ExcludePoleAngle <value>: Exclude reflections near the pole (degrees).
  • RingsToUse <ring_nr>: Ring numbers to include. Can appear multiple times.

Detector:

  • NrPixelsY <value>: Horizontal detector dimension in pixels.
  • NrPixelsZ <value>: Vertical detector dimension in pixels.
  • NrPixels <value>: (deprecated) Square detector shorthand — sets both NrPixelsY and NrPixelsZ.
  • px <value>: Pixel size (microns).
  • tx <value>, ty <value>, tz <value>: Detector tilt angles (degrees).
  • Wedge <value>: Wedge angle correction (degrees). Default: 0.

Output Control:

  • WriteImage <0_or_1>: If 1 (default), write a Zarr/ZIP file with the simulated NF images. If 0, skip image writing.
  • SaveReducedOutput: If present, skip writing the full raw binary output file.

2.3. Input: .mic File Format

The .mic file describes a 2D grain microstructure. It contains a 4-line header (skipped), followed by one line per grain/voxel:

<header line 1>
<header line 2>
<header line 3>
<header line 4>
<dummy> <dummy> <dummy> <xs> <ys> <edgeLen> <ud> <eul1> <eul2> <eul3> <confidence> <dummy>
...
Column Description
xs, ys Grain center position (microns)
edgeLen Grain edge length (microns)
ud Triangle orientation (+1 or −1)
eul1, eul2, eul3 Euler angles (radians)
confidence Reconstruction confidence

Important

Euler angles in the .mic file must be in radians. The code internally converts them to degrees before computing the orientation matrix.

2.4. Output Files

File Format Description
<prefix> Binary (uint16) Raw simulation array: [nLayers × nrFiles × NrPixelsZ × NrPixelsY]. Each value = spot hit count. Preceded by an 8192-byte header.
<prefix>.bin Binary (uint16) Sparse format: 5 values per illuminated pixel [row, col, frame, layer, count].
SpotsInfo.bin Binary (int) Bit array indicating which pixels were illuminated.
SimulatedSpots.csv CSV Per-spot log: VoxRowNr, DistanceNr, FrameNr, HorPx, VerPx, OmegaRaw, YRaw, ZRaw.
<prefix>.zip Zarr/ZIP (Optional, WriteImage 1) 4D Zarr array [nLayers, nrFiles, NrPixelsZ, NrPixelsY], dtype uint16, blosc/zstd compressed.

2.5. Reading the Zarr/ZIP Output in Python

import zarr
import matplotlib.pyplot as plt
import numpy as np

# Open the Zarr/ZIP file
store = zarr.ZipStore('NF_simulation_output.zip', mode='r')
root = zarr.group(store=store)
data = root['exchange']['data']

print("Data shape:", data.shape)  # (nLayers, nrFiles, NrPixelsZ, NrPixelsY)
print("Data type:", data.dtype)   # uint16

# Display frame 0 of layer 0
fig, ax = plt.subplots()
ax.imshow(data[0, 0, :, :], cmap='hot')
ax.set_title("NF Simulation - Layer 0, Frame 0")
plt.colorbar(ax.images[0], label='Spot Hit Count')
plt.show()

store.close()

2.6. Example NF Parameter File

# ==================================
# Example NF Simulation Parameters
# ==================================

# --- Detector Layers ---
nDistances        2
Lsd               5000.0
Lsd               7000.0
BC                1024.0 1024.0
BC                1024.0 1024.0

# --- Omega Scanning ---
StartNr           1
EndNr             360
OmegaStart        0.0
OmegaStep         1.0
OmegaRange        -180.0 180.0
BoxSize           -1500000 1500000 -1500000 1500000

# --- Material ---
SpaceGroup        225
LatticeParameter  3.6 3.6 3.6 90 90 120
Wavelength        0.172979
MaxRingRad        4000.0
ExcludePoleAngle  10.0
RingsToUse        1
RingsToUse        2
RingsToUse        3

# --- Detector ---
px                1.48
NrPixelsY         2048
NrPixelsZ         2048
tx                0.0
ty                0.0
tz                0.0

# --- Output ---
WriteImage        1

2.7. Technical Implementation Details

Parallelization

  • OpenMP: The main voxel simulation loop is parallelized over grains using OpenMP. The number of threads is controlled by the optional nCPUs argument.
  • Timing: Wall-clock time is measured using omp_get_wtime() to accurately report parallel speedup.

Portable Binary Resolution

  • Both simulateNF and compareNF use the midas_paths.h utility to locate the GetHKLList binary at runtime. The binary is found relative to the running executable's location or via the MIDAS_HOME environment variable, so MIDAS does not need to be installed at a fixed path.

2.8. See Also