Auto filter

.gitignore

@@ -6,6 +6,10 @@ build
 # other typical case file names
 *.vti
 *morphology.csv
+*.png
 
 # vscode settings cache
 .vscode
+
+# macOS Finder metadata
+.DS_Store

examples/synthetic_api_example/run_synthetic_rve.py

@@ -13,102 +13,118 @@
 from raptor.api import (
     compute_morphology,
     compute_porosity,
+    compute_spectral_components,
     create_grid,
     create_melt_pool,
     create_path_vectors,
+    visualize,
     write_morphology,
     write_vtk,
 )
-from raptor.utilities import MeltPoolFilter
+from raptor.utilities import (
+    MeltPoolFilter,
+    plot_melt_pool_signal,
+    reconstruct_spectral_signal,
+)
 
 
-# =============================================================================
-# User configuration
-# =============================================================================
-# These are the main parameters to vary in a standard AM optimization study.
-# All lengths are in meters, speed is in m/s, power is in W, and angles are in
-# degrees.
-
-# Process parameters
-# Accepted by the API, but does not affect the result (see note below).
+# 0. User defined process parameters
 LASER_POWER = 370.0
 SCAN_SPEED = 1.7
 HATCH_SPACING = 140.0e-6
 LAYER_HEIGHT = 30.0e-6
 SCAN_ROTATION = 67.0
 
-# Mean melt-pool geometry
 MELT_POOL_WIDTH = 148.0e-6
 MELT_POOL_DEPTH = 118.0e-6
-MELT_POOL_HEIGHT = 59.2e-6
-
-# Stochastic melt-pool model
+MELT_POOL_HEIGHT = 30.0e-6
 MELT_POOL_WIDTH_STD_DEV = 18.0e-6
 MELT_POOL_LENGTH = 300.0e-6
-N_SPECTRAL_MODES = 50
-RANDOM_SEED = 42
 
-# Lamé-curve exponents: 1 = parabola, 2 = ellipse.
 HEIGHT_SHAPE_FACTOR = 1.0
 DEPTH_SHAPE_FACTOR = 1.0
 
-# RVE and output controls
 RVE_MIN_POINT = np.array([0.0, 0.0, 0.0])
 RVE_MAX_POINT = np.array([5.0e-4, 5.0e-4, 5.0e-4])
 VOXEL_RESOLUTION = 5.0e-6
+
 VTK_OUTPUT = "rve.vti"
 MORPHOLOGY_OUTPUT = "rve_morphology.csv"
+WIDTH_DATA_PLOT = "melt_pool_width_timeseries.png"
+ENABLE_VISUALIZATION = False
 
 
 def build_melt_pool():
-    """Create a stochastic melt pool from the top-level configuration."""
-    np.random.seed(RANDOM_SEED)
-
-    # Internal signal-generation settings. These control the numerical
-    # representation of the stochastic melt-pool history.
-    sampling_frequency = 250000.0
-    time_series_duration = 0.08
+    # Create melt pools from convolution filter
 
-    melt_pool_filter = MeltPoolFilter(
+    # Instantiate object
+    mp_filter = MeltPoolFilter(
         MELT_POOL_WIDTH,
         MELT_POOL_WIDTH_STD_DEV,
         SCAN_SPEED,
-        [sampling_frequency, time_series_duration],
+        VOXEL_RESOLUTION,
+        random_seed=42,
+    )
+
+    # Define physical scales
+    mp_filter.add_effect("melt_pool", [MELT_POOL_LENGTH, None, 1])
+
+    # Generate stochastic melt pool
+    mp_filter.initialize()
+    width_data = mp_filter.generate_fluctuations(1.0, mp_filter.n_points, mp_filter.t)
+    width_spectral = compute_spectral_components(width_data, tolerance=VOXEL_RESOLUTION)
+    reconstructed_width = reconstruct_spectral_signal(width_data[:, 0], width_spectral)
+    reconstruction_rmse = np.sqrt(
+        np.mean((width_data[:, 1] - reconstructed_width) ** 2)
+    )
+    print(
+        "Melt-pool signal: "
+        f"duration={mp_filter.duration:.6e} s, "
+        f"samples={mp_filter.n_points}, "
+        f"n_modes={width_spectral.shape[0]}, "
+        f"reconstruction_rmse={reconstruction_rmse:.6e} m"
+    )
+    plot_melt_pool_signal(
+        width_data,
+        width_spectral,
+        MELT_POOL_WIDTH,
+        MELT_POOL_WIDTH_STD_DEV,
+        WIDTH_DATA_PLOT,
+        value_label="Melt-pool width (µm)",
     )
 
-    # MeltPoolFilter converts this length scale to a temporal frequency
-    # using SCAN_SPEED / MELT_POOL_LENGTH.
-    melt_pool_filter.add_effect("melt_pool", [MELT_POOL_LENGTH, None, 1.0])
-    melt_pool_filter.initialize()
-    width_data = melt_pool_filter.generate_fluctuations(noise_scale=1.0)
+    # scale melt pool data by constant factor
+    depth_scale = MELT_POOL_DEPTH / MELT_POOL_WIDTH
+    height_scale = MELT_POOL_HEIGHT / MELT_POOL_WIDTH
 
-    # Depth and height reuse the width history, scaled to their requested means.
+    # assign shape to melt pool and cap (1 = parabola, 2 = ellipse)
     melt_pool_dict = {
-        "width": (width_data, N_SPECTRAL_MODES, 1.0, 1.0),
+        "width": (width_data, None, 1.0, 2.0),
         "depth": (
             width_data,
-            N_SPECTRAL_MODES,
-            MELT_POOL_DEPTH / MELT_POOL_WIDTH,
+            None,
+            depth_scale,
             DEPTH_SHAPE_FACTOR,
         ),
         "height": (
             width_data,
-            N_SPECTRAL_MODES,
-            MELT_POOL_HEIGHT / MELT_POOL_WIDTH,
+            None,
+            height_scale,
             HEIGHT_SHAPE_FACTOR,
         ),
     }
 
-    return create_melt_pool(melt_pool_dict, enable_random_phases=True)
+    return create_melt_pool(
+        melt_pool_dict, enable_random_phases=True, tolerance=VOXEL_RESOLUTION
+    )
 
 
 def main():
+    # 1. Create voxel grid for the representative volume element (RVE)
     bound_box = np.array([RVE_MIN_POINT, RVE_MAX_POINT])
     grid = create_grid(VOXEL_RESOLUTION, bound_box=bound_box)
 
-    # LASER_POWER is accepted by the current scan-path API but is not used by
-    # the porosity calculation. An optimization must supply its relationship to
-    # melt-pool geometry through the parameters.
+    # 2. Create path vectors through the representative volume element (RVE)
     path_vectors = create_path_vectors(
         bound_box,
         LASER_POWER,
@@ -117,21 +133,30 @@ def main():
         LAYER_HEIGHT,
         SCAN_ROTATION,
         scan_extension=max(RVE_MAX_POINT - RVE_MIN_POINT),
-        extra_layers=0,
+        extra_layers=10,
     )
 
+    # 3. Create melt pools from convolution filter
     melt_pool = build_melt_pool()
+
+    # 4. Compute porosity using conic section / superellipse curves for melt pool mask
     porosity = compute_porosity(grid, path_vectors, melt_pool, jit_warmup=True)
 
+    # 5. Write porosity field to .VTI
     write_vtk(grid.origin, grid.resolution, porosity, VTK_OUTPUT)
 
+    # 6. Compute morphology
     morphology = compute_morphology(
         porosity,
         VOXEL_RESOLUTION,
         ["area", "equivalent_diameter_area"],
     )
     write_morphology(morphology, MORPHOLOGY_OUTPUT)
 
+    # 7. Optionally visualize using PyVista (requires a graphical display)
+    if ENABLE_VISUALIZATION:
+        visualize(VTK_OUTPUT)
+
 
 if __name__ == "__main__":
     main()

pyproject.toml

@@ -23,6 +23,8 @@ classifiers = [
 dependencies = [
     "numpy",
     "numba",
+    "scipy",
+    "matplotlib",
     "PyYAML",
     "vtk",
     "scikit-image",

src/raptor/api.py

@@ -9,7 +9,7 @@
 # https://github.com/ORNL-MDF/Raptor/LICENSE
 # =============================================================================
 import time
-from typing import List, Tuple, Optional, Dict, Any
+from typing import Any, Dict, List, Optional
 import numpy as np
 import pandas as pd
 import vtk
@@ -65,72 +65,139 @@ def create_path_vectors(
     return scan_path_builder.process_vectors()
 
 
-def compute_spectral_components(melt_pool_data: np.ndarray, n_modes: int) -> np.ndarray:
+def compute_required_modes(data: np.ndarray, reconstruction_rmse: float) -> int:
+    """Return the minimum number of real Fourier modes for an RMSE target."""
+    return compute_spectral_components(data, tolerance=reconstruction_rmse).shape[0]
 
-    dt = melt_pool_data[1, 0] - melt_pool_data[0, 0]
-    mode0 = melt_pool_data[:, 1].mean()
-    fft_resolution = np.fft.fft(melt_pool_data[:, 1])
-    F = np.zeros_like(fft_resolution)
-    n_fft = len(fft_resolution)
-    if n_modes == 1:
-        spectral_array = np.array([[mode0, 0, 0]])
-    else:
-        for i in range(1, n_modes):
-            F[i] = fft_resolution[i]
-            F[n_fft - i] = fft_resolution[n_fft - i]
 
-        frequencies = np.float64(1 / (dt * n_fft)) * np.arange(
-            n_modes, dtype=np.float64
+def compute_spectral_components(
+    melt_pool_data: np.ndarray,
+    n_modes: Optional[int] = None,
+    tolerance: Optional[float] = None,
+) -> np.ndarray:
+    """Convert a uniformly sampled real signal to a sparse cosine expansion.
+
+    ``n_modes`` includes the mean (DC) mode and retains the corresponding
+    low-frequency prefix.  When ``tolerance`` is used, the smallest set of
+    Fourier bins whose discarded energy satisfies the requested reconstruction
+    RMSE is retained.
+    """
+    data = np.asarray(melt_pool_data, dtype=np.float64)
+    if data.ndim != 2 or data.shape[1] != 2 or data.shape[0] < 2:
+        raise ValueError("melt_pool_data must have shape (n, 2) with n >= 2.")
+    if not np.isfinite(data).all():
+        raise ValueError("melt_pool_data must contain only finite values.")
+    if (n_modes is None) == (tolerance is None):
+        raise ValueError("Provide exactly one of n_modes or tolerance.")
+
+    time_values = data[:, 0]
+    signal = data[:, 1]
+    time_steps = np.diff(time_values)
+    dt = time_steps[0]
+    if dt <= 0.0 or not np.allclose(time_steps, dt, rtol=1.0e-7, atol=0.0):
+        raise ValueError("Time samples must be strictly increasing and uniform.")
+
+    n_samples = signal.size
+    fft_values = np.fft.rfft(signal)
+    frequencies = np.fft.rfftfreq(n_samples, d=dt)
+    candidate_bins = np.arange(1, fft_values.size)
+
+    # Parseval energy represented by each positive-frequency cosine.  Interior
+    # rFFT bins represent a conjugate pair; the Nyquist bin does not.
+    energy_weights = np.full(fft_values.size, 2.0)
+    energy_weights[0] = 1.0
+    if n_samples % 2 == 0:
+        energy_weights[-1] = 1.0
+    if tolerance is not None:
+        if not np.isfinite(tolerance) or tolerance < 0.0:
+            raise ValueError("tolerance must be a finite, non-negative value.")
+        candidate_energy = (
+            energy_weights[candidate_bins] * np.abs(fft_values[candidate_bins]) ** 2
         )
-        phases = np.float64(np.angle(F[:n_modes]))
-        amplitudes = np.float64(np.abs(F[:n_modes]) / n_fft)
-        spectral_array = np.vstack(
-            [
-                np.array([mode0, 0, 0]),
-                np.vstack([amplitudes[1:], frequencies[1:], phases[1:]]).transpose(),
-            ]
+        energy_indices = np.argsort(candidate_energy)[::-1]
+        energy_order = candidate_bins[energy_indices]
+        allowed_energy = (tolerance * n_samples) ** 2
+        ordered_energy = candidate_energy[energy_indices]
+        required_energy = ordered_energy.sum() - allowed_energy
+        retained_count = (
+            0
+            if required_energy <= 0.0
+            else np.searchsorted(
+                np.cumsum(ordered_energy), required_energy, side="left"
+            )
+            + 1
         )
-    return np.float64(spectral_array)
+        selected_bins = energy_order[:retained_count]
+    else:
+        if not isinstance(n_modes, (int, np.integer)) or n_modes < 1:
+            raise ValueError("n_modes must be a positive integer.")
+        if n_modes > fft_values.size:
+            raise ValueError(
+                f"n_modes cannot exceed {fft_values.size} for this time series."
+            )
+        selected_bins = candidate_bins[: n_modes - 1]
+
+    # Frequency order is convenient for evaluation and deterministic output.
+    selected_bins = np.sort(selected_bins)
+    amplitudes = energy_weights[selected_bins] * np.abs(fft_values[selected_bins])
+    amplitudes /= n_samples
+    phases = np.angle(fft_values[selected_bins])
+    # FFT phases are relative to sample zero; account for an absolute time axis.
+    phases -= 2.0 * np.pi * frequencies[selected_bins] * time_values[0]
+
+    spectral_array = np.empty((selected_bins.size + 1, 3), dtype=np.float64)
+    spectral_array[0] = (signal.mean(), 0.0, 0.0)
+    spectral_array[1:, 0] = amplitudes
+    spectral_array[1:, 1] = frequencies[selected_bins]
+    spectral_array[1:, 2] = phases
+    return spectral_array
 
 
 def create_melt_pool(
-    melt_pool_dict: Dict[str, Any], enable_random_phases: bool
+    melt_pool_dict: Dict[str, Any],
+    enable_random_phases: bool,
+    tolerance: Optional[float] = None,
 ) -> MeltPool:
 
-    processed_components: Dict[str, Tuple[np.ndarray, float]] = {}
-    max_modes = 0
-
-    # 1. Determine the maximum number of modes required.
-    for _, nmodes, _, _ in melt_pool_dict.values():
-        max_modes = max(max_modes, nmodes)
+    processed_components: Dict[str, np.ndarray] = {}
 
-    # 2. Process each component into its spectral format
+    # Process every component before padding so results do not depend on
+    # dictionary insertion order.
     for key, (data, n_modes, scale, shape_factor) in melt_pool_dict.items():
         # Option A: Input data is a raw time-series [time, value]
         if data.shape[1] == 2:
-            spectral_array = compute_spectral_components(data, n_modes)
+            component_tolerance = tolerance
+            if tolerance is not None and scale != 0.0:
+                component_tolerance = tolerance / abs(scale)
+            spectral_array = compute_spectral_components(
+                data,
+                n_modes=n_modes,
+                tolerance=component_tolerance if n_modes is None else None,
+            )
             spectral_array[:, 0] *= scale
 
         # Option B: Input data is a spectral array [amplitude, frequency, phase]
         elif data.shape[1] == 3:
             spectral_array = data.copy()
+            spectral_array[:, 0] *= scale
 
         else:
             raise ValueError(
                 f"Unsupported data shape: {data.shape}.  "
                 f"Must be [time, value] or [amplitude, frequency, phase]"
             )
 
-        # Pad the array with zeros if it has fewer modes than the max.
+        processed_components[key] = np.asarray(spectral_array, dtype=np.float64)
+
+    max_modes = max(array.shape[0] for array in processed_components.values())
+    for key, spectral_array in processed_components.items():
         current_modes = spectral_array.shape[0]
-        if current_modes < max_modes:
+        if current_modes != max_modes:
             pad_array = np.zeros(
                 shape=(max_modes - current_modes, spectral_array.shape[1]),
                 dtype=np.float64,
             )
-            spectral_array = np.vstack([spectral_array, pad_array])
-
-        processed_components[key] = spectral_array
+            processed_components[key] = np.vstack([spectral_array, pad_array])
 
     # 3. Create the MeltPool object
     width_oscillations = processed_components["width"]