diff --git a/CHANGELOG.md b/CHANGELOG.md
index 5e16cd70..de73db04 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -12,6 +12,10 @@ and this project adheres to [Semantic Versioning][].
 
 ### Added
 
+-   Pushed `get_extent` functionality upstream to `spatialdata` (#162)
+-   Multiscale image handling: user can specify a scale, else the best scale is selected automatically given the figure size and dpi (#164)
+-   Large images are automatically rasterized to speed up performance (#164)
+
 ### Fixed
 
 ## [0.0.6] - 2023-11-06
diff --git a/src/spatialdata_plot/pl/basic.py b/src/spatialdata_plot/pl/basic.py
index c941cae1..8aa18061 100644
--- a/src/spatialdata_plot/pl/basic.py
+++ b/src/spatialdata_plot/pl/basic.py
@@ -41,6 +41,7 @@
 )
 from spatialdata_plot.pl.utils import (
     _get_cs_contents,
+    _get_valid_cs,
     _maybe_set_colors,
     _mpl_ax_contains_elements,
     _prepare_cmap_norm,
@@ -318,6 +319,7 @@ def render_images(
         palette: str | list[str] | None = None,
         alpha: float = 1.0,
         quantiles_for_norm: tuple[float | None, float | None] = (None, None),
+        scale: str | list[str] | None = None,
         **kwargs: Any,
     ) -> sd.SpatialData:
         """
@@ -340,6 +342,15 @@ def render_images(
             Alpha value for the shapes.
         quantiles_for_norm
             Tuple of (pmin, pmax) which will be used for quantile normalization.
+        scale
+            Influences the resolution of the rendering. Possibilities for setting this parameter:
+                1) None (default). The image is rasterized to fit the canvas size. For multiscale images, the best scale
+                is selected before the rasterization step.
+                2) Name of one of the scales in the multiscale image to be rendered. This scale is rendered as it is
+                (exception: a dpi is specified in `show()`. Then the image is rasterized to fit the canvas and dpi).
+                3) "full": render the full image without rasterization. In the case of a multiscale image, the scale
+                with the highest resolution is selected. This can lead to long computing times for large images!
+                4) List that is matched to the list of elements (can contain `None`, scale names or "full").
         kwargs
             Additional arguments to be passed to cmap and norm.
 
@@ -383,6 +394,7 @@ def render_images(
             palette=palette,
             alpha=alpha,
             quantiles_for_norm=quantiles_for_norm,
+            scale=scale,
         )
 
         return sdata
@@ -401,6 +413,7 @@ def render_labels(
         na_color: str | tuple[float, ...] | None = (0.0, 0.0, 0.0, 0.0),
         outline_alpha: float = 1.0,
         fill_alpha: float = 0.3,
+        scale: str | list[str] | None = None,
         **kwargs: Any,
     ) -> sd.SpatialData:
         """
@@ -433,6 +446,15 @@ def render_labels(
             Color to be used for NAs values, if present.
         alpha
             Alpha value for the labels.
+        scale
+            Influences the resolution of the rendering. Possibilities for setting this parameter:
+                1) None (default). The image is rasterized to fit the canvas size. For multiscale images, the best scale
+                is selected before the rasterization step.
+                2) Name of one of the scales in the multiscale image to be rendered. This scale is rendered as it is
+                (exception: a dpi is specified in `show()`. Then the image is rasterized to fit the canvas and dpi).
+                3) "full": render the full image without rasterization. In the case of a multiscale image, the scale
+                with the highest resolution is selected. This can lead to long computing times for large images!
+                4) List that is matched to the list of elements (can contain `None`, scale names or "full").
         kwargs
             Additional arguments to be passed to cmap and norm.
 
@@ -470,6 +492,7 @@ def render_labels(
             outline_alpha=outline_alpha,
             fill_alpha=fill_alpha,
             transfunc=kwargs.get("transfunc", None),
+            scale=scale,
         )
 
         return sdata
@@ -502,15 +525,22 @@ def show(
 
         Parameters
         ----------
+        coordinate_systems :
+            Name(s) of the coordinate system(s) to be plotted. If None, all coordinate systems are plotted.
+            If a coordinate system doesn't contain any relevant elements (as specified in the render_* calls),
+            it is automatically not plotted.
+        figsize :
+            Size of the figure (width, height) in inches. The size of the actual canvas may deviate from this,
+            depending on the dpi! In matplotlib, the actual figure size (in pixels) is dpi * figsize.
+            If None, the default of matlotlib is used (6.4, 4.8)
+        dpi :
+            Resolution of the plot in dots per inch (as in matplotlib).
+            If None, the default of matplotlib is used (100.0).
         ax :
             Matplotlib axes object to plot on. If None, a new figure is created.
             Works only if there is one image in the SpatialData object.
         ncols :
             Number of columns in the figure. Default is 4.
-        width :
-            Width of each subplot. Default is 4.
-        height :
-            Height of each subplot. Default is 3.
 
         Returns
         -------
@@ -576,6 +606,42 @@ def show(
             if cs not in sdata.coordinate_systems:
                 raise ValueError(f"Unknown coordinate system '{cs}', valid choices are: {sdata.coordinate_systems}")
 
+        # Check if user specified only certain elements to be plotted
+        cs_contents = _get_cs_contents(sdata)
+        elements_to_be_rendered = []
+        for cmd, params in render_cmds.items():
+            if cmd == "render_images" and cs_contents.query(f"cs == '{cs}'")["has_images"][0]:  # noqa: SIM114
+                if params.elements is not None:
+                    elements_to_be_rendered += (
+                        [params.elements] if isinstance(params.elements, str) else params.elements
+                    )
+            elif cmd == "render_shapes" and cs_contents.query(f"cs == '{cs}'")["has_shapes"][0]:  # noqa: SIM114
+                if params.elements is not None:
+                    elements_to_be_rendered += (
+                        [params.elements] if isinstance(params.elements, str) else params.elements
+                    )
+            elif cmd == "render_points" and cs_contents.query(f"cs == '{cs}'")["has_points"][0]:  # noqa: SIM114
+                if params.elements is not None:
+                    elements_to_be_rendered += (
+                        [params.elements] if isinstance(params.elements, str) else params.elements
+                    )
+            elif cmd == "render_labels" and cs_contents.query(f"cs == '{cs}'")["has_labels"][0]:  # noqa: SIM102
+                if params.elements is not None:
+                    elements_to_be_rendered += (
+                        [params.elements] if isinstance(params.elements, str) else params.elements
+                    )
+
+        # filter out cs without relevant elements
+        coordinate_systems = _get_valid_cs(
+            sdata=sdata,
+            coordinate_systems=coordinate_systems,
+            render_images="render_images" in render_cmds,
+            render_labels="render_labels" in render_cmds,
+            render_points="render_points" in render_cmds,
+            render_shapes="render_shapes" in render_cmds,
+            elements=elements_to_be_rendered,
+        )
+
         # set up canvas
         fig_params, scalebar_params = _prepare_params_plot(
             num_panels=len(coordinate_systems),
@@ -616,64 +682,70 @@ def show(
 
             for cmd, params in render_cmds.items():
                 if cmd == "render_images" and has_images:
-                    _render_images(
-                        sdata=sdata,
-                        render_params=params,
-                        coordinate_system=cs,
-                        ax=ax,
-                        fig_params=fig_params,
-                        scalebar_params=scalebar_params,
-                        legend_params=legend_params,
-                    )
                     wants_images = True
                     wanted_images = params.elements if params.elements is not None else list(sdata.images.keys())
-                    wanted_elements.extend(
-                        [
-                            image
-                            for image in wanted_images
-                            if cs in set(get_transformation(sdata.images[image], get_all=True).keys())
-                        ]
-                    )
+                    wanted_images_on_this_cs = [
+                        image
+                        for image in wanted_images
+                        if cs in set(get_transformation(sdata.images[image], get_all=True).keys())
+                    ]
+                    wanted_elements.extend(wanted_images_on_this_cs)
+                    if len(wanted_images_on_this_cs) > 0:
+                        rasterize = (params.scale is None) or (
+                            isinstance(params.scale, str)
+                            and params.scale != "full"
+                            and (dpi is not None or figsize is not None)
+                        )
+                        _render_images(
+                            sdata=sdata,
+                            render_params=params,
+                            coordinate_system=cs,
+                            ax=ax,
+                            fig_params=fig_params,
+                            scalebar_params=scalebar_params,
+                            legend_params=legend_params,
+                            rasterize=rasterize,
+                        )
 
                 elif cmd == "render_shapes" and has_shapes:
-                    _render_shapes(
-                        sdata=sdata,
-                        render_params=params,
-                        coordinate_system=cs,
-                        ax=ax,
-                        fig_params=fig_params,
-                        scalebar_params=scalebar_params,
-                        legend_params=legend_params,
-                    )
                     wants_shapes = True
                     wanted_shapes = params.elements if params.elements is not None else list(sdata.shapes.keys())
-                    wanted_elements.extend(
-                        [
-                            shape
-                            for shape in wanted_shapes
-                            if cs in set(get_transformation(sdata.shapes[shape], get_all=True).keys())
-                        ]
-                    )
+                    wanted_shapes_on_this_cs = [
+                        shape
+                        for shape in wanted_shapes
+                        if cs in set(get_transformation(sdata.shapes[shape], get_all=True).keys())
+                    ]
+                    wanted_elements.extend(wanted_shapes_on_this_cs)
+                    if len(wanted_shapes_on_this_cs) > 0:
+                        _render_shapes(
+                            sdata=sdata,
+                            render_params=params,
+                            coordinate_system=cs,
+                            ax=ax,
+                            fig_params=fig_params,
+                            scalebar_params=scalebar_params,
+                            legend_params=legend_params,
+                        )
 
                 elif cmd == "render_points" and has_points:
-                    _render_points(
-                        sdata=sdata,
-                        render_params=params,
-                        coordinate_system=cs,
-                        ax=ax,
-                        fig_params=fig_params,
-                        scalebar_params=scalebar_params,
-                        legend_params=legend_params,
-                    )
                     wants_points = True
                     wanted_points = params.elements if params.elements is not None else list(sdata.points.keys())
-                    wanted_elements.extend(
-                        [
-                            point
-                            for point in wanted_points
-                            if cs in set(get_transformation(sdata.points[point], get_all=True).keys())
-                        ]
-                    )
+                    wanted_points_on_this_cs = [
+                        point
+                        for point in wanted_points
+                        if cs in set(get_transformation(sdata.points[point], get_all=True).keys())
+                    ]
+                    wanted_elements.extend(wanted_points_on_this_cs)
+                    if len(wanted_points_on_this_cs) > 0:
+                        _render_points(
+                            sdata=sdata,
+                            render_params=params,
+                            coordinate_system=cs,
+                            ax=ax,
+                            fig_params=fig_params,
+                            scalebar_params=scalebar_params,
+                            legend_params=legend_params,
+                        )
 
                 elif cmd == "render_labels" and has_labels:
                     if sdata.table is not None and isinstance(params.color, str):
@@ -685,24 +757,30 @@ def show(
                                 key=params.color,
                                 palette=params.palette,
                             )
-                    _render_labels(
-                        sdata=sdata,
-                        render_params=params,
-                        coordinate_system=cs,
-                        ax=ax,
-                        fig_params=fig_params,
-                        scalebar_params=scalebar_params,
-                        legend_params=legend_params,
-                    )
                     wants_labels = True
                     wanted_labels = params.elements if params.elements is not None else list(sdata.labels.keys())
-                    wanted_elements.extend(
-                        [
-                            label
-                            for label in wanted_labels
-                            if cs in set(get_transformation(sdata.labels[label], get_all=True).keys())
-                        ]
-                    )
+                    wanted_labels_on_this_cs = [
+                        label
+                        for label in wanted_labels
+                        if cs in set(get_transformation(sdata.labels[label], get_all=True).keys())
+                    ]
+                    wanted_elements.extend(wanted_labels_on_this_cs)
+                    if len(wanted_labels_on_this_cs) > 0:
+                        rasterize = (params.scale is None) or (
+                            isinstance(params.scale, str)
+                            and params.scale != "full"
+                            and (dpi is not None or figsize is not None)
+                        )
+                        _render_labels(
+                            sdata=sdata,
+                            render_params=params,
+                            coordinate_system=cs,
+                            ax=ax,
+                            fig_params=fig_params,
+                            scalebar_params=scalebar_params,
+                            legend_params=legend_params,
+                            rasterize=rasterize,
+                        )
 
                 if title is None:
                     t = cs
diff --git a/src/spatialdata_plot/pl/render.py b/src/spatialdata_plot/pl/render.py
index c411bb0e..03d95662 100644
--- a/src/spatialdata_plot/pl/render.py
+++ b/src/spatialdata_plot/pl/render.py
@@ -11,15 +11,14 @@
 import numpy as np
 import pandas as pd
 import scanpy as sc
-import spatial_image
 import spatialdata as sd
 from anndata import AnnData
 from matplotlib.colors import ListedColormap, Normalize
+from multiscale_spatial_image.multiscale_spatial_image import MultiscaleSpatialImage
 from pandas.api.types import is_categorical_dtype
 from scanpy._settings import settings as sc_settings
+from spatialdata._core.data_extent import get_extent
 from spatialdata.models import (
-    Image2DModel,
-    Labels2DModel,
     PointsModel,
 )
 from spatialdata.transformations import (
@@ -43,7 +42,9 @@
     _get_linear_colormap,
     _map_color_seg,
     _maybe_set_colors,
+    _multiscale_to_spatial_image,
     _normalize,
+    _rasterize_if_necessary,
     _set_color_source_vec,
     to_hex,
 )
@@ -317,6 +318,7 @@ def _render_images(
     fig_params: FigParams,
     scalebar_params: ScalebarParams,
     legend_params: LegendParams,
+    rasterize: bool,
 ) -> None:
     elements = render_params.elements
 
@@ -331,11 +333,31 @@ def _render_images(
     if elements is None:
         elements = list(sdata_filt.images.keys())
 
-    for e in elements:
+    for i, e in enumerate(elements):
         img = sdata.images[e]
-        if not isinstance(img, spatial_image.SpatialImage):
-            img = Image2DModel.parse(img["scale0"].ds.to_array().squeeze(axis=0))
-            logger.warning(f"Multi-scale images not yet supported, using scale0 of multi-scale image '{e}'.")
+        extent = get_extent(img, coordinate_system=coordinate_system)
+        scale = render_params.scale[i] if isinstance(render_params.scale, list) else render_params.scale
+
+        # get best scale out of multiscale image
+        if isinstance(img, MultiscaleSpatialImage):
+            img = _multiscale_to_spatial_image(
+                multiscale_image=img,
+                element=e,
+                dpi=fig_params.fig.dpi,
+                width=fig_params.fig.get_size_inches()[0],
+                height=fig_params.fig.get_size_inches()[1],
+                scale=scale,
+            )
+        # rasterize spatial image if necessary to speed up performance
+        if rasterize:
+            img = _rasterize_if_necessary(
+                image=img,
+                dpi=fig_params.fig.dpi,
+                width=fig_params.fig.get_size_inches()[0],
+                height=fig_params.fig.get_size_inches()[1],
+                coordinate_system=coordinate_system,
+                extent=extent,
+            )
 
         if render_params.channel is None:
             channels = img.coords["c"].values
@@ -362,7 +384,7 @@ def _render_images(
             raise ValueError("If 'cmap' is provided, its length must match the number of channels.")
 
         # prepare transformations
-        trans = get_transformation(sdata.images[e], get_all=True)[coordinate_system]
+        trans = get_transformation(img, get_all=True)[coordinate_system]
         affine_trans = trans.to_affine_matrix(input_axes=("x", "y"), output_axes=("x", "y"))
         trans = mtransforms.Affine2D(matrix=affine_trans)
         trans_data = trans + ax.transData
@@ -488,6 +510,7 @@ def _render_labels(
     fig_params: FigParams,
     scalebar_params: ScalebarParams,
     legend_params: LegendParams,
+    rasterize: bool,
 ) -> None:
     elements = render_params.elements
 
@@ -507,11 +530,32 @@ def _render_labels(
     if elements is None:
         elements = list(sdata_filt.labels.keys())
 
-    for e in elements:
+    for i, e in enumerate(elements):
         label = sdata_filt.labels[e]
-        if not isinstance(label, spatial_image.SpatialImage):
-            label = Labels2DModel.parse(label["scale0"].ds.to_array().squeeze(axis=0))
-            logger.warning(f"Multi-scale labels not yet supported, using scale0 of multi-scale label '{e}'.")
+        extent = get_extent(label, coordinate_system=coordinate_system)
+        scale = render_params.scale[i] if isinstance(render_params.scale, list) else render_params.scale
+
+        # get best scale out of multiscale label
+        if isinstance(label, MultiscaleSpatialImage):
+            label = _multiscale_to_spatial_image(
+                multiscale_image=label,
+                element=e,
+                dpi=fig_params.fig.dpi,
+                width=fig_params.fig.get_size_inches()[0],
+                height=fig_params.fig.get_size_inches()[1],
+                scale=scale,
+                is_label=True,
+            )
+        # rasterize spatial image if necessary to speed up performance
+        if rasterize:
+            label = _rasterize_if_necessary(
+                image=label,
+                dpi=fig_params.fig.dpi,
+                width=fig_params.fig.get_size_inches()[0],
+                height=fig_params.fig.get_size_inches()[1],
+                coordinate_system=coordinate_system,
+                extent=extent,
+            )
 
         if sdata.table is None:
             instance_id = np.unique(label)
@@ -525,7 +569,7 @@ def _render_labels(
             # get instance id based on subsetted table
             instance_id = table.obs[instance_key].values
 
-        trans = get_transformation(sdata.labels[e], get_all=True)[coordinate_system]
+        trans = get_transformation(label, get_all=True)[coordinate_system]
         affine_trans = trans.to_affine_matrix(input_axes=("x", "y"), output_axes=("x", "y"))
         trans = mtransforms.Affine2D(matrix=affine_trans)
         trans_data = trans + ax.transData
@@ -533,7 +577,7 @@ def _render_labels(
         # get color vector (categorical or continuous)
         color_source_vector, color_vector, categorical = _set_color_source_vec(
             sdata=sdata_filt,
-            element=sdata_filt.labels[e],
+            element=label,
             element_name=e,
             value_to_plot=render_params.color,
             layer=render_params.layer,
diff --git a/src/spatialdata_plot/pl/render_params.py b/src/spatialdata_plot/pl/render_params.py
index cca7bd58..e82dfc22 100644
--- a/src/spatialdata_plot/pl/render_params.py
+++ b/src/spatialdata_plot/pl/render_params.py
@@ -106,6 +106,7 @@ class ImageRenderParams:
     palette: ListedColormap | str | None = None
     alpha: float = 1.0
     quantiles_for_norm: tuple[float | None, float | None] = (None, None)
+    scale: str | list[str] | None = None
 
 
 @dataclass
@@ -123,3 +124,4 @@ class LabelsRenderParams:
     outline_alpha: float = 1.0
     fill_alpha: float = 0.4
     transfunc: Callable[[float], float] | None = None
+    scale: str | list[str] | None = None
diff --git a/src/spatialdata_plot/pl/utils.py b/src/spatialdata_plot/pl/utils.py
index 5a66b667..60e84c1e 100644
--- a/src/spatialdata_plot/pl/utils.py
+++ b/src/spatialdata_plot/pl/utils.py
@@ -38,6 +38,7 @@
 from matplotlib.figure import Figure
 from matplotlib.gridspec import GridSpec
 from matplotlib_scalebar.scalebar import ScaleBar
+from multiscale_spatial_image.multiscale_spatial_image import MultiscaleSpatialImage
 from numpy.random import default_rng
 from pandas.api.types import CategoricalDtype, is_categorical_dtype
 from scanpy import settings
@@ -48,10 +49,12 @@
 from skimage.morphology import erosion, square
 from skimage.segmentation import find_boundaries
 from skimage.util import map_array
+from spatial_image import SpatialImage
+from spatialdata._core.operations.rasterize import rasterize
 from spatialdata._core.query.relational_query import _locate_value, get_values
 from spatialdata._logging import logger as logging
 from spatialdata._types import ArrayLike
-from spatialdata.models import Image2DModel, SpatialElement
+from spatialdata.models import Image2DModel, Labels2DModel, SpatialElement
 
 from spatialdata_plot.pl.render_params import (
     CmapParams,
@@ -114,6 +117,11 @@ def _prepare_params_plot(
         axs = None
         if ax is None:
             fig, ax = plt.subplots(figsize=figsize, dpi=dpi, constrained_layout=True)
+        elif isinstance(ax, Axes):
+            # needed for rasterization if user provides Axes object
+            fig = ax.get_figure()
+            fig.set_dpi(dpi)
+
     # set scalebar
     if scalebar_dx is not None:
         scalebar_dx, scalebar_units = _get_scalebar(scalebar_dx, scalebar_units, num_panels)
@@ -936,11 +944,17 @@ def _translate_image(
 ) -> spatial_image.SpatialImage:
     shifts: dict[str, int] = {axis: int(translation.translation[idx]) for idx, axis in enumerate(translation.axes)}
     img = image.values.copy()
+    # for yx images (important for rasterized MultiscaleImages as labels)
+    expanded_dims = False
+    if len(img.shape) == 2:
+        img = np.expand_dims(img, axis=0)
+        expanded_dims = True
+
     shifted_channels = []
 
     # split channels, shift axes individually, them recombine
-    if len(image.shape) == 3:
-        for c in range(image.shape[0]):
+    if len(img.shape) == 3:
+        for c in range(img.shape[0]):
             channel = img[c, :, :]
 
             # iterates over [x, y]
@@ -960,6 +974,12 @@ def _translate_image(
 
             shifted_channels.append(channel)
 
+    if expanded_dims:
+        return Labels2DModel.parse(
+            np.array(shifted_channels[0]),
+            dims=["y", "x"],
+            transformations=image.attrs["transform"],
+        )
     return Image2DModel.parse(
         np.array(shifted_channels),
         dims=["c", "y", "x"],
@@ -1031,3 +1051,180 @@ def _mpl_ax_contains_elements(ax: Axes) -> bool:
     return (
         len(ax.lines) > 0 or len(ax.collections) > 0 or len(ax.images) > 0 or len(ax.patches) > 0 or len(ax.tables) > 0
     )
+
+
+def _get_valid_cs(
+    sdata: sd.SpatialData,
+    coordinate_systems: Sequence[str],
+    render_images: bool,
+    render_labels: bool,
+    render_points: bool,
+    render_shapes: bool,
+    elements: list[str],
+) -> Sequence[str]:
+    """Get names of the valid coordinate systems.
+
+    Valid cs are cs that contain elements to be rendered:
+    1. In case the user specified elements:
+        all cs that contain at least one of those elements
+    2. Else:
+        all cs that contain at least one element that should
+        be rendered (depending on whether images/points/labels/...
+        should be rendered)
+    """
+    cs_mapping = _get_coordinate_system_mapping(sdata)
+    valid_cs = []
+    for cs in coordinate_systems:
+        if (len(elements) > 0 and any(e in elements for e in cs_mapping[cs])) or (
+            len(elements) == 0
+            and (
+                (len(sdata.images.keys()) > 0 and render_images)
+                or (len(sdata.labels.keys()) > 0 and render_labels)
+                or (len(sdata.points.keys()) > 0 and render_points)
+                or (len(sdata.shapes.keys()) > 0 and render_shapes)
+            )
+        ):  # not nice, but ruff wants it (SIM114)
+            valid_cs.append(cs)
+        else:
+            logging.info(f"Dropping coordinate system '{cs}' since it doesn't have relevant elements.")
+    return valid_cs
+
+
+def _rasterize_if_necessary(
+    image: SpatialImage,
+    dpi: float,
+    width: float,
+    height: float,
+    coordinate_system: str,
+    extent: dict[str, tuple[float, float]],
+) -> SpatialImage:
+    """Ensure fast rendering by adapting the resolution if necessary.
+
+    A SpatialImage is prepared for plotting. To improve performance, large images are rasterized.
+
+    Parameters
+    ----------
+    image
+        Input spatial image that should be rendered
+    dpi
+        Resolution of the figure
+    width
+        Width (in inches) of the figure
+    height
+        Height (in inches) of the figure
+    coordinate_system
+        name of the coordinate system the image belongs to
+    extent
+        extent of the (full size) image. Must be a dict containing a tuple with min and
+        max extent for the keys "x" and "y".
+
+    Returns
+    -------
+    SpatialImage
+        Spatial image ready for rendering
+    """
+    has_c_dim = len(image.shape) == 3
+    if has_c_dim:
+        y_dims = image.shape[1]
+        x_dims = image.shape[2]
+    else:
+        y_dims = image.shape[0]
+        x_dims = image.shape[1]
+
+    target_y_dims = dpi * height
+    target_x_dims = dpi * width
+
+    # TODO: when exactly do we want to rasterize?
+    do_rasterization = y_dims > target_y_dims + 100 or x_dims > target_x_dims + 100
+    if x_dims < 2000 and y_dims < 2000:
+        do_rasterization = False
+
+    if do_rasterization:
+        # TODO: do we want min here?
+        target_unit_to_pixels = min(target_y_dims / y_dims, target_x_dims / x_dims)
+        image = rasterize(
+            image,
+            ("y", "x"),
+            [extent["y"][0], extent["x"][0]],
+            [extent["y"][1], extent["x"][1]],
+            coordinate_system,
+            target_unit_to_pixels=target_unit_to_pixels,
+        )
+
+    return image
+
+
+def _multiscale_to_spatial_image(
+    multiscale_image: MultiscaleSpatialImage,
+    element: str,
+    dpi: float,
+    width: float,
+    height: float,
+    scale: str | None = None,
+    is_label: bool = False,
+) -> SpatialImage:
+    """Extract the SpatialImage to be rendered from a multiscale image.
+
+    From the `MultiscaleSpatialImage`, the scale that fits the given image size and dpi most is selected
+    and returned. In case the lowest resolution is still too high, a rasterization step is added.
+
+    Parameters
+    ----------
+    multiscale_image
+        `MultiscaleSpatialImage` that should be rendered
+    element
+        name of the multiscale image
+    dpi
+        dpi of the target image
+    width
+        width of the target image in inches
+    height
+        height of the target image in inches
+    scale
+        specific scale that the user chose, if None the heuristic is used
+    is_label
+        When True, the multiscale image contains labels which don't contain the `c` dimension
+
+    Returns
+    -------
+    SpatialImage
+        To be rendered, extracted from the MultiscaleSpatialImage respecting the dpi and size of the target image.
+    """
+    scales = [leaf.name for leaf in multiscale_image.leaves]
+    x_dims = [multiscale_image[scale].dims["x"] for scale in scales]
+    y_dims = [multiscale_image[scale].dims["y"] for scale in scales]
+
+    if isinstance(scale, str):
+        if scale not in scales and scale != "full":
+            raise ValueError(f'Scale {scale} does not exist. Please select one of {scales} or set scale = "full"!')
+        optimal_scale = scale
+        if scale == "full":
+            # use scale with highest resolution
+            optimal_scale = scales[np.argmax(x_dims)]
+    else:
+        # ensure that lists are sorted
+        order = np.argsort(x_dims)
+        scales = [scales[i] for i in order]
+        x_dims = [x_dims[i] for i in order]
+        y_dims = [y_dims[i] for i in order]
+
+        optimal_x = width * dpi
+        optimal_y = height * dpi
+
+        # get scale where the dimensions are close to the optimal values
+        # when possible, pick higher resolution (worst case: downscaled afterwards)
+        optimal_index_y = np.searchsorted(y_dims, optimal_y)
+        if optimal_index_y == len(y_dims):
+            optimal_index_y -= 1
+        optimal_index_x = np.searchsorted(x_dims, optimal_x)
+        if optimal_index_x == len(x_dims):
+            optimal_index_x -= 1
+
+        # pick the scale with higher resolution (worst case: downscaled afterwards)
+        optimal_scale = scales[min(optimal_index_x, optimal_index_y)]
+
+    # NOTE: problematic if there are cases with > 1 data variable
+    data_var_keys = list(multiscale_image[optimal_scale].data_vars)
+    image = multiscale_image[optimal_scale][data_var_keys[0]]
+
+    return Labels2DModel.parse(image) if is_label else Image2DModel.parse(image)
diff --git a/tests/_images/Images_can_do_rasterization.png b/tests/_images/Images_can_do_rasterization.png
new file mode 100644
index 00000000..7212de26
Binary files /dev/null and b/tests/_images/Images_can_do_rasterization.png differ
diff --git a/tests/_images/Images_can_render_given_scale_of_multiscale_image.png b/tests/_images/Images_can_render_given_scale_of_multiscale_image.png
new file mode 100644
index 00000000..e1ba8949
Binary files /dev/null and b/tests/_images/Images_can_render_given_scale_of_multiscale_image.png differ
diff --git a/tests/_images/Images_can_render_multiscale_image.png b/tests/_images/Images_can_render_multiscale_image.png
new file mode 100644
index 00000000..48553078
Binary files /dev/null and b/tests/_images/Images_can_render_multiscale_image.png differ
diff --git a/tests/_images/Images_can_stop_rasterization_with_scale_full.png b/tests/_images/Images_can_stop_rasterization_with_scale_full.png
new file mode 100644
index 00000000..35c7d84c
Binary files /dev/null and b/tests/_images/Images_can_stop_rasterization_with_scale_full.png differ
diff --git a/tests/_images/Labels_can_do_rasterization.png b/tests/_images/Labels_can_do_rasterization.png
new file mode 100644
index 00000000..a4028164
Binary files /dev/null and b/tests/_images/Labels_can_do_rasterization.png differ
diff --git a/tests/_images/Labels_can_render_given_scale_of_multiscale_labels.png b/tests/_images/Labels_can_render_given_scale_of_multiscale_labels.png
new file mode 100644
index 00000000..4ef7366a
Binary files /dev/null and b/tests/_images/Labels_can_render_given_scale_of_multiscale_labels.png differ
diff --git a/tests/_images/Labels_can_render_multiscale_labels.png b/tests/_images/Labels_can_render_multiscale_labels.png
new file mode 100644
index 00000000..125bcbff
Binary files /dev/null and b/tests/_images/Labels_can_render_multiscale_labels.png differ
diff --git a/tests/_images/Labels_can_stop_rasterization_with_scale_full.png b/tests/_images/Labels_can_stop_rasterization_with_scale_full.png
new file mode 100644
index 00000000..95f40fff
Binary files /dev/null and b/tests/_images/Labels_can_stop_rasterization_with_scale_full.png differ
diff --git a/tests/pl/test_render_images.py b/tests/pl/test_render_images.py
index 6189e189..a1b732ff 100644
--- a/tests/pl/test_render_images.py
+++ b/tests/pl/test_render_images.py
@@ -1,6 +1,8 @@
+import dask.array as da
 import matplotlib
 import scanpy as sc
 import spatialdata_plot  # noqa: F401
+from spatial_image import to_spatial_image
 from spatialdata import SpatialData
 
 from tests.conftest import PlotTester, PlotTesterMeta
@@ -49,3 +51,29 @@ def test_plot_can_pass_cmap_to_each_channel(self, sdata_blobs: SpatialData):
 
     def test_plot_can_normalize_image(self, sdata_blobs: SpatialData):
         sdata_blobs.pl.render_images(elements="blobs_image", quantiles_for_norm=(5, 90)).pl.show()
+
+    def test_plot_can_render_multiscale_image(self, sdata_blobs: SpatialData):
+        sdata_blobs.pl.render_images("blobs_multiscale_image").pl.show()
+
+    def test_plot_can_render_given_scale_of_multiscale_image(self, sdata_blobs: SpatialData):
+        sdata_blobs.pl.render_images("blobs_multiscale_image", scale="scale2").pl.show()
+
+    def test_plot_can_do_rasterization(self, sdata_blobs: SpatialData):
+        temp = sdata_blobs["blobs_image"].data.copy()
+        temp = da.concatenate([temp] * 6, axis=1)
+        temp = da.concatenate([temp] * 6, axis=2)
+        img = to_spatial_image(temp, dims=("c", "y", "x"))
+        img.attrs["transform"] = sdata_blobs["blobs_image"].transform
+        sdata_blobs["blobs_giant_image"] = img
+
+        sdata_blobs.pl.render_images("blobs_giant_image").pl.show()
+
+    def test_plot_can_stop_rasterization_with_scale_full(self, sdata_blobs: SpatialData):
+        temp = sdata_blobs["blobs_image"].data.copy()
+        temp = da.concatenate([temp] * 6, axis=1)
+        temp = da.concatenate([temp] * 6, axis=2)
+        img = to_spatial_image(temp, dims=("c", "y", "x"))
+        img.attrs["transform"] = sdata_blobs["blobs_image"].transform
+        sdata_blobs["blobs_giant_image"] = img
+
+        sdata_blobs.pl.render_images("blobs_giant_image", scale="full").pl.show()
diff --git a/tests/pl/test_render_labels.py b/tests/pl/test_render_labels.py
index aec09e4c..52dada28 100644
--- a/tests/pl/test_render_labels.py
+++ b/tests/pl/test_render_labels.py
@@ -1,6 +1,8 @@
+import dask.array as da
 import matplotlib
 import scanpy as sc
 import spatialdata_plot  # noqa: F401
+from spatial_image import to_spatial_image
 from spatialdata import SpatialData
 
 from tests.conftest import PlotTester, PlotTesterMeta
@@ -21,3 +23,39 @@
 class TestLabels(PlotTester, metaclass=PlotTesterMeta):
     def test_plot_can_render_labels(self, sdata_blobs: SpatialData):
         sdata_blobs.pl.render_labels(elements="blobs_labels").pl.show()
+
+    def test_plot_can_render_multiscale_labels(self, sdata_blobs: SpatialData):
+        sdata_blobs.table.obs["region"] = "blobs_multiscale_labels"
+        sdata_blobs.table.uns["spatialdata_attrs"]["region"] = "blobs_multiscale_labels"
+        sdata_blobs.pl.render_labels("blobs_multiscale_labels").pl.show()
+
+    def test_plot_can_render_given_scale_of_multiscale_labels(self, sdata_blobs: SpatialData):
+        sdata_blobs.table.obs["region"] = "blobs_multiscale_labels"
+        sdata_blobs.table.uns["spatialdata_attrs"]["region"] = "blobs_multiscale_labels"
+        sdata_blobs.pl.render_labels("blobs_multiscale_labels", scale="scale1").pl.show()
+
+    def test_plot_can_do_rasterization(self, sdata_blobs: SpatialData):
+        temp = sdata_blobs["blobs_labels"].data.copy()
+        temp = da.concatenate([temp] * 6, axis=0)
+        temp = da.concatenate([temp] * 6, axis=1)
+        img = to_spatial_image(temp, dims=("y", "x"))
+        img.attrs["transform"] = sdata_blobs["blobs_labels"].transform
+        sdata_blobs["blobs_giant_labels"] = img
+
+        sdata_blobs.table.obs["region"] = "blobs_giant_labels"
+        sdata_blobs.table.uns["spatialdata_attrs"]["region"] = "blobs_giant_labels"
+
+        sdata_blobs.pl.render_labels("blobs_giant_labels").pl.show()
+
+    def test_plot_can_stop_rasterization_with_scale_full(self, sdata_blobs: SpatialData):
+        temp = sdata_blobs["blobs_labels"].data.copy()
+        temp = da.concatenate([temp] * 6, axis=0)
+        temp = da.concatenate([temp] * 6, axis=1)
+        img = to_spatial_image(temp, dims=("y", "x"))
+        img.attrs["transform"] = sdata_blobs["blobs_labels"].transform
+        sdata_blobs["blobs_giant_labels"] = img
+
+        sdata_blobs.table.obs["region"] = "blobs_giant_labels"
+        sdata_blobs.table.uns["spatialdata_attrs"]["region"] = "blobs_giant_labels"
+
+        sdata_blobs.pl.render_labels("blobs_giant_labels", scale="full").pl.show()