Support for KPF noise analysis tool

src/eregion/configs/detectors/kpf.yaml

@@ -0,0 +1,106 @@
+description: KPF detectors (2 CCDs with 2 outputs each and 1 CCD with 1 output)
+detector_type: CCD
+detector_output_class: CCDOutput
+
+objects:
+  - name: 'green'
+    class: DetImage
+    header_index: 0
+    filename_format: '*.fits*'
+    properties:
+      x_size: 4188
+      y_size: 4110
+      pixel_size: 0.015  # mm
+    outputs:
+      - id: 'amp1'
+        ext_id: "GREEN_AMP1" # FITS extension ID
+        ext_slice: [!slice [0, 4110], !slice [0, 2094]]
+        data_slice: [!slice [0, 4110], !slice [0, 2094]]
+        serial_prescan: !slice [0, 4] # 4 pix
+        serial_overscan: !slice [2044, 2094] # 50 pix
+        parallel_prescan: !slice [0, 0] # 0 pix
+        parallel_overscan: !slice [4080, 4110] # 30 pix
+        readout_pixel: [0, 0]
+        parallel_axis: 'y'
+        gain: 5.2 # from header
+      - id: 'amp2'
+        ext_id: "GREEN_AMP2"
+        ext_slice: [!slice [0, 4110], !slice [0, 2094]]
+        data_slice: [ !slice [0, 4110], !slice [2094, 4188]]
+        serial_prescan: !slice [4187, 4183] # 4 pix
+        serial_overscan: !slice [2143, 2093] # 50 pix
+        parallel_prescan: !slice [0, 0] # 0 pix
+        parallel_overscan: !slice [4080, 4110] # 30 pix
+        readout_pixel: [0, 4187]
+        parallel_axis: 'y'
+        gain: 5.32  # electrons/ADU
+    focal_plane_position:
+      x_cen: 0.0  # mm
+      y_cen: 0.0  # mm
+      angle: 0.0  # degrees
+      flip_x: false
+      flip_y: false
+
+  - name: 'red'
+    class: DetImage
+    header_index: 0
+    filename_format: '*.fits*'
+    properties:
+      x_size: 4188
+      y_size: 4110
+      pixel_size: 0.015  # mm
+    outputs:
+      - id: 'amp1'
+        ext_id: "RED_AMP1" # FITS extension ID
+        ext_slice: [!slice [0, 4110], !slice [0, 2094]]
+        data_slice: [!slice [0, 4110], !slice [0, 2094]]
+        serial_prescan: !slice [0, 4] # 4 pix
+        serial_overscan: !slice [2044, 2094] # 50 pix
+        parallel_prescan: !slice [0, 0] # 0 pix
+        parallel_overscan: !slice [4080, 4110] # 30 pix
+        readout_pixel: [0, 0]
+        parallel_axis: 'y'
+        gain: 5.0 # from header
+      - id: 'amp2'
+        ext_id: "RED_AMP2"
+        ext_slice: [!slice [0, 4110], !slice [0, 2094]]
+        data_slice: [ !slice [0, 4110], !slice [2094, 4188]]
+        serial_prescan: !slice [4187, 4183] # 4 pix
+        serial_overscan: !slice [2143, 2093] # 50 pix
+        parallel_prescan: !slice [0, 0] # 0 pix
+        parallel_overscan: !slice [4080, 4110] # 30 pix
+        readout_pixel: [0, 4187]
+        parallel_axis: 'y'
+        gain: 5.3  # electrons/ADU
+    focal_plane_position:
+      x_cen: 0.0  # mm
+      y_cen: 0.0  # mm
+      angle: 0.0  # degrees
+      flip_x: false
+      flip_y: false
+
+  - name: 'cahk'
+    class: DetImage
+    header_index: 0
+    filename_format: '*.fits*'
+    properties:
+      x_size: 1024
+      y_size: 255
+      pixel_size: 0.026  # mm
+    outputs:
+      - id: 'amp1'
+        ext_id: "CA_HK" # FITS extension ID
+        ext_slice: [!slice [0, 255], !slice [0, 1024]]
+        data_slice: [!slice [0, 255], !slice [0, 1024]]
+        serial_prescan: !slice [0, 0]
+        serial_overscan: !slice [1023, 1023]
+        parallel_prescan: !slice [0, 0]
+        parallel_overscan: !slice [254, 254]
+        readout_pixel: [0, 0]
+        parallel_axis: 'y'
+    focal_plane_position:
+      x_cen: 0.0  # mm
+      y_cen: 0.0  # mm
+      angle: 0.0  # degrees
+      flip_x: false
+      flip_y: false

src/eregion/core/welch2d.py

@@ -1,6 +1,7 @@
 import numpy as np
 from typing import Callable, Any
-from scipy.signal.windows import boxcar
+from scipy.signal.windows import boxcar, get_window
+from scipy.signal import welch
 
 def split_trim_array_squares(array: np.ndarray, target_size: int) -> list[np.ndarray]:
     """
@@ -75,3 +76,50 @@ def welch2d(array: np.ndarray,
     errpsds = stdev_func(outpsds, axis=0) / np.sqrt(len(outpsds))
 
     return avpsds, errpsds
+
+
+def spatial_to_temporal_freq(fx, fy, hz_x, hz_y):
+    """
+    Convert spatial frequency coordinates from PSD to temporal frequency.
+
+    :param fx: Spatial frequency along the x-axis (in cycles per pixel).
+    :param fy: Spatial frequency along the y-axis (in cycles per pixel).
+    :param hz_x: Temporal frequency corresponding to the x-axis (in Hz), e.g. pixel readout rate.
+    :param hz_y: Temporal frequency corresponding to the y-axis (in Hz),  e.g. 1/(line time).
+    :return: Apparent temporal frequency (Hx)
+    """
+    return np.abs(fx * hz_x + fy * hz_y)
+
+
+def raveled_welch(array: np.ndarray,
+                  target_shape: int,
+                  window_func: str | Callable[..., Any] = 'boxcar',
+                  **kwargs) -> tuple[np.ndarray, ...]:
+    """
+    Ravel a 2D array into a 1D array and compute the PSD using scipy.signal.welch.
+    :param array: Input 2D array to analyze (typically image data).
+    :param target_shape: Segment length for Welch's method (in pixels), passes to `nperseg` kwarg.
+    :param window_func: string or callable, window function to use for Welch's method. If a string, it should be a
+                        valid window name from `scipy.signal.get_window`. If a callable, it should accept the segment
+                        length and return a 1D window array of that length. Default is 'boxcar'.
+    :param kwargs: Additional kwargs to pass to `scipy.signal.welch`, such as `fs` (sampling frequency), `nfft`, etc.
+    :return: Tuple `(f, Pxx)` where `f` is the array of sample frequencies and `Pxx` is the power spectral density of
+            the raveled array.
+    """
+    # Ravel the 2D array into a 1D array
+    raveled_array = array.ravel()
+
+    # Determine the window function
+    if isinstance(window_func, str):
+        window = get_window(window_func, target_shape)
+    elif callable(window_func):
+        window = window_func(target_shape)
+    else:
+        raise ValueError("window_func must be a string or a callable function.")
+
+    # Compute the PSD using scipy.signal.welch
+    kwargs['noverlap'] = kwargs.get('noverlap', 0)  # set to 0 if not provided
+    kwargs['return_onesided'] = kwargs.get('return_onesided', True)  # set to True if not provided
+    f, Pxx = welch(raveled_array, window=window, nperseg=target_shape, **kwargs)
+
+    return f, Pxx

src/eregion/tasks/custom.py

@@ -50,10 +50,34 @@ def process_raw_numbers(data: np.ndarray, offs) -> np.ndarray:
 
 def load_image_fits_DEIMOS(filename: str) -> tuple[list[Any], list[fits.Header]]:
     input_data_array, input_headers = load_image_fits(filename)
-    for i, hdr in enumerate(input_headers):
+    for i, data in enumerate(input_data_array):
+        if data is None:
+            continue
+        hdr = input_headers[i]
         if "TAPOFFS" in hdr:
             input_data_array[i] = process_raw_numbers(input_data_array[i], hdr["TAPOFFS"])
             logger.debug("Processed raw numbers for DEIMOS HDU %d", i)
         else:
             logger.debug("TAPOFFS not found in header for file %s HDU %d", filename, i)
     return input_data_array, input_headers
+
+def load_image_fits_KPF(filename: str) -> tuple[list[Any], list[fits.Header]]:
+    input_data_array, input_headers = load_image_fits(filename)
+    for i, data in enumerate(input_data_array):
+        if data is None:
+            continue
+        hdr = input_headers[i]
+        # Only GREEN/RED image extensions have 32-bit readouts; Ca H&K (16-bit) and tables are left as-is
+        if hdr.get("XTENSION") == "IMAGE" and abs(hdr["BITPIX"]) == 32:
+            offset_key = [key for key in hdr.keys() if 'OFFSET' in key]
+            if len(offset_key) == 0:
+                logger.debug(f"No OFFSET keys found in header for file {filename} HDU {i}")
+                offset = 0.0
+            else:
+                logger.debug(f"Found OFFSET key {offset_key[0]} in header for file {filename} HDU {i}")
+                offset = hdr[offset_key[0]]
+            # GREEN/RED amps store the 32-bit readout as a value (float32 on disk); ADU = value / 2^16.
+            # Convert by value, not view, so byte order and on-disk dtype don't matter.
+            input_data_array[i] = process_raw_numbers(data.astype(np.uint32), offset)
+            logger.debug("Processed raw numbers for KPF HDU %d", i)
+    return input_data_array, input_headers

src/eregion/tasks/imagegen.py

@@ -99,9 +99,11 @@ def __init__(
         watch_mode=False,
         poll_interval=0,
         max_batch_size=10,
+        settle_time: float = 0.0,
         **kwargs,
     ):
         super().__init__(name, watch_mode, poll_interval, max_batch_size, **kwargs)
+        self.settle_time = settle_time
 
         # Load detector configuration
         self.det_config = DetectorConfig(config_input=detector_config)
@@ -190,34 +192,38 @@ def _discover_once() -> list[str]:
                     else:
                         self.logger.info(f"Item {ipath} is a regular path string")
                         list_to_process.append(ipath)
-                return parse_list_of_files(list_to_process)
+                found_files = parse_list_of_files(list_to_process)
 
             elif isinstance(input_path, str):
                 if "*" in input_path:
                     self.logger.info(f"Item {input_path} is a glob pattern")
-                    return parse_list_of_files(glob2.glob(input_path, recursive=True))
+                    found_files = parse_list_of_files(glob2.glob(input_path, recursive=True))
                 else:
                     self.logger.info(f"Item {input_path} is a regular path string")
-                    return parse_list_of_files([input_path])
+                    found_files = parse_list_of_files([input_path])
 
             else:
                 self.logger.error(f"Invalid input")
                 raise ValueError(
                     "input_path must be a string or list of strings representing file paths, directories, or glob patterns."
                 )
+            if self.watch_mode:
+                now = time.time()
+                found_files = [f for f in found_files if now-os.path.getmtime(f) >= self.settle_time]
+            return found_files
 
         if not self.watch_mode:
             batch = _discover_once()
             yield batch
             return
         # Watch mode
-        while True:
+        while not self._stop_event.is_set():
             batch = _discover_once()
             new_files = [f for f in batch if f not in self._seen_files]
             if new_files:
-                self._seen_files.update(new_files)
                 yield new_files
-            time.sleep(self.poll_interval)
+            if self._stop_event.wait(self.poll_interval):
+                return
 
     def from_arrays(
         self, input_arrays: Iterable[np.ndarray]
@@ -416,6 +422,7 @@ def lazy_run(
         :return: {"images": list of DetImage}
             Generator yielding lists of DetImage objects, stored under the key 'images' in the yielded dict.
         """
+        self._stop_event.clear()
         if self.watch_mode:
             self.logger.info(
                 "Running in watch mode, will monitor input source for new data, setting require_data=False"
@@ -462,6 +469,8 @@ def alloftype(arr, tp):
 
             ## split batch into smaller batches of max_batch_size
             for batch in batched(input_batch, self.max_batch_size):
+                if self._stop_event.is_set():
+                    return
                 images_batch = []
                 for inp in batch:
                     if _file_input:
@@ -471,6 +480,8 @@ def alloftype(arr, tp):
                         inp, data_on_demand=data_on_demand
                     )
                     images_batch.extend(images)
+                if _file_input and self.watch_mode:
+                    self._seen_files.update(batch)
                 yield self.task_result(data=images_batch)
 
     @wraps(lazy_run)