diff --git a/doc/source/_workflows/accuracy_config.yaml b/doc/source/_workflows/accuracy_config.yaml
index 39114b7f..7fe837f9 100644
--- a/doc/source/_workflows/accuracy_config.yaml
+++ b/doc/source/_workflows/accuracy_config.yaml
@@ -8,8 +8,8 @@ outputs:
   path: "outputs_accuracy"
 coregistration:
   step_one:
-    method: "LZD"
-    extra_information: {"subsample": 10000}
+    method: NuthKaab
+    extra_information: null
 statistics:
 - median
 - nmad
diff --git a/doc/source/cli.md b/doc/source/cli.md
index fe01ac75..41c359b1 100644
--- a/doc/source/cli.md
+++ b/doc/source/cli.md
@@ -37,6 +37,8 @@ Also, this template can directly be saved in a YAML file like this:
 xdem workflow_name --template-config template_config.yaml
 ```
 
+Optionally, a output path to save the template can be set directly in the command-line using `--output`.
+
 When edited by a user, a configuration file **must contain at minima the input parameters listed as "required"** on the documentation page of the given workflow.
 xDEM then automatically fills in the rest with default settings. Users are free to edit the configuration file to run only the parts they need.
 
diff --git a/doc/source/cli_accuracy.md b/doc/source/cli_accuracy.md
index 3e260489..95192653 100644
--- a/doc/source/cli_accuracy.md
+++ b/doc/source/cli_accuracy.md
@@ -122,7 +122,7 @@ Elevation input information, split between reference and to-be-aligned elevation
 ::::{tab-item} `reference_elev`
 
 :::{table} Inputs parameters for `reference_elev`
-:widths: 20, 35, 17, 18, 10
+:widths: 20, 40, 20, 10, 10
 
 | Name                  | Description                              | Type       | Default | Required |
 |-----------------------|------------------------------------------|------------|---------|----------|
diff --git a/xdem/cli.py b/xdem/cli.py
index 32c6c78b..93538938 100644
--- a/xdem/cli.py
+++ b/xdem/cli.py
@@ -106,8 +106,7 @@ def main() -> None:
         epilog="examples:\n"
         "  xdem accuracy --config config.yaml\n"
         "  xdem accuracy --config config.yaml --output myoutputfolder\n"
-        "  xdem accuracy --template-config\n"
-        "  xdem accuracy --template-config template_config.yaml",
+        "  xdem accuracy --template-config --output template_config.yaml",
         add_help=False,
         formatter_class=argparse.RawTextHelpFormatter,
     )
diff --git a/xdem/workflows/accuracy.py b/xdem/workflows/accuracy.py
index 58eab3a2..325a3ba1 100644
--- a/xdem/workflows/accuracy.py
+++ b/xdem/workflows/accuracy.py
@@ -74,17 +74,26 @@ def __init__(self, config_dem: str | Dict[str, Any], output: str | None = None)
 
         self.config = self.remove_none(self.config)  # type: ignore
 
-    def _load_data(self) -> None:
-        """Load data."""
+    def _load_data(self) -> tuple[float, float]:
+        """
+        Load data
 
+        :return vmin, vmax: to plot elevation data with the same scale
+        """
         self.to_be_aligned_elev, tba_mask, tba_path_mask = self.load_dem(self.config["inputs"]["to_be_aligned_elev"])
         self.reference_elev, ref_mask, ref_mask_path = self.load_dem(self.config["inputs"].get("reference_elev", None))
         if self.reference_elev is None:
             self.reference_elev = self._get_reference_elevation()
+
+        vmin = min(self.reference_elev.get_stats("min"), self.to_be_aligned_elev.get_stats("min"))
+        vmax = max(self.reference_elev.get_stats("max"), self.to_be_aligned_elev.get_stats("max"))
+
         self.generate_plot(
             self.reference_elev,
             title="Reference DEM",
             filename="reference_elev_map",
+            vmin=vmin,
+            vmax=vmax,
             cmap="terrain",
             cbar_title="Elevation (m)",
         )
@@ -92,6 +101,8 @@ def _load_data(self) -> None:
             self.to_be_aligned_elev,
             title="To-be-aligned DEM",
             filename="to_be_aligned_elev_map",
+            vmin=vmin,
+            vmax=vmax,
             cmap="terrain",
             cbar_title="Elevation (m)",
         )
@@ -109,11 +120,15 @@ def _load_data(self) -> None:
                 self.to_be_aligned_elev,
                 title="Masked (inlier) terrain",
                 filename="masked_elev_map",
+                vmin=vmin,
+                vmax=vmax,
                 mask_path=path_mask,
                 cmap="terrain",
                 cbar_title="Elevation (m)",
             )
 
+        return vmin, vmax
+
     def _get_reference_elevation(self) -> float:
         """
         Get reference elevation.
@@ -169,9 +184,12 @@ def _compute_coregistration(self) -> RasterType:
 
         return aligned_elev
 
-    def _prepare_datas(self) -> None:
+    def _prepare_datas(self, vmin: float, vmax: float) -> None:
         """
         Compute reprojection.
+
+        :param vmin: to plot elevation data with the same scale
+        :param vmax: to plot elevation data with the same scale
         """
         sampling_source = self.config["inputs"]["sampling_grid"]
 
@@ -195,21 +213,32 @@ def _prepare_datas(self) -> None:
         # Intersection
         logging.info("Computing intersection")
         coord_intersection = self.reference_elev.intersection(self.to_be_aligned_elev)
+        self.to_be_aligned_elev = self.to_be_aligned_elev.crop(coord_intersection)
+        self.reference_elev = self.reference_elev.crop(coord_intersection)
+        coord_intersection = self.reference_elev.intersection(self.to_be_aligned_elev)
+        self.reference_elev = self.reference_elev.crop(coord_intersection)
+
+        coord_intersection = self.to_be_aligned_elev.intersection(self.reference_elev)
+
         if sampling_source == "reference_elev":
-            self.reference_elev = self.reference_elev.crop(coord_intersection)
+            self.to_be_aligned_elev = self.to_be_aligned_elev.crop(coord_intersection)
             self.generate_plot(
                 self.to_be_aligned_elev,
-                title="Cropped reference DEM",
-                filename="cropped_reference_elev_map",
+                title="Preprocessed to-be-aligned DEM",
+                filename="preprocessed_to_be_aligned_elev_map",
+                vmin=vmin,
+                vmax=vmax,
                 cmap="terrain",
                 cbar_title="Elevation (m)",
             )
         else:
-            self.to_be_aligned_elev = self.to_be_aligned_elev.crop(coord_intersection)
+            self.reference_elev = self.reference_elev.crop(coord_intersection)
             self.generate_plot(
-                self.to_be_aligned_elev,
-                title="Cropped to-be-aligned DEM",
-                filename="cropped_to_be_aligned_elev_map",
+                self.reference_elev,
+                title="Preprocessed reference DEM",
+                filename="preprocessed_reference_elev_map",
+                vmin=vmin,
+                vmax=vmax,
                 cmap="terrain",
                 cbar_title="Elevation (m)",
             )
@@ -305,11 +334,11 @@ def run(self) -> None:
 
         t0 = time.time()
 
-        self._load_data()
+        vmin, vmax = self._load_data()
 
         # Reprojection step
         if "sampling_grid" in self.config["inputs"]:
-            self._prepare_datas()
+            self._prepare_datas(vmin, vmax)
 
         if self.compute_coreg:
             # Coregistration step
@@ -321,31 +350,13 @@ def run(self) -> None:
         output_grid = self.config["outputs"]["output_grid"]
         ref_elev = self.reference_elev if output_grid == "reference_elev" else self.to_be_aligned_elev
 
-        vmin = vmax = None
-
-        if self.compute_coreg:
-            diff_pairs = [("before", self.to_be_aligned_elev), ("after", aligned_elev.reproject(ref_elev))]
-        else:
-            diff_pairs = [("", self.to_be_aligned_elev)]
-
-        for label, dem in diff_pairs:
-            diff = dem - ref_elev
-            stats_keys = ["min", "max", "nmad", "median"]
-            stats = diff.get_stats(stats_keys)
-
-            if label == "before":
-                self.diff_before, self.stats_before = diff, stats
-                vmin, vmax = -(stats["median"] + 3 * stats["nmad"]), stats["median"] + 3 * stats["nmad"]
-            elif label == "after":
-                self.diff_after, self.stats_after = diff, stats
-            else:
-                self.diff = diff
-                vmin, vmax = -(stats["median"] + 3 * stats["nmad"]), (stats["median"] + 3 * stats["nmad"])
+        stats_keys = ["min", "max", "nmad", "median"]
 
+        def generate_plot_diff(label: str, diff: RasterType, vmin: float, vmax: float) -> None:
             suffix = f"_elev_{label}_coreg_map" if label else "_elev"
             self.generate_plot(
                 diff,
-                title=f"Difference\n{label} coregistration",
+                title=f"Difference {label} coregistration",
                 filename=f"diff{suffix}",
                 vmin=vmin,
                 vmax=vmax,
@@ -353,6 +364,32 @@ def run(self) -> None:
                 cbar_title="Elevation differences (m)",
             )
 
+        if self.compute_coreg:
+
+            self.diff_before = self.to_be_aligned_elev - ref_elev
+            self.stats_before = self.diff_before.get_stats(stats_keys)
+
+            self.diff_after = aligned_elev.reproject(ref_elev) - ref_elev
+            self.stats_after = self.diff_after.get_stats(stats_keys)
+
+            vmin_diff = min(
+                -(self.stats_before["median"] + 3 * self.stats_before["nmad"]),
+                -(self.stats_after["median"] + 3 * self.stats_after["nmad"]),
+            )
+            vmax_diff = max(
+                self.stats_before["median"] + 3 * self.stats_before["nmad"],
+                self.stats_after["median"] + 3 * self.stats_after["nmad"],
+            )
+
+            generate_plot_diff("before", self.diff_before, vmin_diff, vmax_diff)
+            generate_plot_diff("after", self.diff_after, vmin_diff, vmax_diff)
+
+        else:
+            self.diff = self.to_be_aligned_elev - ref_elev
+            self.stats = self.diff.get_stats(stats_keys)
+            vmin, vmax = -(self.stats["median"] + 3 * self.stats["nmad"]), self.stats["median"] + 3 * self.stats["nmad"]
+            generate_plot_diff("", self.diff, vmin, vmax)
+
         if self.compute_coreg:
             stat_items = [
                 (self.reference_elev, "reference_elev", "Reference elevation", 2),
@@ -391,6 +428,7 @@ def run(self) -> None:
 
         if len(list_df_var) > 0:
             df_stats = pd.concat(list_df_var)
+            df_stats.set_index("Data", inplace=True)
         else:
             df_stats = None
         self.df_stats = df_stats
@@ -436,14 +474,14 @@ def create_html(self, list_dict: list[tuple[str, dict[str, Any]]]) -> None:
 
         # Plot input elevation data
         html += "<h2>Elevation datasets</h2>\n"
-        html += "<div style='display: flex; gap: 10px;'>\n"
+        html += "<div style='display: flex'>\n"
         html += (
-            "  <img src='plots/reference_elev_map.png' alt='Image PNG' "
-            "style='max-width: 50%; height: auto; width: 50%;'>\n"
+            "<img src='plots/reference_elev_map.png' alt='Image PNG' "
+            "style='max-width: 50%; height: auto; width: 50%;'>"
         )
         html += (
-            "  <img src='plots/to_be_aligned_elev_map.png' alt='Image PNG' style='max-width: "
-            "50%; height: auto; width: 50%;'>\n"
+            "<img src='plots/to_be_aligned_elev_map.png' alt='Image PNG' "
+            "style='max-width: 50%; height: auto; width: 50%;'>"
         )
         html += "</div>\n"
 
@@ -456,23 +494,55 @@ def format_values(val: Any) -> Any:
             else:
                 return str(val)
 
-        # Metadata: Inputs, coregistration
-        for title, dictionary in list_dict:  # type: ignore
-            html += "<div style='clear: both; margin-bottom: 30px;'>\n"
-            html += f"<h2>{title}</h2>\n"
-            html += "<table border='1' cellspacing='0' cellpadding='5'>\n"
-            html += "<tr><th>Information</th><th>Value</th></tr>\n"
+        def print_dict(title: str, dictionary: dict[str, Any]) -> str:
+            div_html = "<div style='clear: both; margin-bottom: 30px;'>\n"
+            div_html += f"<h2>{title}</h2>\n"
+            div_html += "<table border='1' cellspacing='0' cellpadding='5'>\n"
+            div_html += "<tr><th>Information</th><th>Value</th></tr>\n"
             for key, value in dictionary.items():
                 if isinstance(value, dict):
                     value = {k: format_values(v) for k, v in value.items()}
-                html += f"<tr><td>{key}</td><td>{value}</td></tr>\n"
-            html += "</table>\n"
+                div_html += f"<tr><td>{key}</td><td>{value}</td></tr>\n"
+            div_html += "</table>\n"
+            div_html += "</div>\n"
+            return div_html
+
+        # Metadata: Inputs
+        inputs_information = list_dict[0]
+        html += print_dict(inputs_information[0], inputs_information[1])
+
+        # Plot preprocessed data if did
+        if "sampling_grid" in self.config["inputs"] and self.config["inputs"]["sampling_grid"] is not None:
+            if self.config["inputs"]["sampling_grid"] == "reference_elev":
+                preprocessed_data = "plots/preprocessed_to_be_aligned_elev_map.png"
+            else:
+                preprocessed_data = "plots/preprocessed_reference_elev_map.png"
+
+            html += "<h2>Preprocessed DEM</h2>\n"
+            html += "<div style='display: flex'>\n"
+            html += (
+                "  <img src='" + preprocessed_data + "' alt='Image PNG' style='max-width: "
+                "50%; height: auto; width: 50%;'>\n"
+            )
             html += "</div>\n"
 
+        # Metadata: Inputs
+        for title, dictionary in list_dict[1:]:  # type: ignore
+            html += print_dict(title, dictionary)
+
         # Statistics table:
         if self.df_stats is not None:
             html += "<h2>Statistics</h2>\n"
-            html += self.df_stats.to_html(index=False)
+            html += "<table border='1' cellspacing='0' cellpadding='5'>\n"
+
+            # Plot one stat by row
+            df_cols = "".join([f'<td style="font-weight:bold">{col}</td>' for col in self.df_stats.T.columns])
+            html += f'<tr><td style="font-weight:bold">Data</td>{df_cols}</tr>\n'
+
+            for key, value in self.df_stats.T.iterrows():
+                df_values = "".join([f"<td>{str(val)}</td>" for val in value.values])
+                html += f'<tr><td style="font-weight:bold">{key}</td>{df_values}</tr>\n'
+            html += "</table>\n"
 
         # Coregistration: Add elevation difference plot and histograms before/after
         if self.compute_coreg:
diff --git a/xdem/workflows/topo.py b/xdem/workflows/topo.py
index 2722641d..9e992d50 100644
--- a/xdem/workflows/topo.py
+++ b/xdem/workflows/topo.py
@@ -168,6 +168,7 @@ def generate_terrain_attributes_png(self) -> None:
             ax.set_xticks([])
             ax.set_yticks([])
 
+        [fig.delaxes(ax) for ax in axes.flatten() if not ax.has_data()]
         plt.tight_layout()
         plt.savefig(self.outputs_folder / "plots" / "terrain_attributes_map.png", dpi=300)
         plt.close()