Merge branch 'metadata-writer' of github.com:ddasilva/suncet_processi…

…ng_pipeline into metadata-writer

environment.yml

@@ -21,4 +21,4 @@ dependencies:
   - termcolor==2.4.0
   - pytest==7.1.3
   - h5netcdf==1.1.0
-
+  - pykdtree==1.3.13

requirements.txt

@@ -12,3 +12,4 @@ h5netcdf==1.1.0
 bottleneck==1.3.7
 termcolor==2.4.0
 pytest==7.1.3
+pykdtree==1.3.13

suncet_processing_pipeline/make_level1.py

@@ -3,12 +3,14 @@
 """
 import os
 from glob import glob
+import copy
 from pathlib import Path
 import astropy.units as u
 from astropy.io import fits
 import numpy as np
 import sunpy.map
 from scipy import ndimage, interpolate
+from pykdtree.kdtree import KDTree
 from suncet_processing_pipeline import suncet_utility as utilities
 from suncet_processing_pipeline import config_parser
 
@@ -17,6 +19,7 @@ class Level1:
     def __init__(self, config):
         self.config = config
         # self.metadata = self.__load_metadata_from_level0_5()
+        self.cosmic_rays = np.array([]) #TODO: figure out how to save removed cosmic rays.
 
     def __load_metadata_from_level0_5(self):
          pass
@@ -127,7 +130,7 @@ def __badpixel_removal(self, input_map, badpixel_path):
 
     def __cosmic_ray_removal(self, input_map):
         """
-        Interpolates over bad pixels in the input image using bicubic interpolation if the config file has
+        Interpolates over bad pixels in the input image using linear interpolation if the config file has
         cosmic_ray_removal set to True.
 
         Parameters:
@@ -138,7 +141,7 @@ def __cosmic_ray_removal(self, input_map):
         """
         if self.config.cosmic_ray_removal:
 
-            outlier_mask = utilities.detect_outliers(input_map.data, 500)
+            outlier_mask = utilities.detect_outliers(input_map.data, 20)
 
             # Create coordinate grid for good pixels
             y, x = np.indices(input_map.data.shape)
@@ -150,16 +153,33 @@ def __cosmic_ray_removal(self, input_map):
             # Coordinates of bad pixels
             bad_pixel_coords = np.column_stack((y[outlier_mask], x[outlier_mask]))
 
-            # Perform bicubic interpolation over bad pixels
-            interpolated_values = interpolate.griddata(good_pixel_coords, good_pixel_values, bad_pixel_coords, method='cubic')
+            # Build the kdTree
+            tree = KDTree(good_pixel_coords)
+
+            # Query the nearest neighbors for the bad pixel coordinates
+            distances, indices = tree.query(bad_pixel_coords, k=4)  # k=4 nearest neighbors
+
+            # Retrieve the coordinates and values of the 4 nearest neighbors
+            neighbor_coords = good_pixel_coords[indices]
+            neighbor_values = good_pixel_values[indices]
+
+            # Compute bilinear weights based on the inverse of the distances
+            weights = 1 / (distances + 1e-12)  # Avoid division by zero
+            weights /= np.sum(weights, axis=1, keepdims=True)
+
+            # Compute the linear interpolated values as a weighted sums
+            interpolated_values = np.sum(weights * neighbor_values, axis=1)
+
+            # A slow cubic interpolation at the `bad_pixel_coords` points
+            # interpolated_values = interpolate.griddata(good_pixel_coords, good_pixel_values, bad_pixel_coords, method='cubic')
 
             # Replace bad pixel values with interpolated values
-            cosmic_ray_corrected_map = input_map
-            cosmic_ray_corrected_map.data[outlier_mask] = interpolated_values
+            cosmic_ray_corrected_data = copy.deepcopy(input_map.data)
+            cosmic_ray_corrected_data[bad_pixel_coords[:, 0], bad_pixel_coords[:, 1]] = interpolated_values
         else:
-            cosmic_ray_corrected_map = input_map
+            cosmic_ray_corrected_data = input_map.data.copy()
 
-        return cosmic_ray_corrected_map
+        return sunpy.map.Map(cosmic_ray_corrected_data, input_map.meta.copy())
 
     def __coarse_rotate(self, data, telemetry):
         """

suncet_processing_pipeline/suncet_utility.py

@@ -9,13 +9,13 @@
 from scipy.spatial.transform import Rotation as R
 
 
-def detect_outliers(data, threshold=500):
+def detect_outliers(data, threshold=2):
     """
     Detects pixels in the input 2D array that deviate significantly from their 8 nearest neighbors.
 
     Parameters:
         data (numpy.ndarray): Input 2D array (image) containing pixel values.
-        threshold (float): Threshold for deviation from neighbors. Pixels deviating
+        threshold (float): Relative Threshold for deviation from neighbors. Pixels deviating
                           more than this threshold are considered outliers.
 
     Returns:
@@ -28,8 +28,9 @@ def detect_outliers(data, threshold=500):
     # Use a median filter to calculate median value of neighbors
     median = ndimage.median_filter(data, footprint=kernel, mode='reflect')
 
-    # Calculate absolute deviation of each pixel from its 8 neighbors' median
-    deviation = np.abs(data - median)
+    # Calculate relative deviation of each pixel from its 8 neighbors' median
+    # Avoid division by zero by replacing zero median with a small value
+    deviation = np.abs((data - median)/(median + 1e-12))
 
     # Create a boolean mask for pixels deviating more than the threshold
     outlier_mask = deviation > threshold