diff --git a/.github/workflows/deploy-docs.yml b/.github/workflows/deploy-docs.yml
index 19697b69..a77726bd 100644
--- a/.github/workflows/deploy-docs.yml
+++ b/.github/workflows/deploy-docs.yml
@@ -4,6 +4,11 @@ on:
   release:
     types: [published]
   workflow_dispatch:
+    inputs:
+      version:
+        description: "Documentation version to deploy"
+        required: true
+        type: string
 
 permissions:
   contents: write
@@ -38,6 +43,10 @@ jobs:
           git config --global user.email 'github-actions[bot]@users.noreply.github.com'
 
       - name: "Deploy Github Pages"
+        env:
+          DOCS_VERSION: ${{ github.event.release.tag_name || inputs.version }}
         run: |
-          mike deploy --push --update-aliases ${{ github.event.release.tag_name }} latest
-          mike set-default --push latest
+          mike delete latest
+          mike deploy --update-aliases "$DOCS_VERSION" latest
+          mike set-default latest
+          git push origin gh-pages
diff --git a/docs/api/concepts/holistic_craft.md b/docs/api/concepts/holistic_craft.md
index 6d5da5f0..1db037e4 100644
--- a/docs/api/concepts/holistic_craft.md
+++ b/docs/api/concepts/holistic_craft.md
@@ -64,8 +64,10 @@ Holistic CRAFT follows the same core principle as CRAFT but operates on full ima
     Third-party NMF implementations may not have this limitation
     (e.g., the Semi-NMF from the Overcomplete library).
 
-3. **Estimate Importance**: Use any attribution methods available in Xplique (gradient-based, perturbation-based) to rank concept importance
-4. **Visualize**: Generate concept heatmaps overlaid on images to show "what" and "where"
+3. **Measure Concept Activation**: Use concept coefficients to measure how strongly each concept is present and rank representative images
+4. **Estimate Concept Importance**: Attribute task predictions to concept coefficients to measure how much each concept contributes to the model output
+5. **Visualize Coefficients**: Resize latent coefficient maps and overlay them on input images
+6. **Localize Concepts**: Optionally use a black-box attribution method to identify which input regions drive a selected concept score
 
 Like regular CRAFT, Holistic CRAFT requires splitting the model into two parts: $(g, h)$ such that $f(x) = (g \cdot h)(x)$. The model $g$ maps input to latent space (activation maps), and $h$ maps latent space to predictions. Concepts are extracted from these activation maps in latent space.
 
@@ -77,6 +79,37 @@ This split is implemented through three abstractions:
 
 - **`LatentExtractorBuilder`**: A factory that constructs a `LatentExtractor` for a specific model architecture. It handles all the architecture-specific wiring (defining how to split the model, which layer to extract from, and how to format outputs) so that the rest of the CRAFT pipeline remains model-agnostic.
 
+## Concept Activation, Importance, and Localization
+
+Holistic CRAFT exposes three related quantities that answer different questions:
+
+| Quantity | Definition | Question answered |
+|---|---|---|
+| **Concept activation** | The spatial coefficient map $U_k(x)$, or a reduction of it | How strongly is concept $k$ present? |
+| **Concept importance** | Attribution of the final task prediction to concept $k$ | How much does concept $k$ contribute to the prediction? |
+| **Concept localization** | Attribution of the concept score $s_k(x) = R(U_k(x))$ to the input | Which input regions drive the concept score? |
+
+Concept importance and concept localization follow opposite attribution directions:
+
+```text
+concept coefficients -> final prediction -> concept importance
+
+input image -> concept scores -> input attribution -> concept localization
+```
+
+A coefficient heatmap is a latent-space map resized to the input resolution. A localization
+map is an input-space attribution produced by perturbing the input and observing changes in a
+selected concept score. The maps can therefore differ, especially when the latent
+representation is coarse or its positions have large receptive fields.
+
+Top-image ranking always uses concept coefficients because they represent concept presence.
+Attribution-map magnitude is not a replacement for concept activation.
+
+!!!note
+    Coefficient and localization maps are explanations, not segmentation masks. Black-box
+    localization is not automatically more correct than coefficient visualization; it answers
+    the more specific question of which input perturbations change a selected concept score.
+
 
 ## Example
 
@@ -90,8 +123,8 @@ from xplique_adapters.concepts.torch.latent_data_retinanet import RetinanetExtra
 # Build a latent extractor that splits the model into g(.) and h(.)
 # This provides the input_to_latent (g) and latent_to_logit (h) functions
 latent_extractor = RetinanetExtractorBuilder.build(
-    model, 
-    device="cuda", 
+    model,
+    device="cuda",
     nb_classes=91,
     extraction_location='resnet', # Choose 'resnet' or 'fpn'
     extraction_layer=-1  # Extract from last ResNet feature layer
@@ -162,14 +195,167 @@ explanation_vargrad = craft.compute_explanation_per_concept(
     confidence=0.3,
 )
 
-# Reduce the spatial dimension of the explanation 
+# Reduce the spatial dimension of the explanation
 # to compute the final concepts importances
 importances_vargrad = craft.reduce_to_importance(
     explanation=explanation_vargrad,
 )
 ```
 
-### Using a Different NMF Factorizer
+## Localizing Concepts with Black-Box Attribution
+
+`ConceptLocalizer` exposes the fitted encoder and factorizer as a callable that returns one
+scalar score per learned concept. `compute_concept_attributions()` builds this callable,
+constructs one-hot concept targets, and applies a compatible black-box explainer to each
+requested concept.
+
+For an input batch with shape `(N, H, W, C)`, the returned maps have shape
+`(N, H, W, number_of_concepts)`. Channel `k` always corresponds to concept `k`. When only a
+subset is requested, uncomputed channels are filled entirely with `NaN` so that they cannot
+be confused with valid zero attributions.
+
+TensorFlow callers pass channel-last images. PyTorch callers pass their native channel-first
+`(N, C, H, W)` images to the same high-level method; Xplique handles the layout conversion
+for the wrapped localizer.
+
+### Localizing Concepts with RISE
+
+```python
+from xplique.attributions import Rise
+from xplique.concepts import PartialExplainer
+
+rise = PartialExplainer(
+    Rise,
+    nb_samples=2000,
+    grid_size=7,
+    preservation_probability=0.5,
+    mask_value=0.0,
+)
+
+rise_maps = craft.compute_concept_attributions(
+    images,
+    partial_explainer=rise,
+    concept_ids=[0, 3, 7],
+    concept_reducer="mean",
+)
+
+craft.display_images_per_concept(
+    display_images,
+    concept_maps=rise_maps,
+    order=[0, 3, 7],
+)
+```
+
+`concept_reducer="mean"` reduces every spatial coefficient map to the scalar score that RISE
+attributes to the input. It is the default and is consistent with the mean coefficient
+activation used to rank representative images. The localizer preserves signed scores; a
+custom callable reducer can be used when concept magnitude is intended instead.
+
+Do not pass a task `operator` to `PartialExplainer` for concept localization. One-hot targets
+are created internally to select concept scores directly.
+
+### Localizing Concepts with Sobol
+
+The same API works with `SobolAttributionMethod`:
+
+```python
+from xplique.attributions import SobolAttributionMethod
+
+sobol = PartialExplainer(
+    SobolAttributionMethod,
+    grid_size=8,
+    nb_design=32,
+    perturbation_function="inpainting",
+)
+
+sobol_maps = craft.compute_concept_attributions(
+    images,
+    partial_explainer=sobol,
+    concept_ids=[0, 3, 7],
+)
+
+craft.display_images_per_concept(
+    display_images,
+    concept_maps=sobol_maps,
+    order=[0, 3, 7],
+)
+```
+
+`nb_design` must be a nonzero power of two. Keep Sobol's default `nb_channels=1`: concept
+localization computes a separate single-channel attribution map for each selected concept.
+This differs from Sobol-based concept importance, which attributes the final task prediction
+to concept coefficients.
+
+### Coefficient Maps and Localization Maps
+
+Use `coeffs_u` to display latent coefficient maps:
+
+```python
+craft.display_images_per_concept(
+    display_images,
+    coeffs_u=coefficients,
+    order=selected_concepts,
+)
+```
+
+Use `concept_maps` to display input-space localization maps:
+
+```python
+craft.display_images_per_concept(
+    display_images,
+    concept_maps=rise_maps,
+    order=selected_concepts,
+)
+```
+
+When displaying the top images, coefficients still determine the ranking and localization
+maps only determine the overlay:
+
+```python
+craft.display_top_images_per_concept(
+    display_images,
+    coeffs_u=coefficients,
+    concept_maps=rise_maps,
+    order=selected_concepts,
+    topk=3,
+)
+```
+
+If `order` is omitted, display methods show only the concept channels that were computed.
+Explicitly requesting an uncomputed `NaN` channel raises a `ValueError`. Signed localization
+maps are displayed by absolute magnitude; attribution direction is not represented by the
+current renderer.
+
+!!!warning "Computational cost"
+    Black-box localization evaluates the encoder and `factorizer.encode()` for many perturbed
+    inputs and runs a separate attribution pass for every selected concept. Estimate concept
+    importance first, then localize only a small number of important concepts. Reduce the
+    number of images, RISE samples, Sobol designs, or grid resolution for exploratory runs.
+
+!!!warning "Perturbations use preprocessed inputs"
+    Perturbation parameters operate in the model's preprocessed input space. For example,
+    `mask_value=0.0` is neutral or black only if zero has that meaning after preprocessing.
+    With standard ImageNet normalization, zero generally corresponds to the dataset mean
+    rather than a raw black pixel. Sobol inpainting and blurring must be interpreted relative
+    to the same model-ready representation.
+
+!!!warning "Factorizer compatibility"
+    Localization evaluates the fitted factorizer on unseen, perturbed activations. The
+    factorizer must therefore support out-of-sample `encode()`. Activations should not be
+    clipped merely to satisfy a factorizer because that would change the function being
+    explained.
+
+!!!warning "Black-box scope"
+    Input-to-concept localization currently supports black-box attribution methods only.
+    Ordinary `factorizer.encode()` is not guaranteed to be differentiable, so white-box
+    explainers such as Gradient Input or Integrated Gradients are rejected.
+
+!!!warning "Interpretation"
+    Localization maps measure sensitivity of a concept score to input perturbations. They are
+    not object boundaries or segmentation labels. RISE and Sobol can also produce different
+    maps because they use different perturbation and aggregation strategies.
+
+## Using a Different NMF Factorizer
 
 By default, the standard Sklearn NMF is used to factorize the concepts.
 But other types of factorizers are supported, such as the ones provided
@@ -220,7 +406,7 @@ class CustomLatentData(LatentData):
     def get_activations(self, as_numpy: bool = True, keep_gradients: bool = False):
         """Extract activations from the specified layer."""
         activations = self.fpn_outs[self.extraction_layer]
-        
+
         if not keep_gradients:
             activations = activations.detach()
 
@@ -264,7 +450,7 @@ class CustomExtractorBuilder(LatentExtractorBuilder):
         extraction_layer: int = -1,
         batch_size: int = 1
     ) -> TorchLatentExtractor:
-        
+
         # Define g(.) function: input → latent activations
         def g(self, x):
             # Example: extract from backbone/feature pyramid
@@ -297,7 +483,7 @@ class CustomExtractorBuilder(LatentExtractorBuilder):
             batch_size=batch_size,
             device=device
         )
-        
+
         # Store extraction layer for later use
         latent_extractor.extraction_layer = extraction_layer
         return latent_extractor
@@ -344,6 +530,8 @@ craft.display_images_per_concept(input_images[:5])
 
 {{xplique.concepts.holistic_craft.PartialExplainer}}
 
+{{xplique.concepts.holistic_craft.ConceptLocalizer}}
+
 ## References
 
 [^1]: [CRAFT: Concept Recursive Activation FacTorization for Explainability (2023).](https://arxiv.org/pdf/2211.10154.pdf)
diff --git a/tests/concepts/test_holistic_craft_classification_tf.py b/tests/concepts/test_holistic_craft_classification_tf.py
index c1db481b..7297022f 100644
--- a/tests/concepts/test_holistic_craft_classification_tf.py
+++ b/tests/concepts/test_holistic_craft_classification_tf.py
@@ -7,15 +7,34 @@
 from PIL import Image
 
 import xplique
-from xplique.attributions import Saliency
+from xplique.attributions import Rise, Saliency, SobolAttributionMethod
 from xplique.attributions.gradient_input import GradientInput
 from xplique.concepts import HolisticCraftTf as Craft
 from xplique.concepts.holistic_craft import PartialExplainer
-from xplique.concepts.tf.layered_model_latent_extractor import LayeredModelExtractorBuilder
+from xplique.concepts.tf.latent_extractor import TfLatentExtractor
+from xplique.concepts.tf.layered_model_latent_extractor import (
+    LayeredLatentData,
+    LayeredModelExtractorBuilder,
+)
 from xplique.utils_functions.classification.tf.classifier_tensor import TfClassifierTensor
 from xplique.utils_functions.common.tf.gradients_check import check_model_gradients
 
 
+class _IdentityFactorizer:
+    is_fitted = False
+    requires_positive_activations = False
+
+    def fit(self, activations):
+        self.is_fitted = True
+        return np.eye(2, dtype=np.float32), np.asarray(activations, dtype=np.float32)
+
+    def encode(self, activations):
+        return np.asarray(activations, dtype=np.float32)
+
+    def encode_differentiable(self, activations):
+        return activations
+
+
 def test_classifier_tensor_targets_class_and_preserves_batch_shape():
     predictions = TfClassifierTensor(tf.constant([[0.1, 0.2, 0.7], [0.3, 0.6, 0.1]]))
 
@@ -134,8 +153,6 @@ def test_latent_extractor(image_data, latent_extractor_data):
 
 def test_layered_latent_data_getitem():
     """Test that LayeredLatentData supports integer and slice indexing."""
-    from xplique.concepts.tf.layered_model_latent_extractor import LayeredLatentData
-
     activations = tf.random.normal((4, 7, 7, 16))
     latent_data = LayeredLatentData(activations)
 
@@ -199,6 +216,23 @@ def craft_data(image_data, latent_extractor_data, device_param):
         return craft
 
 
+@pytest.fixture
+def tiny_craft_data():
+    """Create a deterministic identity CRAFT pipeline for localization tests."""
+    values = np.arange(2 * 4 * 4 * 2, dtype=np.float32).reshape(2, 4, 4, 2)
+    extractor = TfLatentExtractor(
+        model=lambda inputs: inputs,
+        input_to_latent_model=lambda inputs: LayeredLatentData(inputs),
+        latent_to_logit_model=lambda latent_data: tf.reduce_mean(
+            latent_data.activations, axis=(1, 2)
+        ),
+        batch_size=2,
+    )
+    craft = Craft(extractor, number_of_concepts=2, factorizer=_IdentityFactorizer())
+    craft.fit(tf.constant(values))
+    return craft, values
+
+
 def test_craft_reencode(image_data, craft_data):
     """Test CRAFT encode and decode operations."""
     _, input_tensor = image_data
@@ -340,3 +374,48 @@ def test_craft_sobol_importance(image_data, craft_data):
     order = importances_sobol.argsort()[::-1]
     assert order.shape == (10,), "Should have ordering for all concepts"
     assert len(np.unique(order)) == 10, "All concepts should have unique ordering"
+
+
+def test_craft_make_concept_localizer_matches_reduced_transform(tiny_craft_data):
+    craft, images = tiny_craft_data
+
+    localizer = craft.make_concept_localizer("mean")
+    scores = localizer(images).numpy()
+    coeffs_u = craft.transform(images)
+    expected = np.mean(coeffs_u, axis=(1, 2))
+
+    assert scores.shape == (2, craft.number_of_concepts)
+    assert scores.dtype == np.float32
+    assert np.all(np.isfinite(scores))
+    np.testing.assert_allclose(scores, expected, rtol=1e-5, atol=1e-5)
+
+
+@pytest.mark.parametrize("method", ["rise", "sobol"])
+def test_craft_compute_concept_attributions_black_box_smoke(tiny_craft_data, method):
+    craft, images = tiny_craft_data
+    tf.random.set_seed(1)
+    if method == "rise":
+        explainer = PartialExplainer(
+            Rise,
+            nb_samples=8,
+            grid_size=2,
+            preservation_probability=0.5,
+        )
+    else:
+        explainer = PartialExplainer(
+            SobolAttributionMethod,
+            grid_size=2,
+            nb_design=2,
+            perturbation_function="inpainting",
+        )
+
+    maps = craft.compute_concept_attributions(
+        images[:1],
+        partial_explainer=explainer,
+        concept_ids=[0],
+    )
+
+    assert maps.shape == (1, 4, 4, craft.number_of_concepts)
+    assert maps.dtype == np.float32
+    assert np.all(np.isfinite(maps[..., 0]))
+    assert np.all(np.isnan(np.delete(maps, 0, axis=-1)))
diff --git a/tests/concepts/test_holistic_craft_classification_torch.py b/tests/concepts/test_holistic_craft_classification_torch.py
index 495a29e5..51320737 100644
--- a/tests/concepts/test_holistic_craft_classification_torch.py
+++ b/tests/concepts/test_holistic_craft_classification_torch.py
@@ -3,22 +3,47 @@
 
 import numpy as np
 import pytest
+import tensorflow as tf
 import torch
 import torchvision.transforms as T
 from PIL import Image
 from torchvision import models
 
 import xplique
-from xplique.attributions import Saliency
+from xplique.attributions import Rise, Saliency, SobolAttributionMethod
 from xplique.attributions.gradient_input import GradientInput
 from xplique.concepts import HolisticCraftTorch as Craft
 from xplique.concepts.holistic_craft import PartialExplainer
-from xplique.concepts.torch.layered_model_latent_extractor import LayeredModelExtractorBuilder
+from xplique.concepts.tf.holistic_craft import HolisticCraftTf
+from xplique.concepts.tf.latent_extractor import TfLatentExtractor
+from xplique.concepts.tf.layered_model_latent_extractor import LayeredLatentData as TfLatentData
+from xplique.concepts.torch.latent_extractor import TorchLatentExtractor
+from xplique.concepts.torch.layered_model_latent_extractor import (
+    LayeredLatentData as TorchLatentData,
+)
+from xplique.concepts.torch.layered_model_latent_extractor import (
+    LayeredModelExtractorBuilder,
+)
 from xplique.utils_functions.classification.torch.classifier_tensor import TorchClassifierTensor
 from xplique.utils_functions.common.torch.gradients_check import check_model_gradients
 from xplique.wrappers import TorchWrapper
 
 
+class _IdentityFactorizer:
+    is_fitted = False
+    requires_positive_activations = False
+
+    def fit(self, activations):
+        self.is_fitted = True
+        return np.eye(2, dtype=np.float32), np.asarray(activations, dtype=np.float32)
+
+    def encode(self, activations):
+        return np.asarray(activations, dtype=np.float32)
+
+    def encode_differentiable(self, activations):
+        return activations
+
+
 def test_classifier_tensor_targets_class_and_preserves_batch_shape():
     predictions = TorchClassifierTensor.from_predictions(
         torch.tensor([[0.1, 0.2, 0.7], [0.3, 0.6, 0.1]])
@@ -134,17 +159,15 @@ def test_latent_extractor(image_data, latent_extractor_data):
 
 def test_layered_latent_data_getitem(device_param):
     """Test that LayeredLatentData supports integer and slice indexing."""
-    from xplique.concepts.torch.layered_model_latent_extractor import LayeredLatentData
-
     activations = torch.randn(4, 16, 7, 7).to(device_param)
-    latent_data = LayeredLatentData(activations)
+    latent_data = TorchLatentData(activations)
 
     sliced = latent_data[1:3]
-    assert isinstance(sliced, LayeredLatentData)
+    assert isinstance(sliced, TorchLatentData)
     assert sliced.activations.shape[0] == 2
 
     item = latent_data[0]
-    assert isinstance(item, LayeredLatentData)
+    assert isinstance(item, TorchLatentData)
 
 
 def test_latent_extractor_gradients(image_data, latent_extractor_data):
@@ -209,6 +232,32 @@ def craft_data(image_data, latent_extractor_data, device_param):
     return craft
 
 
+def _make_tiny_torch_craft(device):
+    values = np.arange(2 * 4 * 4 * 2, dtype=np.float32).reshape(2, 4, 4, 2)
+    images = torch.from_numpy(values.transpose(0, 3, 1, 2)).to(device)
+    extractor = TorchLatentExtractor(
+        model=torch.nn.Identity(),
+        input_to_latent_model=lambda inputs: TorchLatentData(inputs),
+        latent_to_logit_model=lambda latent_data: latent_data.activations,
+        device=str(device),
+        batch_size=2,
+    )
+    craft = Craft(
+        extractor,
+        number_of_concepts=2,
+        device=str(device),
+        factorizer=_IdentityFactorizer(),
+    )
+    craft.fit(images)
+    return craft, images, values
+
+
+@pytest.fixture
+def tiny_craft_data(device_param):
+    """Create a deterministic identity CRAFT pipeline for localization tests."""
+    return _make_tiny_torch_craft(device_param)
+
+
 def test_craft_reencode(image_data, craft_data):
     """Test CRAFT encode and decode operations."""
     _, input_tensor = image_data
@@ -356,3 +405,107 @@ def test_craft_sobol_importance(image_data, craft_data):
     order = importances_sobol.argsort()[::-1]
     assert order.shape == (10,), "Should have ordering for all concepts"
     assert len(np.unique(order)) == 10, "All concepts should have unique ordering"
+
+
+def test_craft_make_concept_localizer_matches_reduced_transform(tiny_craft_data):
+    craft, images_nchw, images_nhwc = tiny_craft_data
+
+    localizer = craft.make_concept_localizer("mean")
+    scores = localizer(images_nhwc)
+    scores = scores.numpy()
+    coeffs_u = craft.transform(images_nchw)
+    expected = np.mean(coeffs_u, axis=(1, 2))
+
+    assert scores.shape == (2, craft.number_of_concepts)
+    assert scores.dtype == np.float32
+    assert np.all(np.isfinite(scores))
+    np.testing.assert_allclose(scores, expected, rtol=5e-4, atol=2e-3)
+
+
+def test_compute_concept_attributions_normalizes_native_nchw_inputs(tiny_craft_data):
+    craft, images_nchw, _ = tiny_craft_data
+    observed_shapes = []
+
+    class RecordingExplainer:
+        def __init__(self, model, batch_size):
+            del batch_size
+            self.model = model
+
+        def explain(self, inputs, targets):
+            observed_shapes.append(tuple(inputs.shape))
+            scores = self.model(inputs)
+            assert scores.shape == (len(inputs), craft.number_of_concepts)
+            return np.ones(inputs.shape[:3], dtype=np.float32)
+
+    explainer = PartialExplainer(RecordingExplainer)
+    input_variants = [
+        images_nchw[:1],
+        images_nchw[:1].detach().cpu().numpy(),
+        [images_nchw[0]],
+    ]
+    for inputs in input_variants:
+        maps = craft.compute_concept_attributions(inputs, explainer, concept_ids=[0])
+        assert maps.shape == (1, 4, 4, craft.number_of_concepts)
+
+    assert observed_shapes == [(1, 4, 4, 2)] * len(input_variants)
+
+
+@pytest.mark.parametrize("method", ["rise", "sobol"])
+def test_craft_compute_concept_attributions_black_box_smoke(tiny_craft_data, method):
+    craft, images_nchw, _ = tiny_craft_data
+    tf.random.set_seed(1)
+    if method == "rise":
+        explainer = PartialExplainer(
+            Rise,
+            nb_samples=8,
+            grid_size=2,
+            preservation_probability=0.5,
+        )
+    else:
+        explainer = PartialExplainer(
+            SobolAttributionMethod,
+            grid_size=2,
+            nb_design=2,
+            perturbation_function="inpainting",
+        )
+
+    maps = craft.compute_concept_attributions(
+        images_nchw[:1],
+        partial_explainer=explainer,
+        concept_ids=[0],
+    )
+
+    assert maps.shape == (1, 4, 4, craft.number_of_concepts)
+    assert maps.dtype == np.float32
+    assert np.all(np.isfinite(maps[..., 0]))
+    assert np.all(np.isnan(np.delete(maps, 0, axis=-1)))
+
+
+@pytest.mark.parametrize("reducer", ["mean", "sum"])
+def test_tf_torch_concept_localizer_score_parity(reducer):
+    values = np.arange(2 * 4 * 4 * 2, dtype=np.float32).reshape(2, 4, 4, 2)
+    tf_extractor = TfLatentExtractor(
+        model=lambda inputs: inputs,
+        input_to_latent_model=lambda inputs: TfLatentData(inputs),
+        latent_to_logit_model=lambda latent_data: latent_data.activations,
+        batch_size=2,
+    )
+    tf_craft = HolisticCraftTf(
+        tf_extractor,
+        number_of_concepts=2,
+        factorizer=_IdentityFactorizer(),
+    )
+    torch_craft, torch_images, _ = _make_tiny_torch_craft(torch.device("cpu"))
+    tf_images = tf.constant(values)
+    tf_craft.fit(tf_images)
+
+    tf_scores = tf_craft.make_concept_localizer(reducer)(tf_images).numpy()
+    torch_scores = torch_craft.make_concept_localizer(reducer)(values).numpy()
+
+    np.testing.assert_allclose(tf_scores, torch_scores, rtol=1e-6, atol=1e-6)
+    np.testing.assert_allclose(
+        torch_scores,
+        getattr(np, reducer)(torch_craft.transform(torch_images), axis=(1, 2)),
+        rtol=1e-6,
+        atol=1e-6,
+    )
diff --git a/tests/concepts/test_holistic_craft_regressions.py b/tests/concepts/test_holistic_craft_regressions.py
index 3f77818a..78db2550 100644
--- a/tests/concepts/test_holistic_craft_regressions.py
+++ b/tests/concepts/test_holistic_craft_regressions.py
@@ -2,13 +2,17 @@
 
 from contextlib import contextmanager
 
-import matplotlib.pyplot as plt
+import matplotlib
 import numpy as np
 import pytest
 from sklearn.exceptions import NotFittedError
 
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt  # pylint: disable=wrong-import-position
+
+from xplique.attributions.gradient_input import GradientInput
 from xplique.concepts.craft import Factorization
-from xplique.concepts.holistic_craft import HolisticCraft, PartialExplainer
+from xplique.concepts.holistic_craft import ConceptLocalizer, HolisticCraft, PartialExplainer
 from xplique.concepts.latent_extractor import LatentData
 
 
@@ -60,6 +64,18 @@ def latent_to_logit(self, latent_data):
         return _Prediction()
 
 
+class _SemanticExtractor:
+    batch_size = 2
+
+    def input_to_latent_generator(self, inputs, resize=None, keep_gradients=False):
+        del resize, keep_gradients
+        inputs = np.asarray(inputs)
+        concepts = np.zeros((inputs.shape[0], inputs.shape[1], inputs.shape[2], 2))
+        concepts[:, :2, :2, 0] = inputs[:, :2, :2, 0]
+        concepts[:, 2:, 2:, 1] = inputs[:, 2:, 2:, 1]
+        yield _LatentData(concepts.astype(np.float32))
+
+
 class _Factorizer:
     is_fitted = True
     requires_positive_activations = False
@@ -72,6 +88,11 @@ class _UnfittedFactorizer(_Factorizer):
     is_fitted = False
 
 
+class _NonInductiveFactorizer(_Factorizer):
+    def encode(self, activations):
+        raise NotImplementedError("out-of-sample encoding is not supported")
+
+
 class _ArrayLike:
     def __init__(self, values):
         self.values = values
@@ -95,12 +116,22 @@ class _Framework:
 
 
 class _Craft(HolisticCraft):
-    def __init__(self, latent_data, batch_size=1):
-        factorizer = _Factorizer()
+    def __init__(
+        self,
+        latent_data,
+        batch_size=1,
+        number_of_concepts=2,
+        factorizer=None,
+        extractor=None,
+    ):
+        factorizer = factorizer or _Factorizer()
         super().__init__(
-            _Extractor(latent_data, batch_size), number_of_concepts=2, factorizer=factorizer
+            extractor or _Extractor(latent_data, batch_size),
+            number_of_concepts=number_of_concepts,
+            factorizer=factorizer,
         )
-        self.factorization = Factorization(None, 0, None, factorizer, None, np.eye(2))
+        concept_bank = np.eye(number_of_concepts, dtype=np.float32)
+        self.factorization = Factorization(None, 0, None, factorizer, None, concept_bank)
         self.framework = "tf"
         self._framework_module = _Framework
 
@@ -136,6 +167,18 @@ def to_batched_tensor(self):
         return np.ones((1, 1), dtype=np.float32)
 
 
+def _make_spatial_craft(number_of_concepts=3, factorizer=None):
+    activations = np.arange(2 * 2 * 2 * number_of_concepts, dtype=np.float32).reshape(
+        2, 2, 2, number_of_concepts
+    )
+    return _Craft(
+        [_LatentData(activations)],
+        batch_size=2,
+        number_of_concepts=number_of_concepts,
+        factorizer=factorizer,
+    )
+
+
 def test_factorization_preserves_its_positional_field_order():
     factorization = Factorization("inputs", 3, "crops", "reducer", "crops_u", "concept_bank")
 
@@ -290,3 +333,478 @@ def test_display_validates_concept_order_and_handles_a_single_column():
     figure = craft.display_top_images_per_concept(images, topk=2, coeffs_u=coeffs_u, order=[0])
     assert len(figure.axes) == 2
     plt.close(figure)
+
+
+def test_concept_localizer_reducers_handle_spatial_and_global_coefficients():
+    craft = _make_spatial_craft(number_of_concepts=3)
+    localizer = craft.make_concept_localizer("mean")
+    inputs = np.ones((2, 4, 4, 3), dtype=np.float32)
+
+    assert isinstance(localizer, ConceptLocalizer)
+
+    scores_mean = localizer(inputs)
+    expected_mean = np.mean(craft.transform(inputs), axis=(1, 2))
+    np.testing.assert_allclose(scores_mean, expected_mean)
+    assert scores_mean.dtype == np.float32
+
+    scores_sum = craft.make_concept_localizer("sum")(inputs)
+    expected_sum = np.sum(craft.transform(inputs), axis=(1, 2))
+    np.testing.assert_allclose(scores_sum, expected_sum)
+
+    scores_max = craft.make_concept_localizer("max")(inputs)
+    expected_max = np.max(craft.transform(inputs), axis=(1, 2))
+    np.testing.assert_allclose(scores_max, expected_max)
+
+    tokens = np.arange(2 * 4 * 3, dtype=np.float32).reshape(2, 4, 3)
+    token_craft = _Craft([_LatentData(tokens)], batch_size=2, number_of_concepts=3)
+    token_scores = token_craft.make_concept_localizer()(inputs)
+    np.testing.assert_allclose(token_scores, np.mean(tokens, axis=1))
+
+    globals_only = np.arange(2 * 3, dtype=np.float32).reshape(2, 3)
+    global_craft = _Craft([_LatentData(globals_only)], batch_size=2, number_of_concepts=3)
+    global_scores = global_craft.make_concept_localizer()(inputs)
+    np.testing.assert_allclose(global_scores, globals_only)
+
+
+def test_concept_localizer_reducer_validation_and_shape_errors():
+    craft = _make_spatial_craft(number_of_concepts=3)
+
+    with pytest.raises(ValueError, match="concept_reducer"):
+        craft.make_concept_localizer("median")
+
+    with pytest.raises(ValueError, match="concept_reducer"):
+        craft.make_concept_localizer(42)
+
+    callable_localizer = craft.make_concept_localizer(
+        lambda coeffs: np.mean(np.abs(coeffs), axis=(1, 2), dtype=np.float64)
+    )
+    callable_inputs = np.ones((2, 4, 4, 3), dtype=np.float32)
+    callable_scores = callable_localizer(callable_inputs)
+    assert callable_scores.dtype == np.float32
+    np.testing.assert_allclose(
+        callable_scores,
+        np.mean(np.abs(craft.transform(callable_inputs)), axis=(1, 2)),
+    )
+
+    with pytest.raises(ValueError, match="must have shape"):
+        craft.make_concept_localizer(lambda coeffs: np.mean(coeffs, axis=(0, 1, 2)))(
+            np.ones((2, 4, 4, 3), dtype=np.float32)
+        )
+
+    empty_batch_craft = _Craft(
+        [_LatentData(np.empty((0, 3), dtype=np.float32))],
+        number_of_concepts=3,
+    )
+    with pytest.raises(ValueError, match="empty dimensions"):
+        empty_batch_craft.make_concept_localizer()(np.ones((1, 2, 2, 3), dtype=np.float32))
+
+    empty_spatial_craft = _Craft(
+        [_LatentData(np.empty((2, 0), dtype=np.float32))],
+        number_of_concepts=3,
+    )
+    with pytest.raises(ValueError, match="empty dimensions"):
+        empty_spatial_craft.make_concept_localizer("max")(np.ones((1, 2, 2, 3), dtype=np.float32))
+    with pytest.raises(ValueError, match="finite"):
+        craft.make_concept_localizer(lambda coeffs: np.full((2, 3), np.nan))(
+            np.ones((2, 4, 4, 3), dtype=np.float32)
+        )
+
+    invalid_rank_craft = _Craft(
+        [_LatentData(np.array([1.0, 2.0], dtype=np.float32))],
+        number_of_concepts=2,
+    )
+    with pytest.raises(ValueError, match="at least 2 dimensions"):
+        invalid_rank_craft.make_concept_localizer()(np.ones((1, 2, 2, 1), dtype=np.float32))
+
+
+def test_compute_concept_attributions_orchestration_and_targets():
+    craft = _make_spatial_craft(number_of_concepts=3)
+    images = np.ones((2, 4, 4, 3), dtype=np.float32)
+
+    records = {"batch_sizes": [], "inputs": [], "targets": []}
+
+    class RecordingExplainer:
+        def __init__(self, model, batch_size):
+            self.model = model
+            records["batch_sizes"].append(batch_size)
+
+        def explain(self, inputs, targets):
+            concept_id = int(np.argmax(targets[0]))
+            records["inputs"].append(inputs.copy())
+            records["targets"].append(targets.copy())
+            scores = self.model(inputs)
+            assert scores.shape == (inputs.shape[0], craft.number_of_concepts)
+            return np.full(inputs.shape[:3], fill_value=float(concept_id), dtype=np.float32)
+
+    maps = craft.compute_concept_attributions(
+        images,
+        partial_explainer=PartialExplainer(RecordingExplainer),
+        concept_ids=[2, 0],
+    )
+
+    assert records["batch_sizes"] == [craft.batch_size]
+    assert len(records["inputs"]) == 2
+    np.testing.assert_array_equal(records["inputs"][0], images)
+    np.testing.assert_array_equal(records["inputs"][1], images)
+    assert len(records["targets"]) == 2
+    np.testing.assert_array_equal(records["targets"][0], [[0.0, 0.0, 1.0], [0.0, 0.0, 1.0]])
+    np.testing.assert_array_equal(records["targets"][1], [[1.0, 0.0, 0.0], [1.0, 0.0, 0.0]])
+
+    assert maps.shape == (2, 4, 4, 3)
+    assert np.all(np.isnan(maps[..., 1]))
+    assert np.all(np.isfinite(maps[..., 0]))
+    assert np.all(np.isfinite(maps[..., 2]))
+    np.testing.assert_allclose(maps[..., 2], 2.0)
+    np.testing.assert_allclose(maps[..., 0], 0.0)
+
+
+def test_compute_concept_attributions_validates_inputs_and_explainer_type():
+    craft = _make_spatial_craft(number_of_concepts=3)
+    images = np.ones((2, 4, 4, 3), dtype=np.float32)
+
+    class ShapeExplainer:
+        def __init__(self, model, batch_size):
+            del model, batch_size
+
+        def explain(self, inputs, targets):
+            del targets
+            return np.ones((inputs.shape[0], inputs.shape[1], inputs.shape[2], 2), dtype=np.float32)
+
+    with pytest.raises(TypeError, match="PartialExplainer"):
+        craft.compute_concept_attributions(images, partial_explainer=ShapeExplainer)
+
+    with pytest.raises(ValueError, match="between 0"):
+        craft.compute_concept_attributions(images, PartialExplainer(_Explainer), concept_ids=[3])
+    with pytest.raises(ValueError, match="duplicate"):
+        craft.compute_concept_attributions(images, PartialExplainer(_Explainer), concept_ids=[1, 1])
+    with pytest.raises(ValueError, match="at least one"):
+        craft.compute_concept_attributions(images, PartialExplainer(_Explainer), concept_ids=[])
+    with pytest.raises(ValueError, match="between 0"):
+        craft.compute_concept_attributions(images, PartialExplainer(_Explainer), concept_ids=[True])
+    with pytest.raises(ValueError, match="one channel"):
+        craft.compute_concept_attributions(images, PartialExplainer(ShapeExplainer))
+
+    class NonFiniteExplainer(ShapeExplainer):
+        def explain(self, inputs, targets):
+            del targets
+            return np.full(inputs.shape[:3], np.nan, dtype=np.float32)
+
+    with pytest.raises(ValueError, match="finite"):
+        craft.compute_concept_attributions(images, PartialExplainer(NonFiniteExplainer))
+
+
+@pytest.mark.parametrize(
+    "images",
+    [
+        np.empty((0, 4, 4, 3), dtype=np.float32),
+        np.empty((2, 4), dtype=np.float32),
+        np.empty((2, 0, 4, 3), dtype=np.float32),
+        [],
+    ],
+)
+def test_compute_concept_attributions_rejects_invalid_image_batches(images):
+    craft = _make_spatial_craft(number_of_concepts=3)
+
+    with pytest.raises(ValueError, match="images"):
+        craft.compute_concept_attributions(images, PartialExplainer(_Explainer), concept_ids=[0])
+
+
+def test_compute_concept_attributions_rejects_explicit_operators_before_setup():
+    craft = _make_spatial_craft(number_of_concepts=3)
+    images = np.ones((2, 4, 4, 3), dtype=np.float32)
+    partial_explainer = PartialExplainer(
+        _Explainer,
+        operator=lambda model, inputs, targets: model(inputs),
+    )
+
+    with pytest.raises(ValueError, match="does not accept a custom operator"):
+        craft.compute_concept_attributions(images, partial_explainer, concept_ids=[0])
+
+
+def test_compute_concept_attributions_accepts_explicit_none_operator():
+    craft = _make_spatial_craft(number_of_concepts=3)
+    images = np.ones((2, 4, 4, 3), dtype=np.float32)
+
+    class SingleChannelExplainer(_Explainer):
+        def __init__(self, model, batch_size, operator=None):
+            super().__init__(model, batch_size)
+            assert operator is None
+
+        def explain(self, coeffs_u, targets):
+            del targets
+            return np.ones(coeffs_u.shape[:3] + (1,), dtype=np.float32)
+
+    maps = craft.compute_concept_attributions(
+        images,
+        PartialExplainer(SingleChannelExplainer, operator=None),
+        concept_ids=[0],
+    )
+
+    assert np.all(np.isfinite(maps[..., 0]))
+
+
+def test_compute_concept_attributions_all_concepts_are_finite_by_default():
+    craft = _make_spatial_craft(number_of_concepts=3)
+    images = np.ones((2, 4, 4, 3), dtype=np.float32)
+
+    class ConstantExplainer:
+        def __init__(self, model, batch_size):
+            del model, batch_size
+
+        def explain(self, inputs, targets):
+            concept_id = int(np.argmax(targets[0]))
+            return np.full(inputs.shape[:3], concept_id + 1.0, dtype=np.float32)
+
+    maps = craft.compute_concept_attributions(images, PartialExplainer(ConstantExplainer))
+    assert maps.shape == (2, 4, 4, 3)
+    assert np.all(np.isfinite(maps))
+    np.testing.assert_allclose(maps[..., 0], 1.0)
+    np.testing.assert_allclose(maps[..., 1], 2.0)
+    np.testing.assert_allclose(maps[..., 2], 3.0)
+
+
+def test_compute_concept_attributions_rejects_whitebox_explainers_and_non_inductive_factorizer():
+    craft = _make_spatial_craft(number_of_concepts=3)
+    images = np.ones((2, 4, 4, 3), dtype=np.float32)
+
+    with pytest.raises(ValueError, match="black-box attribution"):
+        craft.compute_concept_attributions(images, PartialExplainer(GradientInput))
+
+    craft_non_inductive = _make_spatial_craft(
+        number_of_concepts=3,
+        factorizer=_NonInductiveFactorizer(),
+    )
+
+    class PassThroughExplainer:
+        def __init__(self, model, batch_size):
+            self.model = model
+
+        def explain(self, inputs, targets):
+            del targets
+            self.model(inputs)
+            return np.ones(inputs.shape[:3], dtype=np.float32)
+
+    with pytest.raises(RuntimeError, match="cannot encode unseen activations"):
+        craft_non_inductive.compute_concept_attributions(
+            images,
+            PartialExplainer(PassThroughExplainer),
+            concept_ids=[0],
+        )
+
+    class ExplainerNotImplemented:
+        def __init__(self, model, batch_size):
+            del model, batch_size
+
+        def explain(self, inputs, targets):
+            del inputs, targets
+            raise NotImplementedError("explainer operation is unavailable")
+
+    with pytest.raises(NotImplementedError, match="explainer operation"):
+        craft.compute_concept_attributions(
+            images,
+            PartialExplainer(ExplainerNotImplemented),
+            concept_ids=[0],
+        )
+
+
+def test_display_accepts_concept_maps_and_preserves_ranking_behavior():
+    craft = _make_spatial_craft(number_of_concepts=3)
+    images = np.ones((2, 4, 4, 3), dtype=np.float32)
+    coeffs_u = np.zeros((2, 2, 2, 3), dtype=np.float32)
+    coeffs_u[0, :, :, 0] = 10.0
+    coeffs_u[1, :, :, 0] = 1.0
+
+    concept_maps = np.full((2, 4, 4, 3), np.nan, dtype=np.float32)
+    concept_maps[..., 0] = 0.0
+    concept_maps[..., 2] = 2.0
+
+    displayed_maps = []
+
+    def fake_display_concept_heatmap(image, concept_heatmap, concept_idx, ax, **kwargs):
+        del ax, kwargs
+        displayed_maps.append((image.copy(), concept_idx, np.array(concept_heatmap)))
+
+    craft.display_concept_heatmap = fake_display_concept_heatmap
+
+    legacy_figure = craft.display_images_per_concept(images, coeffs_u=coeffs_u, order=[0])
+    assert [entry[1] for entry in displayed_maps] == [0, 0]
+    np.testing.assert_allclose([entry[2].max() for entry in displayed_maps], [10.0, 1.0])
+    plt.close(legacy_figure)
+    displayed_maps.clear()
+
+    def fail_if_transform_called(_):
+        raise AssertionError("map-only display must not recompute concept coefficients")
+
+    craft.transform = fail_if_transform_called
+    figure = craft.display_images_per_concept(images, coeffs_u=coeffs_u, concept_maps=concept_maps)
+    assert len(figure.axes) == 4
+    assert [entry[1] for entry in displayed_maps] == [0, 0, 2, 2]
+    np.testing.assert_allclose([entry[2].max() for entry in displayed_maps], [0.0, 0.0, 2.0, 2.0])
+    plt.close(figure)
+
+    with pytest.raises(ValueError, match="not available"):
+        craft.display_images_per_concept(
+            images,
+            coeffs_u=coeffs_u,
+            concept_maps=concept_maps,
+            order=[1],
+        )
+
+    partially_invalid_maps = concept_maps.copy()
+    partially_invalid_maps[0, 0, 0, 0] = np.nan
+    with pytest.raises(ValueError, match="entirely finite or entirely NaN"):
+        craft.display_images_per_concept(
+            images,
+            concept_maps=partially_invalid_maps,
+            order=[0],
+        )
+
+    rank_calls = []
+
+    def recording_topk(captured_coeffs, topk):
+        rank_calls.append(np.array(captured_coeffs))
+        return np.array([[0], [0], [0]], dtype=int)
+
+    craft.get_topk_images_per_concept = recording_topk
+    figure = craft.display_top_images_per_concept(
+        images,
+        topk=1,
+        coeffs_u=coeffs_u,
+        concept_maps=concept_maps,
+        order=[0],
+    )
+    assert len(rank_calls) == 1
+    np.testing.assert_array_equal(rank_calls[0], coeffs_u)
+    assert displayed_maps[-1][1] == 0
+    np.testing.assert_allclose(displayed_maps[-1][2], 0.0)
+    plt.close(figure)
+
+
+def test_display_concept_heatmap_resizing_and_finite_validation(monkeypatch):
+    craft = _make_spatial_craft(number_of_concepts=3)
+    image = np.ones((4, 4, 3), dtype=np.float32)
+    figure, ax = plt.subplots(1, 1)
+
+    resize_calls = {"count": 0}
+
+    def recording_resize(*args, **kwargs):
+        resize_calls["count"] += 1
+        return np.ones((4, 4, 1), dtype=np.float32)
+
+    monkeypatch.setattr("xplique.concepts.holistic_craft.cv2.resize", recording_resize)
+
+    craft.display_concept_heatmap(image, np.ones((4, 4), dtype=np.float32), concept_idx=0, ax=ax)
+    assert resize_calls["count"] == 0
+
+    craft.display_concept_heatmap(image, np.ones((2, 2), dtype=np.float32), concept_idx=0, ax=ax)
+    assert resize_calls["count"] == 1
+
+    craft.display_concept_heatmap(
+        image,
+        np.ones((4, 4, 1), dtype=np.float32),
+        concept_idx=0,
+        ax=ax,
+    )
+
+    with pytest.raises(ValueError, match="only finite"):
+        craft.display_concept_heatmap(
+            image,
+            np.where(np.indices((4, 4))[0] == 0, np.nan, 1.0).astype(np.float32),
+            concept_idx=0,
+            ax=ax,
+        )
+    plt.close(figure)
+
+
+def test_display_concept_heatmap_uses_absolute_magnitude(monkeypatch):
+    craft = _make_spatial_craft(number_of_concepts=3)
+    image = np.ones((2, 2, 3), dtype=np.float32)
+    heatmap = np.array([[-10.0, 1.0], [2.0, 3.0]], dtype=np.float32)
+    displayed = []
+
+    monkeypatch.setattr(
+        "xplique.concepts.holistic_craft.show_ax",
+        lambda img, ax, **kwargs: displayed.append(np.asarray(img)),
+    )
+    monkeypatch.setattr(
+        "xplique.concepts.holistic_craft._clip_percentile",
+        lambda values, percentile: values,
+    )
+
+    figure, ax = plt.subplots(1, 1)
+    craft.display_concept_heatmap(
+        image,
+        heatmap,
+        concept_idx=0,
+        ax=ax,
+        filter_percentile=75,
+    )
+
+    overlay = displayed[1][..., 0]
+    assert overlay[0, 0] == 10.0
+    assert np.count_nonzero(overlay) == 1
+    plt.close(figure)
+
+    monkeypatch.setattr(
+        "xplique.concepts.holistic_craft.cv2.resize",
+        lambda *args, **kwargs: np.full((4, 4), -1.0, dtype=np.float32),
+    )
+    displayed.clear()
+    figure, ax = plt.subplots(1, 1)
+    craft.display_concept_heatmap(
+        np.ones((4, 4, 3), dtype=np.float32),
+        heatmap,
+        concept_idx=0,
+        ax=ax,
+        filter_percentile=75,
+    )
+    assert np.all(displayed[1] >= 0.0)
+    plt.close(figure)
+
+
+def test_semantic_localization_follows_spatial_concept_dependencies():
+    craft = _Craft(
+        latent_data=None,
+        number_of_concepts=2,
+        extractor=_SemanticExtractor(),
+    )
+    images = np.ones((1, 4, 4, 2), dtype=np.float32)
+
+    def fail_if_differentiable_encoding_called(_):
+        raise AssertionError("black-box localization must not use encode_differentiable")
+
+    craft.factorizer.encode_differentiable = fail_if_differentiable_encoding_called
+
+    class PixelOcclusionExplainer:
+        def __init__(self, model, batch_size):
+            del batch_size
+            self.model = model
+
+        def explain(self, inputs, targets):
+            concept_id = int(np.argmax(targets[0]))
+            base_scores = self.model(inputs)
+            maps = np.zeros(inputs.shape[:3], dtype=np.float32)
+            for row in range(inputs.shape[1]):
+                for column in range(inputs.shape[2]):
+                    perturbed = inputs.copy()
+                    perturbed[:, row, column, :] = 0.0
+                    perturbed_scores = self.model(perturbed)
+                    maps[:, row, column] = (
+                        base_scores[:, concept_id] - perturbed_scores[:, concept_id]
+                    )
+            return maps
+
+    maps = craft.compute_concept_attributions(
+        images,
+        partial_explainer=PartialExplainer(PixelOcclusionExplainer),
+    )
+
+    assert maps.shape == (1, 4, 4, 2)
+    assert np.all(np.isfinite(maps))
+
+    concept_0_inside = maps[0, :2, :2, 0].sum()
+    concept_0_outside = maps[0, :, :, 0].sum() - concept_0_inside
+    concept_1_inside = maps[0, 2:, 2:, 1].sum()
+    concept_1_outside = maps[0, :, :, 1].sum() - concept_1_inside
+    assert concept_0_inside > concept_0_outside
+    assert concept_1_inside > concept_1_outside
diff --git a/xplique/concepts/holistic_craft.py b/xplique/concepts/holistic_craft.py
index 06090c0c..a42fdf2d 100644
--- a/xplique/concepts/holistic_craft.py
+++ b/xplique/concepts/holistic_craft.py
@@ -12,6 +12,7 @@
 from matplotlib.figure import Figure
 from sklearn.exceptions import NotFittedError
 
+from xplique.attributions.base import WhiteBoxExplainer
 from xplique.attributions.global_sensitivity_analysis.sobol_attribution_method import (
     SobolAttributionMethod,
 )
@@ -58,7 +59,7 @@ class PartialExplainer:
 
     This class stores an explainer class and its configuration kwargs, allowing
     the explainer to be instantiated later when the model and batch_size become
-    available during concept importance estimation.
+    available during concept importance estimation or concept localization.
 
     Parameters
     ----------
@@ -82,7 +83,7 @@ def __init__(self, explainer_class, **kwargs):
         if "model" in kwargs or "batch_size" in kwargs:
             raise ValueError(
                 "PartialExplainer should not receive 'model' or 'batch_size' arguments. "
-                "These will be provided automatically during importance estimation."
+                "These will be provided automatically during explanation."
             )
 
         self.explainer_class = explainer_class
@@ -107,6 +108,129 @@ def __call__(self, model, batch_size):
         return self.explainer_class(model=model, batch_size=batch_size, **self.kwargs)
 
 
+class ConceptLocalizer:
+    """Map input samples to one scalar score per learned concept.
+
+    Parameters
+    ----------
+    parent_craft
+        Fitted :class:`HolisticCraft` instance used to extract coefficients.
+    concept_reducer
+        Reduction applied over non-batch, non-concept coefficient dimensions.
+        Supported strings are ``"mean"``, ``"sum"``, and ``"max"``. A callable
+        may be supplied for custom reduction semantics. Callable reducers receive
+        the complete coefficient array with shape
+        ``(batch_size, *spatial_dimensions, number_of_concepts)`` and must return
+        an array with shape ``(batch_size, number_of_concepts)``. Coefficients are
+        reduced as-is, so signed factorizers retain their sign; use a callable
+        reducer such as a mean absolute value when magnitude is desired.
+
+    Returns
+    -------
+    concept_scores
+        Float32 concept scores with shape
+        ``(batch_size, number_of_concepts)``.
+
+    Notes
+    -----
+    The parent CRAFT instance must be fitted. Localization repeatedly evaluates
+    ``factorizer.encode()`` on perturbed inputs, so the factorizer must support
+    out-of-sample encoding.
+
+    Raises
+    ------
+    ValueError
+        If the reducer, concept coefficients, or reduced scores are invalid.
+    RuntimeError
+        If the factorizer cannot encode unseen activations.
+    """
+
+    parent_craft: "HolisticCraft"
+    concept_reducer: Union[str, Callable]
+
+    def __init__(
+        self,
+        parent_craft: "HolisticCraft",
+        concept_reducer: Union[str, Callable] = "mean",
+    ):
+        self.parent_craft = parent_craft
+        self.concept_reducer = concept_reducer
+        self._validate_reducer()
+
+    def _validate_reducer(self) -> None:
+        if isinstance(self.concept_reducer, str):
+            allowed_reducers = {"mean", "sum", "max"}
+            if self.concept_reducer not in allowed_reducers:
+                raise ValueError(
+                    "concept_reducer must be one of {'mean', 'sum', 'max'} "
+                    "or a callable returning shape (batch_size, number_of_concepts)."
+                )
+        elif not callable(self.concept_reducer):
+            raise ValueError(
+                "concept_reducer must be one of {'mean', 'sum', 'max'} "
+                "or a callable returning shape (batch_size, number_of_concepts)."
+            )
+
+    def _reduce_coefficients(self, coeffs_u: np.ndarray) -> np.ndarray:
+        """Reduce concept coefficients to one scalar score per concept and sample."""
+        coeffs_u = np.asarray(coeffs_u)
+        if coeffs_u.ndim < 2:
+            raise ValueError(
+                "Concept coefficients must have at least 2 dimensions: "
+                "(batch_size, number_of_concepts)."
+            )
+        if any(size == 0 for size in coeffs_u.shape):
+            raise ValueError(
+                f"Concept coefficients must not contain empty dimensions, got {coeffs_u.shape}."
+            )
+        if coeffs_u.shape[-1] != self.parent_craft.number_of_concepts:
+            raise ValueError(
+                f"Concept coefficients contain {coeffs_u.shape[-1]} concepts, expected "
+                f"{self.parent_craft.number_of_concepts}."
+            )
+
+        spatial_axes = tuple(range(1, coeffs_u.ndim - 1))
+        if not spatial_axes:
+            scores = coeffs_u
+        elif isinstance(self.concept_reducer, str):
+            reducers = {
+                "mean": np.mean,
+                "sum": np.sum,
+                "max": np.max,
+            }
+            scores = reducers[self.concept_reducer](coeffs_u, axis=spatial_axes)
+        else:
+            scores = self.concept_reducer(coeffs_u)
+
+        scores = np.asarray(scores)
+        expected_shape = (coeffs_u.shape[0], self.parent_craft.number_of_concepts)
+        if scores.shape != expected_shape:
+            raise ValueError(
+                f"Reduced concept scores must have shape {expected_shape}, got {scores.shape}."
+            )
+
+        scores = scores.astype(np.float32, copy=False)
+        if not np.all(np.isfinite(scores)):
+            raise ValueError("Reduced concept scores must contain only finite values.")
+        return scores
+
+    def _compute_scores(self, inputs: Any) -> np.ndarray:
+        """Encode inputs and reduce their concept coefficients."""
+        try:
+            coeffs_u = self.parent_craft.transform(inputs)
+        except NotImplementedError as error:
+            raise RuntimeError(
+                "The selected factorizer cannot encode unseen activations. "
+                "Input-to-concept localization requires factorizer.encode() "
+                "because attribution methods evaluate perturbed inputs."
+            ) from error
+        return self._reduce_coefficients(coeffs_u)
+
+    def __call__(self, inputs: Any) -> np.ndarray:
+        """Return reduced concept scores for a batch of inputs."""
+        return self._compute_scores(inputs)
+
+
 class HolisticCraft(ABC):
     """
     Framework-agnostic CRAFT implementation for holistic model explanations.
@@ -125,8 +249,14 @@ class HolisticCraft(ABC):
     The workflow involves:
     1. Extracting latent activations from a computer vision model
     2. Factorizing activations into interpretable concepts using NMF
-    3. Computing concept importance using gradient-based attribution methods
-    4. Visualizing concepts as spatial heatmaps overlaid on input images
+    3. Computing concept importance by attributing task predictions to concepts
+    4. Visualizing coefficient activation maps overlaid on input images
+    5. Optionally localizing concept scores to input regions with black-box attribution
+
+    Concept activation, concept importance, and concept localization answer different
+    questions. Coefficients describe how strongly a concept is present, importance
+    describes how much it contributes to a task prediction, and localization describes
+    which input regions drive a selected concept score.
 
     Parameters
     ----------
@@ -649,6 +779,298 @@ def make_concept_decoder(self, latent_data: LatentData) -> Any:
         """
         raise NotImplementedError
 
+    def make_concept_localizer(
+        self,
+        concept_reducer: Union[str, Callable] = "mean",
+    ) -> ConceptLocalizer:
+        """Create a callable mapping input images to one score per concept.
+
+        Parameters
+        ----------
+        concept_reducer
+            Reduction used to turn each concept coefficient map into a scalar
+            concept score. Supported strings are ``"mean"``, ``"sum"``, and
+            ``"max"``. A callable must return shape
+            ``(batch_size, number_of_concepts)``.
+
+        Returns
+        -------
+        localizer
+            Callable suitable for Xplique black-box explainers, returning
+            float32 scores with shape ``(batch_size, number_of_concepts)``.
+
+        Raises
+        ------
+        ValueError
+            If ``concept_reducer`` is invalid.
+        """
+        return ConceptLocalizer(self, concept_reducer)
+
+    def _validate_concept_ids(
+        self,
+        concept_ids: Optional[List[int]],
+        parameter_name: str = "concept_ids",
+    ) -> List[int]:
+        """Validate concept id selections and return them as a list."""
+        if concept_ids is None:
+            return list(range(self.number_of_concepts))
+
+        try:
+            concept_ids = list(concept_ids)
+        except TypeError as error:
+            raise ValueError(f"{parameter_name} must be an iterable of concept IDs") from error
+
+        if not concept_ids:
+            raise ValueError(f"{parameter_name} must contain at least one concept ID")
+        if len(concept_ids) > self.number_of_concepts:
+            raise ValueError(f"{parameter_name} cannot contain more IDs than number_of_concepts")
+        if any(
+            not isinstance(concept_id, (int, np.integer))
+            or isinstance(concept_id, bool)
+            or not 0 <= concept_id < self.number_of_concepts
+            for concept_id in concept_ids
+        ):
+            raise ValueError(
+                f"{parameter_name} concept IDs must be integers between 0 and "
+                f"{self.number_of_concepts - 1}"
+            )
+        if len(set(concept_ids)) != len(concept_ids):
+            raise ValueError(f"{parameter_name} cannot contain duplicate concept IDs")
+
+        return [int(concept_id) for concept_id in concept_ids]
+
+    def _to_numpy_image_batch(self, images: Union[np.ndarray, List[Any]]) -> np.ndarray:
+        """Convert image inputs to a validated rank-4 NumPy batch."""
+        if isinstance(images, list):
+            image_batch = []
+            for img in images:
+                image = np.asarray(self._to_numpy(img))
+                if image.ndim == 4 and image.shape[0] == 1:
+                    image = image[0]
+                if image.ndim != 3:
+                    expected_shape = (
+                        "(C, H, W) or (1, C, H, W)"
+                        if self.framework == "torch"
+                        else "(H, W, C) or (1, H, W, C)"
+                    )
+                    raise ValueError(f"images list entries must have shape {expected_shape}.")
+                image_batch.append(image)
+            if not image_batch:
+                raise ValueError("images list must contain at least one image.")
+            images_np = np.stack(image_batch, axis=0)
+        else:
+            images_np = self._to_numpy(images)
+            images_np = np.asarray(images_np)
+
+        if images_np.ndim == 3:
+            images_np = np.expand_dims(images_np, axis=0)
+
+        if images_np.ndim != 4:
+            expected_shape = "(N, C, H, W)" if self.framework == "torch" else "(N, H, W, C)"
+            raise ValueError(f"images must be a non-empty image batch with shape {expected_shape}.")
+        if images_np.shape[0] == 0 or min(images_np.shape[1:]) <= 0:
+            raise ValueError(
+                "images must contain at least one image with positive spatial dimensions."
+            )
+
+        return images_np
+
+    def _normalize_image_batch_to_nhwc(self, images: Union[np.ndarray, List[Any]]) -> np.ndarray:
+        """Convert a framework-native image batch to NHWC float32."""
+        images_np = self._to_numpy_image_batch(images)
+
+        if self.framework == "torch":
+            images_np = np.moveaxis(images_np, 1, -1)
+
+        return images_np.astype(np.float32, copy=False)
+
+    def compute_concept_attributions(
+        self,
+        images,
+        partial_explainer: PartialExplainer,
+        concept_ids: Optional[List[int]] = None,
+        concept_reducer: Union[str, Callable] = "mean",
+    ) -> np.ndarray:
+        """Compute input-space localization maps for selected learned concepts.
+
+        The fitted encoder and factorizer are exposed as a multi-output callable
+        whose output is one scalar score per concept. A standard black-box
+        attribution method then attributes each selected score to the input.
+
+        Parameters
+        ----------
+        images
+            Non-empty batch of images. TensorFlow inputs use ``(N, H, W, C)``;
+            PyTorch inputs use ``(N, C, H, W)``. A single image or a list of
+            images in the corresponding framework layout is also accepted.
+        partial_explainer
+            Deferred black-box Xplique explainer configuration. Omit ``operator``
+            or leave it as ``None``; one-hot concept targets select the score directly.
+        concept_ids
+            Optional iterable of unique concept IDs to localize. Uncomputed
+            result channels are filled entirely with ``NaN``. If omitted, all
+            concepts are localized.
+        concept_reducer
+            Reduction from coefficient maps to scalar concept scores.
+        Returns
+        -------
+        concept_maps
+            Float32 maps with shape ``(N, H, W, number_of_concepts)``.
+            Channel ``k`` always corresponds to concept ``k``.
+
+        Notes
+        -----
+        Localization evaluates the fitted factorizer on perturbed inputs, so the
+        factorizer must support out-of-sample ``encode()``. The localizer keeps
+        signed concept scores unchanged; use a custom reducer when concept
+        magnitude rather than signed activation is intended. White-box explainers
+        are rejected when ``explainer_class`` is a class. Callable factories are
+        allowed and are responsible for producing a compatible black-box explainer.
+
+        Raises
+        ------
+        NotFittedError
+            If CRAFT has not been fitted.
+        TypeError
+            If ``partial_explainer`` is not a :class:`PartialExplainer`.
+        ValueError
+            If inputs, concept IDs, reducer, explainer type, or explainer output
+            shape is invalid.
+        RuntimeError
+            If the factorizer cannot encode perturbed, unseen activations.
+        """
+        self.check_if_fitted()
+
+        if not isinstance(partial_explainer, PartialExplainer):
+            raise TypeError(
+                f"partial_explainer must be a PartialExplainer instance, got "
+                f"{type(partial_explainer).__name__}."
+            )
+
+        if partial_explainer.kwargs.get("operator") is not None:
+            raise ValueError(
+                "Concept localization uses one-hot concept targets and does not accept a "
+                "custom operator. Omit the operator argument."
+            )
+
+        explainer_class = partial_explainer.explainer_class
+        # Class-based white-box explainers can be rejected explicitly. Callable factories
+        # are intentionally left to their own validation because their class is unknown.
+        if isinstance(explainer_class, type) and issubclass(explainer_class, WhiteBoxExplainer):
+            raise ValueError(
+                "Input-to-concept localization currently supports black-box attribution "
+                "methods only. The ConceptLocalizer uses factorizer.encode(), which is "
+                "not guaranteed to be differentiable. Use Rise, SobolAttributionMethod, "
+                "Occlusion, Lime, KernelShap, or another compatible black-box explainer."
+            )
+
+        selected_concepts = self._validate_concept_ids(concept_ids, parameter_name="concept_ids")
+        attribution_inputs = self._normalize_image_batch_to_nhwc(images)
+        localizer = self.make_concept_localizer(concept_reducer)
+
+        explainer_instance = partial_explainer(model=localizer, batch_size=self.batch_size)
+
+        num_images, height, width, _ = attribution_inputs.shape
+        concept_maps = np.full(
+            (num_images, height, width, self.number_of_concepts),
+            np.nan,
+            dtype=np.float32,
+        )
+
+        for concept_id in selected_concepts:
+            targets = np.zeros((num_images, self.number_of_concepts), dtype=np.float32)
+            targets[:, concept_id] = 1.0
+
+            single_concept_map = explainer_instance.explain(attribution_inputs, targets)
+            single_concept_map = self._to_numpy(single_concept_map)
+            if single_concept_map.ndim == 4:
+                if single_concept_map.shape[-1] != 1:
+                    raise ValueError(
+                        "Concept attribution maps must have one channel when returned as "
+                        f"rank-4 tensors, got {single_concept_map.shape}."
+                    )
+                single_concept_map = single_concept_map[..., 0]
+            elif single_concept_map.ndim != 3:
+                raise ValueError(
+                    "Concept attribution maps must have shape (N, H, W) or (N, H, W, 1), "
+                    f"got {single_concept_map.shape}."
+                )
+
+            expected_shape = (num_images, height, width)
+            if single_concept_map.shape != expected_shape:
+                raise ValueError(
+                    "Concept attribution map shape must match attribution inputs spatial "
+                    f"shape {expected_shape}, got {single_concept_map.shape}."
+                )
+            if not np.all(np.isfinite(single_concept_map)):
+                raise ValueError(
+                    f"Concept attribution map for concept {concept_id} must contain "
+                    "only finite values."
+                )
+
+            concept_maps[..., concept_id] = single_concept_map.astype(np.float32, copy=False)
+
+        return concept_maps
+
+    def _prepare_concept_maps(
+        self,
+        concept_maps: np.ndarray,
+    ) -> Tuple[np.ndarray, List[int]]:
+        """Validate concept maps and identify which concept channels are available."""
+        concept_maps = np.asarray(concept_maps, dtype=np.float32)
+
+        if concept_maps.ndim != 4:
+            raise ValueError("concept_maps must have shape (N, H, W, n_concepts).")
+        if any(size == 0 for size in concept_maps.shape):
+            raise ValueError(
+                f"concept_maps must not contain empty dimensions, got {concept_maps.shape}."
+            )
+        if concept_maps.shape[-1] != self.number_of_concepts:
+            raise ValueError(
+                f"concept_maps contains {concept_maps.shape[-1]} concepts, expected "
+                f"{self.number_of_concepts}."
+            )
+
+        channel_is_uncomputed = np.all(np.isnan(concept_maps), axis=(0, 1, 2))
+        channel_is_computed = np.all(np.isfinite(concept_maps), axis=(0, 1, 2))
+        if not np.all(channel_is_uncomputed | channel_is_computed):
+            raise ValueError("Each concept map channel must be entirely finite or entirely NaN.")
+
+        available_concepts = np.flatnonzero(channel_is_computed).tolist()
+        return concept_maps, available_concepts
+
+    def _resolve_concept_map_source(
+        self,
+        concept_maps: np.ndarray,
+        images_np: np.ndarray,
+        requested_concepts: Optional[List[int]],
+    ) -> Tuple[np.ndarray, List[int]]:
+        """Validate display maps and resolve the concept columns to display."""
+        concept_maps, available_concepts = self._prepare_concept_maps(concept_maps)
+        if concept_maps.shape[0] != images_np.shape[0]:
+            raise ValueError(
+                "concept_maps and images must contain the same number of samples, "
+                f"got {concept_maps.shape[0]} and {images_np.shape[0]}."
+            )
+
+        if requested_concepts is None:
+            if not available_concepts:
+                raise ValueError("No computed concept maps are available to display.")
+            concepts_id = available_concepts
+        else:
+            unavailable = [
+                concept_id
+                for concept_id in requested_concepts
+                if concept_id not in available_concepts
+            ]
+            if unavailable:
+                raise ValueError(
+                    f"Requested concept maps are not available for concept IDs: {unavailable}."
+                )
+            concepts_id = requested_concepts
+
+        return concept_maps, concepts_id
+
     def _prepare_display_concept_inputs(
         self,
         images: Union[np.ndarray, List[Any]],
@@ -705,36 +1127,8 @@ def _prepare_display_concept_inputs(
                 f"{self.number_of_concepts}"
             )
 
-        # convert images to HWC numpy format for display
-        if self.framework == "torch":
-            # channel first (C, H, W) -> channel last (H, W, C) for each image
-            images_np = np.stack([self._to_numpy(img.squeeze().permute(1, 2, 0)) for img in images])
-        else:
-            images_np = np.stack([self._to_numpy(img) for img in images])
-
-        if order is None:
-            concepts_id = list(range(self.number_of_concepts))
-        else:
-            try:
-                concepts_id = list(order)
-            except TypeError as error:
-                raise ValueError("order must be an iterable of concept IDs") from error
-            if not concepts_id:
-                raise ValueError("order must contain at least one concept ID")
-            if len(concepts_id) > self.number_of_concepts:
-                raise ValueError("order cannot contain more IDs than number_of_concepts")
-            if any(
-                not isinstance(concept_id, (int, np.integer))
-                or isinstance(concept_id, bool)
-                or not 0 <= concept_id < self.number_of_concepts
-                for concept_id in concepts_id
-            ):
-                raise ValueError(
-                    f"order concept IDs must be integers between 0 and "
-                    f"{self.number_of_concepts - 1}"
-                )
-            if len(set(concepts_id)) != len(concepts_id):
-                raise ValueError("order cannot contain duplicate concept IDs")
+        images_np = self._normalize_image_batch_to_nhwc(images)
+        concepts_id = self._validate_concept_ids(order, parameter_name="order")
 
         return images_np, coeffs_u, concepts_id
 
@@ -749,16 +1143,18 @@ def display_concept_heatmap(
     ) -> None:
         """Overlay a single concept heatmap on a single image.
 
-        Displays the image on the given axis, then overlays the concept activation
-        heatmap after filtering low activations, resizing to image resolution, and
-        clipping outlier values.
+        Displays the image on the given axis, then overlays the concept heatmap
+        after filtering low-magnitude activations and clipping outlier values.
+        The heatmap may be either a latent concept activation map or an
+        input-space attribution map. Absolute magnitudes are displayed, so signed
+        attribution direction is intentionally not represented.
 
         Parameters
         ----------
         image
             Single image as HWC numpy array, shape (H, W, C)
         concept_heatmap
-            Raw concept activation map, shape (H', W')
+            Concept heatmap, shape (H', W') or (H', W', 1)
         concept_idx
             Index of the concept, used to select the colormap
         ax
@@ -772,17 +1168,33 @@ def display_concept_heatmap(
             This parameter allows to avoid outliers in case of too extreme values.
             Default to 5.
         """
-        dsize = (image.shape[1], image.shape[0])  # cv2 expects (width, height)
+        concept_heatmap = np.asarray(concept_heatmap, dtype=np.float32)
+        if concept_heatmap.ndim == 3 and concept_heatmap.shape[-1] == 1:
+            concept_heatmap = concept_heatmap[..., 0]
+        elif concept_heatmap.ndim != 2:
+            raise ValueError(
+                f"concept_heatmap must have shape (H, W) or (H, W, 1), got {concept_heatmap.shape}."
+            )
+
+        if not np.all(np.isfinite(concept_heatmap)):
+            raise ValueError("concept_heatmap must contain only finite values.")
+
+        concept_heatmap = np.abs(concept_heatmap)
 
         # Display the image
         show_ax(image, ax=ax)
 
         # only show concept if excess N-th percentile
-        sigma = np.percentile(concept_heatmap.flatten(), filter_percentile)
+        sigma = np.percentile(concept_heatmap, filter_percentile)
         heatmap = concept_heatmap * (concept_heatmap > sigma)
 
-        # resize the heatmap before clipping
-        heatmap = cv2.resize(heatmap[:, :, None], dsize=dsize, interpolation=cv2.INTER_CUBIC)
+        # resize the heatmap before clipping when needed
+        if heatmap.shape[:2] != image.shape[:2]:
+            dsize = (image.shape[1], image.shape[0])  # cv2 expects (width, height)
+            heatmap = cv2.resize(heatmap[:, :, None], dsize=dsize, interpolation=cv2.INTER_CUBIC)
+            heatmap = np.maximum(heatmap, 0.0)
+        else:
+            heatmap = heatmap[:, :, None]
         heatmap = _clip_percentile(heatmap, clip_percentile)
 
         # Display the heatmap overlay
@@ -796,6 +1208,7 @@ def display_images_per_concept(
         filter_percentile: int = 80,
         clip_percentile: int = 5,
         order: Optional[List[int]] = None,
+        concept_maps: Optional[np.ndarray] = None,
     ) -> Figure:
         """
         Display concept heatmaps overlaid on images.
@@ -823,15 +1236,32 @@ def display_images_per_concept(
         order
             Optional list of concept IDs to specify display order. If None,
             concepts are shown in sequential order
+        concept_maps
+            Optional input-space concept attribution maps of shape
+            (N, H, W, n_concepts). When provided, maps are used for overlays
+            while `coeffs_u` keeps its original meaning. Each concept channel
+            must be either entirely finite or entirely ``NaN``. If ``order`` is
+            None, only finite channels are displayed. Signed maps are displayed
+            by absolute magnitude.
 
         Returns
         -------
         fig
             matplotlib figure with len(images) rows and number_of_concepts columns
         """
-        images_np, coeffs_u, concepts_id = self._prepare_display_concept_inputs(
-            images, coeffs_u, order
-        )
+        if concept_maps is None:
+            images_np, coeffs_u, concepts_id = self._prepare_display_concept_inputs(
+                images, coeffs_u, order
+            )
+            heatmap_source = coeffs_u
+        else:
+            images_np = self._normalize_image_batch_to_nhwc(images)
+            requested_concepts = (
+                None if order is None else self._validate_concept_ids(order, parameter_name="order")
+            )
+            heatmap_source, concepts_id = self._resolve_concept_map_source(
+                concept_maps, images_np, requested_concepts
+            )
 
         nb_cols = len(concepts_id)
         nb_rows = len(images_np)
@@ -847,7 +1277,7 @@ def display_images_per_concept(
             for image_id, image in enumerate(images_np):
                 self.display_concept_heatmap(
                     image=image,
-                    concept_heatmap=coeffs_u[image_id, :, :, c_i],
+                    concept_heatmap=heatmap_source[image_id, :, :, c_i],
                     concept_idx=c_i,
                     ax=axs[image_id, i],
                     filter_percentile=filter_percentile,
@@ -891,6 +1321,7 @@ def display_top_images_per_concept(
         clip_percentile: int = 5,
         order: Optional[List[int]] = None,
         coeffs_u: Optional[np.ndarray] = None,
+        concept_maps: Optional[np.ndarray] = None,
     ) -> Figure:
         """Display top N images per concept ranked by average activation.
 
@@ -910,6 +1341,13 @@ def display_top_images_per_concept(
             Optional pre-computed concept coefficients. If None, will call
             self.transform(images) to compute them. Use this to pass the
             coefficients stored in factorization.coeffs_u after fit().
+        concept_maps
+            Optional input-space concept attribution maps of shape
+            (N, H, W, n_concepts) used only for display overlays. Top-image
+            ranking remains based on mean concept coefficients. Each concept
+            channel must be either entirely finite or entirely ``NaN``. If
+            ``order`` is None, only finite channels are displayed. Signed maps
+            are displayed by absolute magnitude.
 
         Returns
         -------
@@ -920,6 +1358,13 @@ def display_top_images_per_concept(
             images, coeffs_u, order
         )
 
+        heatmap_source = coeffs_u
+        if concept_maps is not None:
+            requested_concepts = None if order is None else concepts_id
+            heatmap_source, concepts_id = self._resolve_concept_map_source(
+                concept_maps, images_np, requested_concepts
+            )
+
         nb_rows = topk
         nb_cols = len(concepts_id)
         fig, axs = plt.subplots(nb_rows, nb_cols, figsize=(2 * nb_cols, 2 * nb_rows))
@@ -935,7 +1380,7 @@ def display_top_images_per_concept(
             for j, image_id in enumerate(topk_images_ids[c_i]):
                 self.display_concept_heatmap(
                     image=images_np[image_id],
-                    concept_heatmap=coeffs_u[image_id, :, :, c_i],
+                    concept_heatmap=heatmap_source[image_id, :, :, c_i],
                     concept_idx=c_i,
                     ax=axs[j, i],
                     filter_percentile=filter_percentile,
diff --git a/xplique/concepts/tf/holistic_craft.py b/xplique/concepts/tf/holistic_craft.py
index b3b41521..b893879b 100644
--- a/xplique/concepts/tf/holistic_craft.py
+++ b/xplique/concepts/tf/holistic_craft.py
@@ -1,6 +1,6 @@
 """TensorFlow-specific wrapper for HolisticCraft."""
 
-from typing import Any, Optional, Union
+from typing import Any, Callable, Optional, Union
 
 import numpy as np
 import tensorflow as tf
@@ -10,7 +10,7 @@
     _pad_and_stack_box_predictions,
 )
 
-from ..holistic_craft import ConceptDecoder, HolisticCraft
+from ..holistic_craft import ConceptDecoder, ConceptLocalizer, HolisticCraft
 from ..latent_extractor import LatentData
 from .factorizer import TfSklearnNMFFactorizer
 from .latent_extractor import TfLatentExtractor as LatentExtractor
@@ -173,6 +173,33 @@ def make_concept_decoder(self, latent_data: LatentData) -> tf.keras.layers.Layer
 
         return ConceptDecoderTf(self, latent_data)
 
+    def make_concept_localizer(
+        self,
+        concept_reducer: Union[str, Callable] = "mean",
+    ) -> ConceptLocalizer:
+        """Create a TensorFlow concept localizer for black-box attribution.
+
+        Parameters
+        ----------
+        concept_reducer
+            Reduction from coefficient maps to one scalar score per concept.
+
+        Returns
+        -------
+        localizer
+            Callable returning a tensor with shape ``(batch_size, K)``.
+        """
+        return ConceptLocalizerTf(self, concept_reducer)
+
+
+class ConceptLocalizerTf(ConceptLocalizer):
+    """TensorFlow concept localizer callable."""
+
+    def __call__(self, inputs: Any) -> tf.Tensor:
+        """Return reduced concept scores for a batch of inputs."""
+        scores = self._compute_scores(inputs)
+        return tf.convert_to_tensor(scores, dtype=tf.float32)
+
 
 class ConceptDecoderTf(tf.keras.layers.Layer, ConceptDecoder):
     """
diff --git a/xplique/concepts/torch/holistic_craft.py b/xplique/concepts/torch/holistic_craft.py
index 9c718d1a..dfd57411 100644
--- a/xplique/concepts/torch/holistic_craft.py
+++ b/xplique/concepts/torch/holistic_craft.py
@@ -1,6 +1,6 @@
 """PyTorch-specific wrapper for HolisticCraft."""
 
-from typing import Any, Optional, Union
+from typing import Any, Callable, Optional, Union
 
 import numpy as np
 import torch
@@ -12,7 +12,7 @@
 )
 from xplique.wrappers import TorchWrapper
 
-from ..holistic_craft import ConceptDecoder, HolisticCraft
+from ..holistic_craft import ConceptDecoder, ConceptLocalizer, HolisticCraft
 from ..latent_extractor import LatentData
 from .factorizer import TorchSklearnNMFFactorizer
 from .latent_extractor import TorchLatentExtractor as LatentExtractor
@@ -190,6 +190,48 @@ def make_concept_decoder(self, latent_data: LatentData) -> TorchWrapper:
         )
         return wrapped_decoder
 
+    def make_concept_localizer(
+        self,
+        concept_reducer: Union[str, Callable] = "mean",
+    ) -> TorchWrapper:
+        """Create a PyTorch concept localizer for black-box attribution.
+
+        Parameters
+        ----------
+        concept_reducer
+            Reduction from coefficient maps to one scalar score per concept.
+
+        Returns
+        -------
+        localizer
+            Xplique ``TorchWrapper`` returning a tensor with shape
+            ``(batch_size, K)`` and gradients disabled.
+        """
+        torch_localizer = ConceptLocalizerTorch(self, concept_reducer).eval()
+        return TorchWrapper(
+            torch_localizer,
+            device=self.device,
+            is_channel_first=True,
+            requires_grad=False,
+        )
+
+
+class ConceptLocalizerTorch(nn.Module, ConceptLocalizer):
+    """PyTorch concept localizer module."""
+
+    def __init__(
+        self,
+        parent_craft: HolisticCraft,
+        concept_reducer: Union[str, Callable] = "mean",
+    ) -> None:
+        super().__init__()
+        ConceptLocalizer.__init__(self, parent_craft, concept_reducer)
+
+    def forward(self, inputs: torch.Tensor) -> torch.Tensor:
+        """Return reduced concept scores for native NCHW inputs."""
+        scores = self._compute_scores(inputs)
+        return torch.as_tensor(scores, dtype=torch.float32, device=inputs.device)
+
 
 class ConceptDecoderTorch(nn.Module, ConceptDecoder):
     """