FIX AFRI formula and reconcile the bare soil index with its reference

docs/spectral-indices.md

@@ -23,10 +23,10 @@ By calculating normalized differences, empirical scaling offsets, and non-linear
          ┌──────────────────────────┴──────────────────────────┐
          ▼                                                     ▼
  [Single-Band Processing Arrays]                       [Multi-Band Raster Math Engine]
- ├─ NDVICalculator (NIR, Red)                           ├─ BICalculator (NIR, Red, Green)
- ├─ NDWICalculator (Green, NIR)                         └─ SAVICalculator (NIR, Red, $L=0.5$)
- ├─ UICalculator (SWIR2, NIR)                          
- └─ AFRICalculator (NIR, SWIR1)                        
+ ├─ NDVICalculator (NIR, Red)                           ├─ BICalculator / BSI
+ ├─ NDWICalculator (Green, NIR)                         │    (SWIR1, Red, NIR, Blue)
+ ├─ UICalculator (SWIR2, NIR)                          └─ SAVICalculator (NIR, Red, $L=0.5$)
+ └─ AFRICalculator (NIR, SWIR1 | SWIR2)                
                                     │
                                     ▼
                     ┌────────────────────────────────┐
@@ -81,23 +81,27 @@ In sparse landscapes, variations in soil moisture, organic matter, and roughness
 
 The $L = 0.5$ adjustment factor shifts the intersection of the soil line back to the coordinate origin, minimizing the effect of background soil brightness. The multiplicative scaler $(1 + L) = 1.5$ ensures the final output remains comparable to the standard $[-1.0, 1.0]$ range.
 
-### AFRICalculator` (Aerosol Free Vegetation Index)
+### `AFRICalculator` (Aerosol Free Vegetation Index)
 
 #### Scientific and Physical Objective
 
 The AFRI is designed to map dense forest canopies and high-biomass woody vegetation while providing a path to bypass atmospheric aerosol scattering (like smoke, haze, or dust). It enhances structural forest signatures while reducing sensitivity to variations in solar illumination, terrain shadowing, and background soil signatures by utilizing the short-wave infrared band instead of visible red.
 
 #### Mathematical Formulation
 
-The calculation separates the input into two distinct, interacting factors:
+Karnieli et al. (2001) define two formulations, selected with the `variant` parameter:
 
-$$\text{AFRI} = (\text{NIR} - 0.66) \times \left( \frac{\text{SWIR1}}{\text{NIR} + 0.66 \times \text{SWIR1}} \right)$$
+$$\text{AFRI}_{1.6} = \frac{\text{NIR} - 0.66 \times \text{SWIR1.6}}{\text{NIR} + 0.66 \times \text{SWIR1.6}} \qquad \text{AFRI}_{2.1} = \frac{\text{NIR} - 0.50 \times \text{SWIR2.1}}{\text{NIR} + 0.50 \times \text{SWIR2.1}}$$
 
 #### Biophysical Interaction Properties
 
-- **NIR Offset Constant ($\text{NIR} - 0.66$):** Dense, healthy forest canopies consistently exhibit high near-infrared reflectance. The empirical constant **0.66** serves as a structural threshold; pixels with low near-infrared values (such as open water, shadows, or asphalt) produce negative or near-zero results, effectively suppressing non-vegetated features.
-
-- **Non-Linear Modulation Ratio:** The second term uses the short-wave infrared band ($\text{SWIR1}$) to modulate the index response. Because moisture-rich forest leaf canopies absorb $\text{SWIR1}$ energy while reflecting $\text{NIR}$, this ratio stays small but positive for healthy forests. This dampens variations caused by topographic shadows, ensuring consistent canopy mapping across rugged terrain.
+- **SWIR Substitution:** AFRI is structurally NDVI with a SWIR band substituted for the visible red. Aerosol scattering is strongly wavelength dependent and falls off toward longer wavelengths, so a SWIR-based index sees through smoke, haze and dust that would depress a red-based index. This is what makes AFRI usable over biomass-burning plumes and dust-laden atmospheres where NDVI fails.
+
+- **The Coefficients ($0.66$ and $0.50$):** These are **scaling factors applied to the SWIR band**, not thresholds. They come from the empirical reflectance relationships $\rho_{0.645} \approx 0.66 \times \rho_{1.6}$ and $\rho_{0.469} \approx 0.50 \times \rho_{2.1}$ reported in the source paper, and their role is to place the SWIR reflectance on the scale of the visible band it replaces. Applying the coefficient anywhere other than to the SWIR reflectance makes the expression dimensionally incoherent.
+
+- **Consistency With NDVI:** Under clear-sky conditions Karnieli et al. report that AFRI and NDVI values are almost identical. That equivalence is a practical validation check: on a haze-free scene, an AFRI implementation that does **not** correlate strongly and positively with NDVI is computing something else.
+
+> **Changed in 1.4.0.** Earlier releases computed $(\text{NIR} - 0.66) \times \text{SWIR1} / (\text{NIR} + 0.66 \times \text{SWIR1})$, which subtracts a bare dimensionless constant from a reflectance and then multiplies by a band ratio. It is a different quantity from Karnieli's index, correlating with it at only 0.17, and it produced negative values over 96% of the bundled example while this page claimed a $[0, +1]$ range. AFRI results from earlier versions are not comparable with current output.
 
 ### `NDWICalculator` (Normalized Difference Water Index)
 
@@ -123,17 +127,29 @@ In contrast, land features like healthy vegetation or dry soils reflect much mor
 
 #### Scientific and Physical Objective
 
-The Bare Soil Index isolates exposed soil surfaces, agricultural fallow fields, mining areas, and bare rock outcroppings by contrasting visible light combinations with near-infrared reflectance.
+The Bare Soil Index isolates exposed soil surfaces, agricultural fallow fields, mining areas, and bare rock outcroppings by contrasting short-wave infrared and red brightness against near-infrared and blue reflectance.
 
 #### Mathematical Formulation
 
-The index is calculated using a specialized normalized difference layout:
+The default formulation is the Bare Soil Index (BSI), selected automatically when `swir1_path` and `blue_path` are supplied:
 
-$$\text{BI} = \frac{(\text{NIR} - \text{Green}) - \text{Red}}{(\text{NIR} + \text{Green} + \text{Red})} = \frac{\text{NIR} - \text{Green} - \text{Red}}{\text{NIR} + \text{Green} + \text{Red}}$$
+$$\text{BSI} = \frac{(\text{SWIR1} + \text{Red}) - (\text{NIR} + \text{Blue})}{(\text{SWIR1} + \text{Red}) + (\text{NIR} + \text{Blue})}$$
 
 #### Biophysical Interaction Properties
 
- **Functional Inversion Reference:** Many published geological studies invert this formula to produce positive values for exposed soils. In this specific implementation, the layout uses $\text{NIR} - (\text{Green} + \text{Red})$ in the numerator. As a result, dense vegetation yields positive values, while bare soils—which show similar reflectance values across the visible green, red, and near-infrared bands—cluster near zero or drop into negative values. Concrete and asphalt structures typically produce strong negative values due to low near-infrared reflectance relative to visible wavelengths.
+- **Four-Band Contrast:** Bare soil and exposed rock are bright in $\text{SWIR1}$ and $\text{Red}$ and comparatively dark in $\text{NIR}$ and $\text{Blue}$. Vegetation is the exact inverse — a strong $\text{NIR}$ plateau against low visible reflectance. Pairing the bands this way makes the numerator change sign between the two surface types, so **BSI is positive over exposed ground and negative over canopy**, which is what allows lithological exposure to be separated from vegetation cover.
+
+- **Why SWIR Is Required:** The $\text{SWIR1}$ term carries the soil-moisture and clay-mineral response that no visible-band combination reproduces. An index built only from visible and NIR bands cannot distinguish dry bare soil from a sparsely vegetated surface of similar visible brightness.
+
+#### Legacy Formulation
+
+Previous releases computed a different expression, which remains reachable by passing `green_path` instead of `swir1_path`/`blue_path`:
+
+$$\text{BI}_{\text{legacy}} = \frac{\text{NIR} - \text{Green} - \text{Red}}{\text{NIR} + \text{Green} + \text{Red}}$$
+
+This is **not a published bare-soil index**. It subtracts two visible bands from NIR, so it responds primarily to vegetation brightness — dense vegetation yields positive values and bare soils cluster near zero or negative, the opposite of what the name implies. It is retained so existing results stay reproducible, and it emits a `DeprecationWarning`. Use the BSI formulation for new work.
+
+> **Citation note.** Earlier documentation attributed this module to As-syakur et al. (2012), which defines EBBI, $(\text{SWIR} - \text{NIR}) / (10\sqrt{\text{SWIR} + \text{TIR}})$ — a different formula requiring a thermal band this tool does not accept. That citation has been removed. The legacy expression above has no published source and should be treated as an original formulation.
 
 ### `UICalculator` (Urban Index)
 
@@ -163,9 +179,9 @@ The following matrix cross-references the required sensor channels, target range
 |---|---|---|---|---|---|---|
 |**`NDVICalculator`**|Canopy Health & Density|$\frac{\text{NIR} - \text{Red}}{\text{NIR} + \text{Red}}$|$[-1.0, \,\, +1.0]$|B5, B4|B8, B4|`'RdYlGn'` / `'YlGn'`|
 |**`SAVICalculator`**|Sparse / Arid Shrublands|$\frac{\text{NIR} - \text{Red}}{\text{NIR} + \text{Red} + 0.5} \times 1.5$|$[-1.0, \,\, +1.0]$|B5, B4|B8, B4|`'RdYlGn'` / `'YlGn'`|
-|**`AFRICalculator`**|Dense / Woody Forests|$(\text{NIR} - 0.66) \times \frac{\text{SWIR1}}{\text{NIR} + 0.66 \cdot \text{SWIR1}}$|$[0.0, \,\, +1.0]$|B5, B6|B8, B11|`'YlGn'`|
+|**`AFRICalculator`**|Dense / Woody Forests, Hazy Scenes|$\frac{\text{NIR} - 0.66 \cdot \text{SWIR1}}{\text{NIR} + 0.66 \cdot \text{SWIR1}}$|$[-1.0, \,\, +1.0]$|B5, B6|B8, B11|`'YlGn'`|
 |**`NDWICalculator`**|Water Bodies & Hydrology|$\frac{\text{Green} - \text{NIR}}{\text{Green} + \text{NIR}}$|$[-1.0, \,\, +1.0]$|B3, B5|B3, B8|`'Blues'`|
-|**`BICalculator`**|Bare Soil & Exposed Rock|$\frac{\text{NIR} - \text{Green} - \text{Red}}{\text{NIR} + \text{Green} + \text{Red}}$|$[-1.0, \,\, +1.0]$|B5, B4, B3|B8, B4, B3|`'inferno'` / `'hot'`|
+|**`BICalculator`**|Bare Soil & Exposed Rock|$\frac{(\text{SWIR1} + \text{Red}) - (\text{NIR} + \text{Blue})}{(\text{SWIR1} + \text{Red}) + (\text{NIR} + \text{Blue})}$|$[-1.0, \,\, +1.0]$|B6, B4, B5, B2|B11, B4, B8, B2|`'inferno'` / `'hot'`|
 |**`UICalculator`**|Built-up Urban Areas|$\frac{\text{SWIR2} - \text{NIR}}{\text{SWIR2} + \text{NIR}}$|$[-1.0, \,\, +1.0]$|B7, B5|B12, B8|`'coolwarm'`|
 
 ## Operational Implementation Examples

fezrs/tools/spectral_indices/afri_calculator.py

@@ -1,28 +1,91 @@
-# Import packages and libraries
-
 # Import module and files
 from fezrs.base import BaseTool
 from fezrs.tools.spectral_indices._division import divide_with_nan
-from fezrs.utils.type_handler import BandPathType
+from fezrs.utils.type_handler import AFRIVariantType, BandPathType
+
+
+# Coefficients from Karnieli et al. (2001), doi:10.1016/S0034-4257(01)00190-0.
+# Each derives from an empirical relationship between visible and SWIR
+# reflectance over vegetated surfaces: rho_0.645 ~= 0.66 * rho_1.6 and
+# rho_0.469 ~= 0.50 * rho_2.1. The coefficient only has meaning as a scaling
+# factor applied to the SWIR band.
+AFRI_COEFFICIENTS = {"1.6": 0.66, "2.1": 0.50}
+
+AFRI_BANDS = {"1.6": "swir1", "2.1": "swir2"}
 
 
 # Calculator class
 class AFRICalculator(BaseTool):
-    def __init__(self, nir_path: BandPathType, swir1_path: BandPathType):
-        super().__init__(nir_path=nir_path, swir1_path=swir1_path)
+    """
+    Aerosol Free Vegetation Index (Karnieli et al., 2001).
+
+    AFRI substitutes a SWIR band for the visible red used by NDVI. Aerosol
+    scattering is strongly wavelength dependent and falls off toward longer
+    wavelengths, so a SWIR-based index penetrates smoke, haze and dust that
+    would otherwise depress a red-based index. Under clear-sky conditions
+    Karnieli et al. report that AFRI and NDVI are almost identical, which is a
+    useful check on any implementation.
+
+    Two formulations are defined:
+
+        AFRI_1.6 = (NIR - 0.66 * SWIR1.6) / (NIR + 0.66 * SWIR1.6)
+        AFRI_2.1 = (NIR - 0.50 * SWIR2.1) / (NIR + 0.50 * SWIR2.1)
+    """
+
+    def __init__(
+        self,
+        nir_path: BandPathType,
+        swir1_path: BandPathType | None = None,
+        swir2_path: BandPathType | None = None,
+        variant: AFRIVariantType = "1.6",
+    ):
+        """
+        Args:
+            nir_path: Near-infrared band.
+            swir1_path: SWIR ~1.6 um band, required for the ``"1.6"`` variant.
+            swir2_path: SWIR ~2.1 um band, required for the ``"2.1"`` variant.
+            variant: Which AFRI formulation to compute.
+        """
+        if variant not in AFRI_COEFFICIENTS:
+            raise ValueError(
+                f"Invalid AFRI variant: {variant!r}. "
+                f"Must be one of {sorted(AFRI_COEFFICIENTS)}."
+            )
+
+        self.variant: AFRIVariantType = variant
+        self.coefficient = AFRI_COEFFICIENTS[variant]
+        self.swir_band = AFRI_BANDS[variant]
+
+        required_path = swir1_path if variant == "1.6" else swir2_path
+        if required_path is None:
+            raise ValueError(
+                f"AFRI variant {variant!r} requires "
+                f"{'swir1_path' if variant == '1.6' else 'swir2_path'}."
+            )
+
+        band_paths = {"nir_path": nir_path}
+        if swir1_path is not None:
+            band_paths["swir1_path"] = swir1_path
+        if swir2_path is not None:
+            band_paths["swir2_path"] = swir2_path
+
+        super().__init__(**band_paths)
+
         self.normalized_bands = self.files_handler.get_normalized_bands(
-            requested_bands=["nir", "swir1"]
+            requested_bands=["nir", self.swir_band]
         )
 
     def _validate(self):
         pass
 
     def process(self):
-        nir, swir1 = (self.normalized_bands[band] for band in ("nir", "swir1"))
+        nir, swir = (
+            self.normalized_bands[band] for band in ("nir", self.swir_band)
+        )
 
-        self._output = (nir - 0.66) * divide_with_nan(
-            swir1,
-            nir + (0.66 * swir1),
+        self._output = divide_with_nan(
+            nir - self.coefficient * swir,
+            nir + self.coefficient * swir,
         )
         return self._output
 

fezrs/tools/spectral_indices/bi_calculator.py

@@ -1,34 +1,126 @@
+import warnings
+
 # Import module and files
 from fezrs.base import BaseTool
 from fezrs.tools.spectral_indices._division import divide_with_nan
-from fezrs.utils.type_handler import BandPathType
+from fezrs.utils.type_handler import BandPathType, BIFormulationType
 
 
 # Calculator class
 class BICalculator(BaseTool):
+    """
+    Bare soil / exposed rock index.
+
+    Two formulations are available.
+
+    ``"bsi"`` (default when SWIR1 and Blue are supplied) is the Bare Soil Index:
+
+        BSI = ((SWIR1 + Red) - (NIR + Blue)) / ((SWIR1 + Red) + (NIR + Blue))
+
+    It works because bare soil and exposed rock are bright in SWIR1 and Red and
+    dark in NIR and Blue, which is the inverse of the vegetation response. That
+    four-band contrast is what separates lithological exposure from canopy, and
+    it is the standard instrument for the target this tool documents.
+
+    ``"legacy"`` is the expression previous releases computed:
+
+        BI_legacy = (NIR - Green - Red) / (NIR + Green + Red)
+
+    It is retained so existing workflows keep running and remain reproducible,
+    but it is not a published bare-soil index: it subtracts two visible bands
+    from NIR, so it responds primarily to vegetation brightness rather than to
+    soil or rock exposure. Prefer ``"bsi"`` for new work.
+    """
+
     def __init__(
         self,
         nir_path: BandPathType,
         red_path: BandPathType,
-        green_path: BandPathType,
+        green_path: BandPathType | None = None,
+        swir1_path: BandPathType | None = None,
+        blue_path: BandPathType | None = None,
+        formulation: BIFormulationType | None = None,
     ):
-        super().__init__(nir_path=nir_path, red_path=red_path, green_path=green_path)
+        """
+        Args:
+            nir_path: Near-infrared band.
+            red_path: Red band.
+            green_path: Green band, required by the ``"legacy"`` formulation.
+            swir1_path: SWIR ~1.6 um band, required by ``"bsi"``.
+            blue_path: Blue band, required by ``"bsi"``.
+            formulation: ``"bsi"`` or ``"legacy"``. Inferred from the supplied
+                bands when omitted.
+        """
+        if formulation is None:
+            formulation = (
+                "bsi" if swir1_path is not None and blue_path is not None else "legacy"
+            )
+
+        if formulation not in ("bsi", "legacy"):
+            raise ValueError(
+                f"Invalid BI formulation: {formulation!r}. "
+                "Must be 'bsi' or 'legacy'."
+            )
+
+        self.formulation: BIFormulationType = formulation
+
+        if formulation == "bsi":
+            if swir1_path is None or blue_path is None:
+                raise ValueError(
+                    "The 'bsi' formulation requires both swir1_path and blue_path."
+                )
+            self._required_bands = ("swir1", "red", "nir", "blue")
+        else:
+            if green_path is None:
+                raise ValueError(
+                    "The 'legacy' formulation requires green_path."
+                )
+            warnings.warn(
+                "BICalculator's 'legacy' formulation "
+                "(NIR - Green - Red) / (NIR + Green + Red) is not a published "
+                "bare-soil index and responds mainly to vegetation brightness. "
+                "Pass swir1_path and blue_path to compute BSI instead.",
+                DeprecationWarning,
+                stacklevel=2,
+            )
+            self._required_bands = ("nir", "red", "green")
+
+        band_paths = {"nir_path": nir_path, "red_path": red_path}
+        if green_path is not None:
+            band_paths["green_path"] = green_path
+        if swir1_path is not None:
+            band_paths["swir1_path"] = swir1_path
+        if blue_path is not None:
+            band_paths["blue_path"] = blue_path
+
+        super().__init__(**band_paths)
+
         self.normalized_bands = self.files_handler.get_normalized_bands(
-            requested_bands=["nir", "red", "green"]
+            requested_bands=list(self._required_bands)
         )
 
     def _validate(self):
         pass
 
     def process(self):
-        nir, red, green = (
-            self.normalized_bands[band] for band in ("nir", "red", "green")
-        )
+        if self.formulation == "bsi":
+            swir1, red, nir, blue = (
+                self.normalized_bands[band]
+                for band in ("swir1", "red", "nir", "blue")
+            )
+            self._output = divide_with_nan(
+                (swir1 + red) - (nir + blue),
+                (swir1 + red) + (nir + blue),
+            )
+        else:
+            nir, red, green = (
+                self.normalized_bands[band] for band in ("nir", "red", "green")
+            )
+            self._output = divide_with_nan(
+                (nir - green) - red,
+                (nir + green) + red,
+            )
 
-        self._output = divide_with_nan(
-            (nir - green) - red,
-            (nir + green) + red,
-        )
         return self._output
 
     def execute(

fezrs/utils/type_handler.py

@@ -101,3 +101,10 @@ class BandPathsType(TypedDict, total=False):
     "direction",
 ]
 """Type alias for magnitude/direction change detection."""
+
+
+AFRIVariantType = Literal["1.6", "2.1"]
+"""Type alias for the two AFRI formulations defined by Karnieli et al. (2001)."""
+
+BIFormulationType = Literal["bsi", "legacy"]
+"""Type alias for the supported bare-soil index formulations."""