diff --git a/src/astroimred/fitsmgmt/table.py b/src/astroimred/fitsmgmt/table.py
index 066c663..fba77ea 100644
--- a/src/astroimred/fitsmgmt/table.py
+++ b/src/astroimred/fitsmgmt/table.py
@@ -749,7 +749,7 @@ def _check_mismatch(row, keys, values):
         summary_table = inputs.to_pandas()
         fitslist = summary_table[table_filecol].to_list()
     elif isinstance(inputs, pd.DataFrame):
-        summary_table = inputs
+        summary_table = inputs.copy()
         fitslist = summary_table[table_filecol].to_list()
     else:
         # No need to sort here because the real "sort" will be done later in fits_summary
@@ -788,12 +788,13 @@ def _check_mismatch(row, keys, values):
     if selecting:
         for k, v in zip(type_key, type_val, strict=False):
             if isinstance(v, str):
-                match_mask = summary_table[k].str.match(v)
-                summary_table = summary_table[match_mask]
-                fitslist = np.array(fitslist)[match_mask].tolist()
-                # NOTE: Is there a better way to do this?
-                with contextlib.suppress(ValueError):
-                    summary_table.reset_index(inplace=True, drop=True)
+                match_mask = summary_table[k].astype("string").str.match(v, na=False)
+                summary_table = summary_table[match_mask].reset_index(drop=True)
+                fitslist = [
+                    item
+                    for item, keep in zip(fitslist, match_mask, strict=False)
+                    if keep
+                ]
             else:  # not used as regex
                 _type_key.append(k)
                 _type_val.append(v)
diff --git a/src/astroimred/imutil/_util_fits.py b/src/astroimred/imutil/_util_fits.py
index 20a07a4..baf077a 100644
--- a/src/astroimred/imutil/_util_fits.py
+++ b/src/astroimred/imutil/_util_fits.py
@@ -456,7 +456,7 @@ def extract_stack_metadata(
     hdr0 = _parse_imc_data_header(items[0], extension=extension, parse_data=False)[1]
     if hdr0 is None:
         raise ValueError("Could not read header from the first input image.")
-    ndim = hdr0["NAXIS"]
+    ndim = items[0].data.ndim if isinstance(items[0], CCDData) else hdr0["NAXIS"]
     # N x ndim. sizes[i, :] = images[i].shape
     shapes = np.ones((ncombine, ndim), dtype=int)
     raw_shapes = np.ones((ncombine, ndim), dtype=int)
@@ -527,7 +527,8 @@ def extract_stack_metadata(
             if extract_snoise:
                 sns[i] = float(hdr[snoise])
 
-            if hdr["NAXIS"] != ndim:
+            item_ndim = item.data.ndim if isinstance(item, CCDData) else hdr["NAXIS"]
+            if item_ndim != ndim:
                 raise ValueError(
                     "All FITS files must have the identical ndim, "
                     + "though they can have different sizes."
@@ -553,7 +554,11 @@ def extract_stack_metadata(
                 )
 
             # NOTE: the indexing in python is [z, y, x] order!!
-            raw_shape = tuple(int(hdr[f"NAXIS{i}"]) for i in range(ndim, 0, -1))
+            raw_shape = (
+                item.data.shape
+                if isinstance(item, CCDData)
+                else tuple(int(hdr[f"NAXIS{j}"]) for j in range(ndim, 0, -1))
+            )
             raw_shapes[i,] = raw_shape
             shapes[i,] = _trimmed_shape(raw_shape, trimsec)
         else:
@@ -571,6 +576,10 @@ def extract_stack_metadata(
                 )[0]
                 if data is None:
                     raise ValueError(f"Could not read data from input {i}.") from None
+                if data.ndim != ndim:
+                    raise ValueError(
+                        "All input images must have the identical ndim."
+                    ) from None
                 raw_shape = data.shape
             else:
                 _, hdr = _parse_imc_data_header(
diff --git a/src/astroimred/imutil/imstat.py b/src/astroimred/imutil/imstat.py
index ef2f292..4fca1a0 100644
--- a/src/astroimred/imutil/imstat.py
+++ b/src/astroimred/imutil/imstat.py
@@ -5,7 +5,7 @@
 
 import numpy as np
 import reducers.lowlevel as rdl
-from astro_ndslice import slicefy
+from astro_ndslice import listify, slicefy
 from astropy import units as u
 from astropy.io import fits
 from astropy.visualization import ZScaleInterval
@@ -200,17 +200,38 @@ def give_stats(
     num_extrema : 2-tuple of int or `None`, optional
         Number of low and high extreme values to report as ``(n_lo, n_hi)``. If
         `None`, extrema calculation is skipped.
+
+    Notes
+    -----
+    Sections may have different sizes. Overlapping sections count their shared
+    pixels once per section. Masks and non-finite pixels are excluded.
+
+    Raises
+    ------
+    ValueError
+        If no finite unmasked pixels remain in the selected sections.
     """
     data, hdr = _data_header_from_array_or_path(item, extension=extension)
     if mask is not None:
-        data = np.array(data, copy=True)
-        data[mask] = np.nan
+        mask = np.broadcast_to(np.asarray(mask, dtype=bool), data.shape)
 
     if statsecs is not None:
-        statsecs = [statsecs] if isinstance(statsecs, str) else list(statsecs)
-        data = np.array([data[slicefy(sec)] for sec in statsecs])
+        if isinstance(statsecs, tuple) and all(isinstance(s, slice) for s in statsecs):
+            statsecs = [statsecs]
+        else:
+            statsecs = listify(statsecs)
+        samples = []
+        for sec in statsecs:
+            section = slicefy(sec, ndim=data.ndim)
+            sample = data[section]
+            samples.append(sample.ravel() if mask is None else sample[~mask[section]])
+        data = np.concatenate(samples) if samples else np.empty(0, dtype=data.dtype)
+    elif mask is not None:
+        data = data[~mask]
 
     data = _finite_reducer_values(data)
+    if data.size == 0:
+        raise ValueError("No finite unmasked pixels remain in the selected sections.")
 
     std, mean = rdl.std_mean_valid(data, ddof=1)
     d_min, d_max = rdl.minmax_valid(data)
diff --git a/src/astroimred/phot/aperture.py b/src/astroimred/phot/aperture.py
index fdd6d45..ef9f50b 100644
--- a/src/astroimred/phot/aperture.py
+++ b/src/astroimred/phot/aperture.py
@@ -36,8 +36,8 @@ def cutout_from_ap(
 
     Parameters
     ----------
-    ap : `~astroapers.Aperture`
-        The aperture object.
+    ap : astroapers aperture object
+        The aperture object, such as `~astroapers.CircAp` or `~astroapers.EllipAp`.
     ccd : `~astropy.nddata.CCDData` or `numpy.ndarray`
         The CCD data.
     method : str, optional
@@ -48,10 +48,11 @@ def cutout_from_ap(
     Returns
     -------
     `~astropy.nddata.Cutout2D` or list of `~astropy.nddata.Cutout2D`
-        The cutout objects.
+        Cutouts trimmed to the image boundary for every method, with matching
+        data, coordinates, and WCS when supplied by `ccd`.
     """
     data = ccd.data if isinstance(ccd, CCDData) else np.asarray(ccd)
-    positions = np.asarray(ap.positions, dtype=np.float64).reshape(-1, 2)
+    wcs = getattr(ccd, "wcs", None) if isinstance(ccd, CCDData) else None
     if method not in {"bbox", "center", "exact"}:
         raise ValueError(f"Unsupported aperture method: {method!r}")
     boxes = ap.bboxes()
@@ -61,10 +62,11 @@ def cutout_from_ap(
     elif method == "exact":
         cutout_data = ap.weighted_cutout(data, fill_value=fill_value)
     cuts = []
-    for idx, (pos, box) in enumerate(zip(positions, boxes, strict=True)):
-        cut = Cutout2D(data, position=pos, size=box.shape)
+    for idx, (pos, box) in enumerate(zip(ap.positions, boxes, strict=True)):
+        cut = Cutout2D(data, position=pos, size=box.shape, mode="trim", wcs=wcs)
         if method != "bbox":
-            cut.data = cutout_data[idx]
+            sl_cut = box.overlap_slices(data.shape)[1]
+            cut.data = cutout_data[idx][sl_cut]
         cuts.append(cut)
     return cuts[0] if len(cuts) == 1 else cuts
 
diff --git a/src/astroimred/phot/center.py b/src/astroimred/phot/center.py
index c2ff2e3..5431411 100644
--- a/src/astroimred/phot/center.py
+++ b/src/astroimred/phot/center.py
@@ -2,6 +2,7 @@
 import logging
 
 import numpy as np
+import reducers as rd
 import reducers.lowlevel as rdl
 from astropy.modeling import Fittable2DModel, Parameter
 from astropy.modeling.fitting import LevMarLSQFitter
@@ -717,19 +718,19 @@ def _fit_2dgaussian(data, error=None, mask=None):
         )
 
     # assign zero weight to masked pixels
-    if data.mask is not np.ma.nomask:
-        weights[data.mask] = 0.0
+    mask = np.ma.getmaskarray(data)
+    weights[mask] = 0.0
 
-    mask = data.mask
-    data.fill_value = 0.0
-    data = data.filled()
+    data = data.filled(0.0)
 
     # Subtract the minimum of the data as a rough background estimate.
     # This will also make the data values positive, preventing issues with
     # the moment estimation in data_properties. Moments from negative data
     # values can yield undefined Gaussian parameters, e.g., x/y_stddev.
+    data_min, data_max = rd.minmax(data)
+
     props = sep_extract(
-        data - np.min(data),
+        data - data_min,
         thresh=0.0,  # Use all data points
         mask=mask,
         filter_kernel=None,  # No convolution
@@ -737,8 +738,13 @@ def _fit_2dgaussian(data, error=None, mask=None):
         clean=False,  # No cleaning
     )[0]
 
+    if len(props) == 0:
+        raise ValueError(
+            "No source detected in the cutout data to initialize 2D Gaussian fit."
+        )
+
     init_const = 0.0  # subtracted data minimum above
-    init_amplitude = np.ptp(data)
+    init_amplitude = data_max - data_min
     g_init = GaussianConst2D(
         constant=init_const,
         amplitude=init_amplitude,
diff --git a/src/astroimred/phot/radprof.py b/src/astroimred/phot/radprof.py
index 462c05c..9be0ade 100644
--- a/src/astroimred/phot/radprof.py
+++ b/src/astroimred/phot/radprof.py
@@ -6,10 +6,11 @@
 import astroapers.kernels as aapk
 import numpy as np
 import pandas as pd
+import reducers as rd
 from astropy.nddata import CCDData
 
 from .background import sky_fit
-from .center import circular_bbox_cut
+from .center import circular_bbox_cut, circular_slice
 
 __all__ = [
     "moffat_r",
@@ -367,6 +368,9 @@ def radprof_pix(img, pos, mask=None, rmax=10, sort_dist=False, fitfunc=None, ref
         The image to be profiled.
     pos : array_like
         The xy coordinates of the center of the object (0-indexing).
+    mask : array-like or bool, optional
+        Pixels to exclude. May match the full image or broadcast to the local
+        cutout. Pass ``mask=ccd.mask`` explicitly to use a CCDData mask.
     rmax : int, optional
         The maximum radius to be profiled.
         Default is ``10``.
@@ -392,7 +396,13 @@ def radprof_pix(img, pos, mask=None, rmax=10, sort_dist=False, fitfunc=None, ref
         raise TypeError(f"img must be a CCDData or ndarray (now {type(img) = })")
 
     cut, _, _, dists = circular_bbox_cut(img, pos, radius=rmax, return_dists=True)
-    mask = (dists > rmax) if mask is None else ((dists > rmax) | mask)
+    if mask is not None:
+        mask = np.asarray(mask, dtype=bool)
+        if mask.shape == img.shape:
+            mask = mask[circular_slice(img.shape, pos, rmax)]
+        mask = (dists > rmax) | np.broadcast_to(mask, cut.shape)
+    else:
+        mask = dists > rmax
 
     # cut = Cutout2D(img, pos, 2*rmax + 1).data
     # pos_cut = cut.
@@ -411,21 +421,22 @@ def radprof_pix(img, pos, mask=None, rmax=10, sort_dist=False, fitfunc=None, ref
 
         _r = dists[~mask]
         _i = cut[~mask]
-        _imin, _imax = _i.min(), _i.max()
+        _imin, _imax = rd.minmax(_i)
+        _rmin = rd.min(_r)
         if fitfunc == "gauss":
             fitter = gauss_r
-            p0 = [_i[_r == _r.min()][0] - _imin, max(1, rmax / 6), _imin]
+            p0 = [_i[_r == _rmin][0] - _imin, max(1, rmax / 6), _imin]
             bounds = np.array([[0, 0, _imin - 1], [np.inf, rmax / 2, _imax + 1]])
             # assuming the user gave ~ (2-3)x FWHM, and we want sigma ~ 0.4xFWHM
         elif fitfunc == "moffat":
             fitter = moffat_r
-            p0 = [_i[_r == _r.min()][0] - _imin, 1, 2.5, _imin]
+            p0 = [_i[_r == _rmin][0] - _imin, 1, 2.5, _imin]
             bounds = np.array(
                 [[0, 0.1, 1.0, _imin - 1], [np.inf, rmax, np.inf, _imax + 1]]
             )
         elif fitfunc == "bivt":
             fitter = bivt_r
-            p0 = [_i[_r == _r.min()][0] - _imin, 1, max(1, rmax / 6), _imin]
+            p0 = [_i[_r == _rmin][0] - _imin, 1, max(1, rmax / 6), _imin]
             bounds = np.array(
                 [[0, -2, 1.0e-10, _imin - 1], [np.inf, 15, 2 * rmax, _imax + 1]]
             )
diff --git a/tests/imutil/test_imstat.py b/tests/imutil/test_imstat.py
index d8c339a..8e310ba 100644
--- a/tests/imutil/test_imstat.py
+++ b/tests/imutil/test_imstat.py
@@ -168,6 +168,53 @@ def test_stats_mask_does_not_mutate_input(self):
 
         np.testing.assert_array_equal(arr, original)
 
+    def test_integer_mask_and_unequal_sections(self) -> None:
+        data = np.arange(12, dtype=np.int16).reshape(3, 4)
+        original = data.copy()
+        mask = np.zeros(data.shape, dtype=bool)
+        mask[0, 1] = True
+        result = imstat.give_stats(data, mask=mask, statsecs=["[1:2,1:1]", "[3:4,2:3]"])
+        # Selected pixels: 0, 6, 7, 10, 11 (pixel 1 is masked).
+        assert result["num"] == 5
+        assert result["avg"] == pytest.approx(6.8)
+        assert result["med"] == 7
+        assert result["std"] == pytest.approx(np.sqrt(18.7))
+        np.testing.assert_array_equal(data, original)
+        assert mask[0, 1]
+        assert np.count_nonzero(mask) == 1
+
+    def test_unequal_sections(self) -> None:
+        result = imstat.give_stats(
+            np.arange(12.0).reshape(3, 4),
+            statsecs=["[1:2,1:1]", "[3:4,2:3]"],
+        )
+        assert result["num"] == 6
+        assert result["avg"] == pytest.approx(35 / 6)
+
+    @pytest.mark.parametrize("section", [slice(1, 3), (slice(1, 3), slice(0, 2))])
+    def test_single_slice_section(self, section: object) -> None:
+        data = np.arange(12.0).reshape(3, 4)
+        expected = (
+            data[section, section] if isinstance(section, slice) else data[section]
+        )
+        result = imstat.give_stats(data, statsecs=section)
+        assert result["num"] == expected.size
+        assert result["avg"] == pytest.approx(np.mean(expected))
+
+    @pytest.mark.parametrize(
+        "data, kwargs",
+        [
+            (np.array([], dtype=float), {}),
+            (np.array([np.nan, np.inf]), {}),
+            (np.array([1, 2]), {"mask": True}),
+            (np.ones((2, 2)), {"statsecs": []}),
+            (np.ones((2, 2)), {"statsecs": "[5:6,5:6]"}),
+        ],
+    )
+    def test_no_finite_selected_pixels(self, data: np.ndarray, kwargs: dict) -> None:
+        with pytest.raises(ValueError, match="No finite unmasked pixels"):
+            imstat.give_stats(data, **kwargs)
+
     def test_stats_path_input(self, temp_fits_file):
         """Test statistics on a path-like FITS input."""
         result = imstat.give_stats(temp_fits_file)
diff --git a/tests/phot/test_aperture.py b/tests/phot/test_aperture.py
index 1754f9e..423ecf4 100644
--- a/tests/phot/test_aperture.py
+++ b/tests/phot/test_aperture.py
@@ -121,3 +121,58 @@ def test_cutout_from_ap_rejects_unknown_method():
 
     with pytest.raises(ValueError):
         cutout_from_ap(ap, data, method="not-a-method")
+
+
+def test_cutout_from_ap_edge_coordinates_and_shape():
+    """Cutouts crossing image boundaries must maintain trimmed shape and native coordinates."""
+    from astropy.wcs import WCS
+
+    wcs = WCS(naxis=2)
+    wcs.wcs.crpix = [50, 50]
+    wcs.wcs.cdelt = [0.1, 0.1]
+    wcs.wcs.crval = [10.0, 20.0]
+    wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"]
+
+    data = CCDData(np.ones((100, 100)), wcs=wcs, unit="adu")
+    ap = aap.EllipAp((5.3, 5.7), a=10.2, b=4.6, theta=0.3)
+    pos = ap.positions[0]
+    box = ap.bboxes()[0]
+    sl_img, _ = box.overlap_slices(data.shape)
+    expected_shape = (
+        sl_img[0].stop - sl_img[0].start,
+        sl_img[1].stop - sl_img[1].start,
+    )
+
+    for method in ("bbox", "center", "exact"):
+        cut = cutout_from_ap(ap, data, method=method, fill_value=np.nan)
+        assert cut.shape == expected_shape
+        assert cut.data.shape == expected_shape
+        cut_pos = cut.to_cutout_position(pos)
+        assert_allclose(cut.input_position_original, pos)
+        assert_allclose(cut.input_position_cutout, cut_pos)
+        orig_pos = cut.to_original_position(cut_pos)
+        assert_allclose(orig_pos, pos, atol=1e-12)
+        assert_allclose(cut_pos, (pos[0] - sl_img[1].start, pos[1] - sl_img[0].start))
+        assert_allclose(
+            cut.wcs.pixel_to_world_values(*cut_pos),
+            wcs.pixel_to_world_values(*pos),
+            rtol=0,
+            atol=1e-10,
+        )
+        if method == "bbox":
+            assert_allclose(cut.data, data.data[sl_img])
+        else:
+            weights = ap.weights_center() if method == "center" else ap.weights_exact()
+            local_weights = weights[0][box.overlap_slices(data.shape)[1]]
+            np.testing.assert_allclose(
+                cut.data[local_weights > 0], local_weights[local_weights > 0]
+            )
+
+
+def test_cutout_from_ap_integer_preserves_dtype():
+    """Integer images with bbox method should preserve their integer dtype without padding error."""
+    int_data = np.arange(100, dtype=np.int32).reshape(10, 10)
+    ap = aap.CircAp((0.5, 0.5), r=3)
+    cut = cutout_from_ap(ap, int_data, method="bbox")
+    assert cut.data.dtype == np.int32
+    assert cut.shape == cut.data.shape
diff --git a/tests/phot/test_radprof.py b/tests/phot/test_radprof.py
index ff638b2..d4c7a6b 100644
--- a/tests/phot/test_radprof.py
+++ b/tests/phot/test_radprof.py
@@ -6,6 +6,7 @@
 
 import astroapers as aap
 import numpy as np
+import pytest
 from numpy.testing import assert_allclose
 
 from astroimred.phot.radprof import (
@@ -389,3 +390,63 @@ def test_radprof_pix_at_edge(self, uniform_100x100):
         # Should still return values
         assert len(vals) > 0
         assert_allclose(vals, 10.0, rtol=1e-10)
+
+    def test_radprof_pix_with_full_frame_mask(self):
+        """radprof_pix must accept full-frame masks and exclude masked pixels."""
+        from astropy.nddata import CCDData
+
+        img = np.ones((100, 100), dtype=float)
+        mask = np.zeros((100, 100), dtype=bool)
+        mask[51, 50] = True
+
+        r, vals = radprof_pix(img, pos=(50, 50), mask=mask, rmax=10)
+        assert len(r) > 0
+        assert_allclose(vals, 1.0)
+
+        ccd = CCDData(img, mask=mask, unit="adu")
+        r_ccd, vals_ccd = radprof_pix(ccd, pos=(50, 50), mask=ccd.mask, rmax=10)
+        assert len(r_ccd) == len(r)
+        assert_allclose(r_ccd, r)
+        assert_allclose(vals_ccd, vals)
+
+        r_edge, vals_edge = radprof_pix(img, pos=(5, 5), mask=mask, rmax=10)
+        assert len(r_edge) > 0
+        assert_allclose(vals_edge, 1.0)
+
+    def test_radprof_pix_fit_with_full_frame_mask(self):
+        """radprof_pix with fitfunc must run with full-frame mask."""
+        img = np.ones((100, 100), dtype=float)
+        mask = np.zeros((100, 100), dtype=bool)
+        mask[51, 50] = True
+
+        r, vals, fitter, popt, fwhm = radprof_pix(
+            img, pos=(50, 50), mask=mask, rmax=10, fitfunc="gauss"
+        )
+        assert len(r) > 0
+        assert np.isfinite(fwhm)
+
+    def test_radprof_pix_with_scalar_mask(self):
+        """radprof_pix must accept scalar boolean masks."""
+        img = np.ones((50, 50), dtype=float)
+        r_unmasked, vals_unmasked = radprof_pix(img, pos=(25, 25), mask=False, rmax=5)
+        assert len(r_unmasked) > 0
+        assert_allclose(vals_unmasked, 1.0)
+
+        r_masked, vals_masked = radprof_pix(img, pos=(25, 25), mask=True, rmax=5)
+        assert len(r_masked) == 0
+
+    @pytest.mark.parametrize("shape", [(5,), (1, 5), (5, 1)])
+    def test_radprof_pix_broadcast_cutout_mask(self, shape: tuple[int, ...]) -> None:
+        mask = np.array([False, False, True, False, False]).reshape(shape)
+        radii, values = radprof_pix(np.ones((11, 11)), (5, 5), mask=mask, rmax=2)
+        assert radii.size == 8
+        assert_allclose(values, 1)
+
+    def test_radprof_pix_mask_is_explicit_for_ccddata(self) -> None:
+        from astropy.nddata import CCDData
+
+        ccd = CCDData(np.ones((11, 11)), mask=np.ones((11, 11), dtype=bool), unit="adu")
+        radii, _ = radprof_pix(ccd, (5, 5), rmax=2)
+        assert radii.size == 13
+        masked_radii, _ = radprof_pix(ccd, (5, 5), mask=ccd.mask, rmax=2)
+        assert masked_radii.size == 0