Feature: Haralick Texture Analysis

docs/source/conf.py

@@ -35,6 +35,7 @@
 author = "Benjamin Smith"
 copyright = f"{datetime.now():%Y}, {author}"
 
+
 # Attempt to get version from package; fall back gracefully.
 def _read_version():
     try:
@@ -44,6 +45,7 @@ def _read_version():
         # On RTD before build of extension module this can fail.
         return "0.0.0"
 
+
 version = _read_version()
 release = version
 
@@ -156,6 +158,7 @@ def _read_version():
 .. |NumPy| replace:: **NumPy**
 """
 
+
 # ---------------------------------------------------------------------------
 # Custom hook: ensure README & QUICKSTART included without duplication
 # ---------------------------------------------------------------------------
@@ -174,11 +177,15 @@ def _link_external_markdown():
             try:
                 text = src_path.read_text(encoding="utf-8")
                 # Avoid writing if unchanged
-                if not dest_path.exists() or dest_path.read_text(encoding="utf-8") != text:
+                if (
+                    not dest_path.exists()
+                    or dest_path.read_text(encoding="utf-8") != text
+                ):
                     dest_path.write_text(text, encoding="utf-8")
             except Exception as e:
                 print(f"[conf.py] Warning: could not process {fname}: {e}")
 
+
 _link_external_markdown()
 
 # ---------------------------------------------------------------------------
@@ -213,6 +220,7 @@ def _link_external_markdown():
 distributed chunking with Rust's local parallel execution.
 """
 
+
 def setup(app):
     # Make the note available as a config value
     app.add_config_value("rust_arch_note", RUST_ARCH_NOTE, "env")

examples/basic_usage.py

@@ -32,7 +32,7 @@
 print("Example 2: NDVI with 2D arrays (100x100 image)")
 print("-" * 40)
 nir_2d = np.random.rand(100, 100) * 0.8 + 0.2  # NIR values between 0.2 and 1.0
-red_2d = np.random.rand(100, 100) * 0.4        # Red values between 0.0 and 0.4
+red_2d = np.random.rand(100, 100) * 0.4  # Red values between 0.0 and 0.4
 ndvi_2d_result = ndvi(nir_2d, red_2d)
 print(f"NIR shape:  {nir_2d.shape}")
 print(f"Red shape:  {red_2d.shape}")
@@ -56,11 +56,13 @@
 # Example 4: NDWI with 2D arrays
 print("Example 4: NDWI with 2D arrays (50x50)")
 print("-" * 40)
-green_2d = np.random.rand(50, 50) * 0.5 + 0.1   # Green between 0.1 and 0.6
-nir_2d = np.random.rand(50, 50) * 0.4 + 0.1     # NIR between 0.1 and 0.5
+green_2d = np.random.rand(50, 50) * 0.5 + 0.1  # Green between 0.1 and 0.6
+nir_2d = np.random.rand(50, 50) * 0.4 + 0.1  # NIR between 0.1 and 0.5
 ndwi_2d = ndwi(green_2d, nir_2d)
 print(f"NDWI shape: {ndwi_2d.shape}")
-print(f"NDWI stats -> min: {ndwi_2d.min():.4f} max: {ndwi_2d.max():.4f} mean: {ndwi_2d.mean():.4f}")
+print(
+    f"NDWI stats -> min: {ndwi_2d.min():.4f} max: {ndwi_2d.max():.4f} mean: {ndwi_2d.mean():.4f}"
+)
 print()
 
 # Example 5: Enhanced Vegetation Index (EVI) 1D
@@ -79,12 +81,14 @@
 # Example 6: Enhanced Vegetation Index (EVI) 2D
 print("Example 6: Enhanced Vegetation Index (EVI) 2D (60x60)")
 print("-" * 40)
-nir_evi_2d = np.random.rand(60, 60) * 0.6 + 0.2   # 0.2 - 0.8
-red_evi_2d = np.random.rand(60, 60) * 0.4 + 0.1   # 0.1 - 0.5
-blue_evi_2d = np.random.rand(60, 60) * 0.2 + 0.05 # 0.05 - 0.25
+nir_evi_2d = np.random.rand(60, 60) * 0.6 + 0.2  # 0.2 - 0.8
+red_evi_2d = np.random.rand(60, 60) * 0.4 + 0.1  # 0.1 - 0.5
+blue_evi_2d = np.random.rand(60, 60) * 0.2 + 0.05  # 0.05 - 0.25
 evi_2d = evi(nir_evi_2d, red_evi_2d, blue_evi_2d)
 print(f"EVI shape: {evi_2d.shape}")
-print(f"EVI stats -> min: {evi_2d.min():.4f} max: {evi_2d.max():.4f} mean: {evi_2d.mean():.4f}")
+print(
+    f"EVI stats -> min: {evi_2d.min():.4f} max: {evi_2d.max():.4f} mean: {evi_2d.mean():.4f}"
+)
 print()
 
 # Example 7: Generic normalized difference
@@ -102,6 +106,7 @@
 print("Example 8: Performance (1000x1000 NDVI)")
 print("-" * 40)
 import time
+
 size = 1000
 nir_large = np.random.rand(size, size)
 red_large = np.random.rand(size, size)
@@ -114,7 +119,9 @@
 ndvi_numpy = (nir_large - red_large) / (nir_large + red_large)
 t_numpy = time.time() - t0
 
-print(f"Rust:  {t_rust*1000:.2f} ms  NumPy: {t_numpy*1000:.2f} ms  Speedup: {t_numpy/t_rust:.2f}x  Match: {np.allclose(ndvi_rust, ndvi_numpy, rtol=1e-10)}")
+print(
+    f"Rust:  {t_rust * 1000:.2f} ms  NumPy: {t_numpy * 1000:.2f} ms  Speedup: {t_numpy / t_rust:.2f}x  Match: {np.allclose(ndvi_rust, ndvi_numpy, rtol=1e-10)}"
+)
 print()
 
 # Example 9: Temporal statistics & median composite
@@ -127,7 +134,9 @@
 composite_img = composite(ts, method="median")
 
 print(f"mean shape:    {mean_img.shape}, std shape: {std_img.shape}")
-print(f"median shape:  {median_img.shape}, composite (median) identical: {np.allclose(median_img, composite_img)}")
+print(
+    f"median shape:  {median_img.shape}, composite (median) identical: {np.allclose(median_img, composite_img)}"
+)
 print(f"mean range:    [{mean_img.min():.4f}, {mean_img.max():.4f}]")
 print(f"std range:     [{std_img.min():.4f}, {std_img.max():.4f}]")
 print(f"median range:  [{median_img.min():.4f}, {median_img.max():.4f}]")

examples/map_blocks.py

@@ -29,12 +29,15 @@
     print("  pip install eo-processor[dask] xarray dask")
     raise
 
+
 # Helper to report timing and basic stats
 def report(name: str, result_array):
     """Print name, time and some basic stats of the computed result."""
     if isinstance(result_array, xr.DataArray):
         arr = result_array
-        data = arr.values if not isinstance(arr.data, da.Array) else arr.compute().values
+        data = (
+            arr.values if not isinstance(arr.data, da.Array) else arr.compute().values
+        )
     elif isinstance(result_array, da.Array):
         data = result_array.compute()
     else:
@@ -85,8 +88,12 @@ def example_map_blocks_vs_apply_ufunc():
     nir_dask = da.from_array(nir_np)
     red_dask = da.from_array(red_np)
 
-    nir_xr = xr.DataArray(nir_dask, dims=("y", "x"), coords={"y": np.arange(size), "x": np.arange(size)})
-    red_xr = xr.DataArray(red_dask, dims=("y", "x"), coords={"y": np.arange(size), "x": np.arange(size)})
+    nir_xr = xr.DataArray(
+        nir_dask, dims=("y", "x"), coords={"y": np.arange(size), "x": np.arange(size)}
+    )
+    red_xr = xr.DataArray(
+        red_dask, dims=("y", "x"), coords={"y": np.arange(size), "x": np.arange(size)}
+    )
 
     # Use apply_ufunc which can run the underlying Rust UDF in parallel.
     print("Timing: xarray.apply_ufunc (dask='parallelized') ...")
@@ -100,7 +107,7 @@ def example_map_blocks_vs_apply_ufunc():
         dask="parallelized",
         vectorize=False,
         output_dtypes=[float],
-        dask_gufunc_kwargs={'allow_rechunk': True}
+        dask_gufunc_kwargs={"allow_rechunk": True},
     )
     # Trigger computation and measure
     ndvi_xr_ufunc_computed = ndvi_xr_ufunc.compute()
@@ -112,6 +119,7 @@ def example_map_blocks_vs_apply_ufunc():
     # xarray.map_blocks applies a function block-by-block. The function should accept
     # numpy arrays (blocks) and return a numpy array or xarray DataArray for that block.
     print("Timing: xarray.map_blocks ...")
+
     # define block function that uses the Rust ndvi on numpy blocks
     def block_ndvi(darr_chunk: xr.DataArray):
         # ds will have 'nir' and 'red' DataArrays for the current block
@@ -132,7 +140,6 @@ def block_ndvi(darr_chunk: xr.DataArray):
 
     print("Preparing stacked xarray for map_blocks...")
 
-
     # Build xarray.DataArray inputs (already defined as nir_xr/red_xr)
     start = time.time()
     stacked_xr = xr.concat(
@@ -158,6 +165,7 @@ def block_ndvi(darr_chunk: xr.DataArray):
 
     # 4) dask.array.map_blocks on the underlying dask arrays
     print("Timing: dask.array.map_blocks ...")
+
     def dask_block_ndvi(nir_block, red_block):
         # this will be called with numpy arrays per block
         res_arr = ndvi(nir_block, red_block)
@@ -179,9 +187,20 @@ def dask_block_ndvi(nir_block, red_block):
 
     # Quick consistency checks (allow tiny floating differences)
     print("Sanity checks (all_close to NumPy baseline):")
-    print("apply_ufunc ~ baseline:", np.allclose(ndvi_xr_ufunc_computed, ndvi_numpy, equal_nan=True, atol=1e-10))
-    print("xarray.map_blocks ~ baseline:", np.allclose(ndvi_xr_mapblocks_computed, ndvi_numpy, equal_nan=True, atol=1e-10))
-    print("dask.map_blocks ~ baseline:", np.allclose(ndvi_dask_mapblocks_computed, ndvi_numpy, equal_nan=True, atol=1e-10))
+    print(
+        "apply_ufunc ~ baseline:",
+        np.allclose(ndvi_xr_ufunc_computed, ndvi_numpy, equal_nan=True, atol=1e-10),
+    )
+    print(
+        "xarray.map_blocks ~ baseline:",
+        np.allclose(ndvi_xr_mapblocks_computed, ndvi_numpy, equal_nan=True, atol=1e-10),
+    )
+    print(
+        "dask.map_blocks ~ baseline:",
+        np.allclose(
+            ndvi_dask_mapblocks_computed, ndvi_numpy, equal_nan=True, atol=1e-10
+        ),
+    )
     print()
 
     # Summary timings

examples/morphology_example.py

@@ -2,10 +2,12 @@
 import eo_processor
 from eo_processor import binary_dilation, binary_erosion, binary_opening, binary_closing
 
+
 def print_grid(arr):
     for row in arr:
         print(" ".join(str(x) for x in row))
 
+
 def main():
     print("EO Processor - Morphology Example")
     print("=================================")
@@ -15,8 +17,8 @@ def main():
     # A 3x3 block in the center with some noise
     input_arr = np.zeros((7, 7), dtype=np.uint8)
     input_arr[2:5, 2:5] = 1
-    input_arr[0, 0] = 1 # Noise pixel
-    input_arr[3, 3] = 0 # Hole in the center
+    input_arr[0, 0] = 1  # Noise pixel
+    input_arr[3, 3] = 0  # Hole in the center
 
     print("\nInput Pattern:")
     print_grid(input_arr)
@@ -51,5 +53,6 @@ def main():
     print_grid(closed)
     print("(Note: Hole at [3,3] is filled, noise at [0,0] remains)")
 
+
 if __name__ == "__main__":
     main()

examples/processes_examples.py

@@ -108,7 +108,9 @@ def example_4d():
 def example_chain():
     cube = np.random.rand(12, 64, 64)
     ma = moving_average_temporal(cube, window=5, mode="same")
-    stretched = pixelwise_transform(ma, scale=1.2, offset=-0.1, clamp_min=0.0, clamp_max=1.0)
+    stretched = pixelwise_transform(
+        ma, scale=1.2, offset=-0.1, clamp_min=0.0, clamp_max=1.0
+    )
     print("\nChaining example:")
     print("Input cube shape:", cube.shape)
     print("Moving average shape:", ma.shape)
@@ -150,7 +152,9 @@ def example_perf():
     assert np.allclose(rust_out, naive_out, atol=1e-12)
     print("\nPerformance (1D series length 200k, window=21):")
     print(f"Rust moving_average_temporal: {rust_t:.4f}s")
-    print(f"Naive Python version        : {naive_t:.4f}s  (speedup ~{naive_t / rust_t:.2f}x)")
+    print(
+        f"Naive Python version        : {naive_t:.4f}s  (speedup ~{naive_t / rust_t:.2f}x)"
+    )
 
 
 # ---------------------------------------------------------------------------
@@ -194,7 +198,9 @@ def example_dask():
 # ---------------------------------------------------------------------------
 def example_pixelwise():
     arr = np.array([[0.05, 0.5, 1.2], [0.8, -0.3, 0.4]])
-    scaled = pixelwise_transform(arr, scale=1.5, offset=-0.1, clamp_min=0.0, clamp_max=1.0)
+    scaled = pixelwise_transform(
+        arr, scale=1.5, offset=-0.1, clamp_min=0.0, clamp_max=1.0
+    )
     print("\nPixelwise transform:")
     print("Input:\n", arr)
     print("Scaled & clamped:\n", scaled)

examples/random_forest_example.py

@@ -6,9 +6,10 @@
 from eo_processor import random_forest_predict
 
 # Add the parent directory to the Python path to allow importing from `tests`
-sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..')))
+sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), "..")))
 from tests.utils import sklearn_to_json
 
+
 def main():
     """Main function to run the random forest example."""
     # Generate synthetic data
@@ -37,5 +38,6 @@ def main():
     accuracy = np.mean(predictions == sklearn_predictions)
     print(f"Accuracy compared to scikit-learn: {accuracy * 100:.2f}%")
 
+
 if __name__ == "__main__":
     main()