diff --git a/docs/src/advanced/feature.md b/docs/src/advanced/feature.md
index 98877862..54f8e199 100644
--- a/docs/src/advanced/feature.md
+++ b/docs/src/advanced/feature.md
@@ -11,29 +11,42 @@ with Reader("my-tif.tif") as cog:
     img: ImageData = cog.feature(geojson_feature, max_size=1024)  # we limit the max_size to 1024
 ```
 
-Under the hood, the `.feature` method uses `GDALWarpVRT`'s `cutline` option and
-the `.part()` method. The below process is roughly what `.feature` does for you.
+Under the hood, the `.feature` method uses rasterio's [`rasterize`](https://rasterio.readthedocs.io/en/latest/api/rasterio.features.html#rasterio.features.rasterize)
+function and the `.part()` method. The below process is roughly what `.feature` does for you.
 
 ```python
+from rasterio.features import rasterize, bounds as featureBounds
+
 from rio_tiler.io import Reader
-from rio_tiler.utils import create_cutline
-from rasterio.features import bounds as featureBounds
 
 # Use Reader to open and read the dataset
 with Reader("my_tif.tif") as cog:
-    # Create WTT Cutline
-    cutline = create_cutline(cog.dataset, feat, geometry_crs="epsg:4326")
 
     # Get BBOX of the polygon
     bbox = featureBounds(feat)
 
-    # Read part of the data (bbox) and use the cutline to mask the data
-    data, mask = cog.part(bbox, vrt_options={'cutline': cutline}, max_size=1024)
+    # Read part of the data overlapping with the geometry bbox
+    # assuming that the geometry coordinates are in web mercator
+    img = cog.part(bbox, bounds_crs=f"EPSG:3857", max_size=1024)
+
+    # Rasterize geometry using the same geotransform parameters
+    cutline = rasterize(
+        [feat],
+        out_shape=(img.height, img.width),
+        transform=img.transform,
+        ...
+    )
+
+    # Apply geometry mask to imagery
+    img.array.mask = numpy.where(~cutline, img.array.mask, True)
 ```
 
-Another interesting fact about the `cutline` option is that it can be used with other methods:
+Another interesting way to cut features is to use the GDALWarpVRT's `cutline`
+option with the .part(), .preview(), or .tile() methods:
 
 ```python
+from rio_tiler.utils import create_cutline
+
 bbox = featureBounds(feat)
 
 # Use Reader to open and read the dataset
@@ -41,9 +54,13 @@ with Reader("my_tif.tif") as cog:
     # Create WTT Cutline
     cutline = create_cutline(cog.dataset, feat, geometry_crs="epsg:4326")
 
+    # Get a part of the geotiff but use the cutline to mask the data
+    bbox = featureBounds(feat)
+    img = cog.part(bbox, vrt_options={'cutline': cutline})
+
     # Get a preview of the whole geotiff but use the cutline to mask the data
-    data, mask = cog.preview(vrt_options={'cutline': cutline})
+    img = cog.preview(vrt_options={'cutline': cutline})
 
     # Read a mercator tile and use the cutline to mask the data
-    data, mask = cog.tile(1, 1, 1, vrt_options={'cutline': cutline})
+    img = cog.tile(1, 1, 1, vrt_options={'cutline': cutline})
 ```
diff --git a/rio_tiler/io/rasterio.py b/rio_tiler/io/rasterio.py
index 561670f5..ba20e7e3 100644
--- a/rio_tiler/io/rasterio.py
+++ b/rio_tiler/io/rasterio.py
@@ -13,7 +13,7 @@
 from rasterio import transform
 from rasterio.crs import CRS
 from rasterio.features import bounds as featureBounds
-from rasterio.features import geometry_mask
+from rasterio.features import geometry_mask, rasterize
 from rasterio.io import DatasetReader, DatasetWriter, MemoryFile
 from rasterio.rio.overview import get_maximum_overview_level
 from rasterio.transform import from_bounds as transform_from_bounds
@@ -28,13 +28,14 @@
     ExpressionMixingWarning,
     NoOverviewWarning,
     PointOutsideBounds,
+    RioTilerError,
     TileOutsideBounds,
 )
 from rio_tiler.expression import parse_expression
 from rio_tiler.io.base import BaseReader
 from rio_tiler.models import BandStatistics, ImageData, Info, PointData
 from rio_tiler.types import BBox, Indexes, NumType, RIOResampling
-from rio_tiler.utils import create_cutline, has_alpha_band, has_mask_band
+from rio_tiler.utils import _validate_shape_input, has_alpha_band, has_mask_band
 
 
 @attr.s
@@ -501,6 +502,7 @@ def feature(
         height: Optional[int] = None,
         width: Optional[int] = None,
         buffer: Optional[NumType] = None,
+        all_touched: Optional[bool] = False,
         **kwargs: Any,
     ) -> ImageData:
         """Read part of a Dataset defined by a geojson feature.
@@ -515,23 +517,24 @@ def feature(
             height (int, optional): Output height of the array.
             width (int, optional): Output width of the array.
             buffer (int or float, optional): Buffer on each side of the given aoi. It must be a multiple of `0.5`. Output **image size** will be expanded to `output imagesize + 2 * buffer` (e.g 0.5 = 257x257, 1.0 = 258x258).
+            all_touched (bool, optional): Shape rasterization parameter for mask. If True, all pixels touched the input shape will not be masked. If false, only pixels whose center is within the input shape or that are selected by Bresenham's line algorithm will not be masked.
             kwargs (optional): Options to forward to the `Reader.part` method.
 
         Returns:
             rio_tiler.models.ImageData: ImageData instance with data, mask and input spatial info.
 
         """
+        shape = _validate_shape_input(shape)
+
         if not dst_crs:
             dst_crs = shape_crs
 
         # Get BBOX of the polygon
         bbox = featureBounds(shape)
 
-        cutline = create_cutline(self.dataset, shape, geometry_crs=shape_crs)
         vrt_options = kwargs.pop("vrt_options", {})
-        vrt_options.update({"cutline": cutline})
 
-        return self.part(
+        img = self.part(
             bbox,
             dst_crs=dst_crs,
             bounds_crs=shape_crs,
@@ -545,6 +548,20 @@ def feature(
             **kwargs,
         )
 
+        cutline_mask = rasterize(
+            [shape],
+            out_shape=(img.height, img.width),
+            transform=img.transform,
+            all_touched=all_touched,
+            default_value=0,
+            fill=1,
+            dtype="uint8",
+        ).astype("bool")
+        img.cutline_mask = cutline_mask
+        img.array.mask = numpy.where(~cutline_mask, img.array.mask, True)
+
+        return img
+
     def read(
         self,
         indexes: Optional[Indexes] = None,
diff --git a/rio_tiler/io/xarray.py b/rio_tiler/io/xarray.py
index a62aeef6..b6e3879f 100644
--- a/rio_tiler/io/xarray.py
+++ b/rio_tiler/io/xarray.py
@@ -8,17 +8,17 @@
 from morecantile import Tile, TileMatrixSet
 from rasterio.crs import CRS
 from rasterio.enums import Resampling
-from rasterio.features import is_valid_geom
 from rasterio.rio.overview import get_maximum_overview_level
 from rasterio.transform import from_bounds, rowcol
 from rasterio.warp import calculate_default_transform
 from rasterio.warp import transform as transform_coords
 
 from rio_tiler.constants import WEB_MERCATOR_TMS, WGS84_CRS
-from rio_tiler.errors import PointOutsideBounds, RioTilerError, TileOutsideBounds
+from rio_tiler.errors import PointOutsideBounds, TileOutsideBounds
 from rio_tiler.io.base import BaseReader
 from rio_tiler.models import BandStatistics, ImageData, Info, PointData
 from rio_tiler.types import BBox, WarpResampling
+from rio_tiler.utils import _validate_shape_input
 
 try:
     import xarray
@@ -373,11 +373,7 @@ def feature(
         if not dst_crs:
             dst_crs = shape_crs
 
-        if "geometry" in shape:
-            shape = shape["geometry"]
-
-        if not is_valid_geom(shape):
-            raise RioTilerError("Invalid geometry")
+        shape = _validate_shape_input(shape)
 
         ds = self.input.rio.clip([shape], crs=shape_crs)
 
diff --git a/rio_tiler/models.py b/rio_tiler/models.py
index 1ec276ae..ea937ce4 100644
--- a/rio_tiler/models.py
+++ b/rio_tiler/models.py
@@ -334,6 +334,7 @@ class ImageData:
     dataset_statistics: Optional[Sequence[Tuple[float, float]]] = attr.ib(
         default=None, kw_only=True
     )
+    cutline_mask: Optional[numpy.ndarray] = attr.ib(default=None)
 
     @band_names.default
     def _default_names(self):
@@ -429,12 +430,15 @@ def create_from_list(cls, data: Sequence["ImageData"]) -> "ImageData":
             max_h, max_w = max(h), max(w)
             for img in data:
                 if img.height == max_h and img.width == max_w:
+                    cutline_mask = img.cutline_mask
                     continue
                 arr = numpy.ma.MaskedArray(
                     resize_array(img.array.data, max_h, max_w),
                     mask=resize_array(img.array.mask * 1, max_h, max_w).astype("bool"),
                 )
                 img.array = arr
+        else:
+            cutline_mask = data[0].cutline_mask
 
         arr = numpy.ma.concatenate([img.array for img in data])
 
@@ -472,6 +476,7 @@ def create_from_list(cls, data: Sequence["ImageData"]) -> "ImageData":
             bounds=bounds,
             band_names=band_names,
             dataset_statistics=dataset_statistics,
+            cutline_mask=cutline_mask,
         )
 
     def as_masked(self) -> numpy.ma.MaskedArray:
diff --git a/rio_tiler/mosaic/methods/base.py b/rio_tiler/mosaic/methods/base.py
index fcb4601c..269ad15f 100644
--- a/rio_tiler/mosaic/methods/base.py
+++ b/rio_tiler/mosaic/methods/base.py
@@ -13,6 +13,7 @@ class MosaicMethodBase(abc.ABC):
 
     mosaic: Optional[numpy.ma.MaskedArray] = field(default=None, init=False)
     exit_when_filled: bool = field(default=False, init=False)
+    cutline_mask: Optional[numpy.ndarray] = field(default=None, init=False)
 
     @property
     def is_done(self) -> bool:
@@ -25,8 +26,15 @@ def is_done(self) -> bool:
         if self.mosaic is None:
             return False
 
-        if self.exit_when_filled and not numpy.ma.is_masked(self.mosaic):
-            return True
+        if self.exit_when_filled:
+            if (
+                self.cutline_mask is not None
+                and numpy.sum(numpy.where(~self.cutline_mask, self.mosaic.mask, False))
+                == 0
+            ):
+                return True
+            elif not numpy.ma.is_masked(self.mosaic):
+                return True
 
         return False
 
diff --git a/rio_tiler/mosaic/reader.py b/rio_tiler/mosaic/reader.py
index b96e798b..a37a8dfd 100644
--- a/rio_tiler/mosaic/reader.py
+++ b/rio_tiler/mosaic/reader.py
@@ -93,6 +93,8 @@ def mosaic_reader(
             bounds = img.bounds
             band_names = img.band_names
 
+            pixel_selection.cutline_mask = img.cutline_mask
+
             assets_used.append(asset)
             pixel_selection.feed(img.array)
 
diff --git a/rio_tiler/utils.py b/rio_tiler/utils.py
index fa93b25c..f43aec2f 100644
--- a/rio_tiler/utils.py
+++ b/rio_tiler/utils.py
@@ -537,12 +537,7 @@ def create_cutline(
         str: WKT geometry in form of `POLYGON ((x y, x y, ...)))
 
     """
-    if "geometry" in geometry:
-        geometry = geometry["geometry"]
-
-    if not is_valid_geom(geometry):
-        raise RioTilerError("Invalid geometry")
-
+    geometry = _validate_shape_input(geometry)
     geom_type = geometry["type"]
 
     if geometry_crs:
@@ -636,3 +631,15 @@ def normalize_bounds(bounds: BBox) -> BBox:
         max(bounds[0], bounds[2]),
         max(bounds[1], bounds[3]),
     )
+
+
+def _validate_shape_input(shape: Dict) -> Dict:
+    """Ensure input shape is valid and reduce features to geometry"""
+
+    if "geometry" in shape:
+        shape = shape["geometry"]
+
+    if not is_valid_geom(shape):
+        raise RioTilerError("Invalid geometry")
+
+    return shape
diff --git a/tests/test_mosaic.py b/tests/test_mosaic.py
index f6c75692..b16948e9 100644
--- a/tests/test_mosaic.py
+++ b/tests/test_mosaic.py
@@ -6,6 +6,7 @@
 import numpy
 import pytest
 import rasterio
+from rasterio.crs import CRS
 from rasterio.warp import transform_bounds
 
 from rio_tiler import mosaic
@@ -62,6 +63,12 @@ def _read_point(src_path: str, *args, **kwargs) -> PointData:
         return src.point(*args, **kwargs)
 
 
+def _read_feature(src_path: str, *args, **kwargs) -> ImageData:
+    """Read feature from an asset"""
+    with Reader(src_path) as src:
+        return src.feature(*args, **kwargs)
+
+
 def test_mosaic_tiler():
     """Test mosaic tiler."""
     # test with default and full covered tile and default options
@@ -577,3 +584,123 @@ def test_mosaic_methods_last(m):
     assert img.mask.all()
     assert img.crs == WEB_MERCATOR_TMS.crs
     assert img.bounds == WEB_MERCATOR_TMS.xy_bounds(x, y, z)
+
+
+def test_mosaic_feature():
+    """Test mosaic feature."""
+    cog1 = {
+        "type": "Feature",
+        "geometry": {
+            "type": "Polygon",
+            "coordinates": [
+                [
+                    [479966, 5076955],
+                    [485507, 5083109],
+                    [490381, 5078849],
+                    [488014, 5073641],
+                    [479966, 5076955],
+                ]
+            ],
+        },
+    }
+    cog2 = {
+        "type": "Feature",
+        "geometry": {
+            "type": "Polygon",
+            "coordinates": [
+                [
+                    [700570, 5042202],
+                    [706585, 5046212],
+                    [713045, 5040643],
+                    [709035, 5036410],
+                    [700570, 5042202],
+                ]
+            ],
+        },
+    }
+    both = {
+        "type": "Polygon",
+        "coordinates": [
+            [
+                [575815, 5049776],
+                [584726, 5055123],
+                [589405, 5046212],
+                [581830, 5046212],
+                [575815, 5049776],
+            ]
+        ],
+    }
+    edge = {
+        "type": "Polygon",
+        "coordinates": [
+            [
+                [614608, 5021770],
+                [633407, 5021770],
+                [633407, 5030735],
+                [614608, 5030735],
+                [614608, 5021770],
+            ]
+        ],
+    }
+    away = {
+        "type": "Polygon",
+        "coordinates": [
+            [
+                [900570, 5042202],
+                [906585, 5046212],
+                [913045, 5040643],
+                [909035, 5036410],
+                [900570, 5042202],
+            ]
+        ],
+    }
+    crs = CRS.from_epsg(32618)
+
+    dat, _ = mosaic.mosaic_reader(assets, _read_feature, shape=cog1, shape_crs=crs)
+    assert dat.data.shape == (3, 32, 35)
+    assert list(numpy.unique(dat.array.mask)) == [False, True]
+    assert dat.assets == [asset1]
+    assert dat.crs == crs
+
+    dat, _ = mosaic.mosaic_reader(
+        assets_order, _read_feature, shape=cog2, shape_crs=crs
+    )
+    assert dat.data.shape == (3, 33, 42)
+    assert list(numpy.unique(dat.array.mask)) == [False, True]
+    assert dat.assets == [asset2]
+    assert dat.crs == crs
+
+    dat, _ = mosaic.mosaic_reader(assets, _read_feature, shape=both, shape_crs=crs)
+    assert dat.data.shape == (3, 30, 45)
+    assert list(numpy.unique(dat.array.mask)) == [False, True]
+    assert dat.assets == [asset1]  # Covers both but finishes early
+    assert dat.crs == crs
+
+    dat, _ = mosaic.mosaic_reader(assets, _read_feature, shape=edge, shape_crs=crs)
+    assert dat.data.shape == (3, 30, 63)
+    assert list(numpy.unique(dat.array.mask)) == [
+        False
+    ]  # Squared polygon, fills all pixels
+    assert dat.assets == [
+        asset1,
+        asset2,
+    ]  # At edge of asset 1, will need asset 2 to complete
+    assert dat.crs == crs
+
+    dat, _ = mosaic.mosaic_reader(assets, _read_feature, shape=away, shape_crs=crs)
+    assert dat.data.shape == (3, 33, 42)
+    assert list(numpy.unique(dat.array.mask)) == [True]
+    assert dat.assets == [asset1, asset2]
+    assert dat.crs == crs
+
+    dat, _ = mosaic.mosaic_reader(
+        assets,
+        _read_feature,
+        shape=cog1,
+        shape_crs=crs,
+        pixel_selection=defaults.HighestMethod,
+    )
+    assert dat.data.shape == (3, 32, 35)
+    assert list(numpy.unique(dat.array.mask)) == [False, True]
+    assert dat.assets == [asset1, asset2]
+    assert dat.crs == crs