diff --git a/docs/analyze_radial_percentile.md b/docs/analyze_radial_percentile.md
new file mode 100644
index 000000000..dcb516b8d
--- /dev/null
+++ b/docs/analyze_radial_percentile.md
@@ -0,0 +1,48 @@
+## Analyze radial percentile
+
+This function calculates the average value of pixels within a cutoff threshold from the center of an object and writes 
+the values out to the [Outputs class](outputs.md).  
+
+
+**plantcv.analyze.radial_percentile**(*img, mask, roi=None, percentile=50, label=None*)
+
+**returns** List of average values for either grayscale or RGB
+
+- **Parameters:**
+    - img - RGB or grayscale image.
+    - mask - Binary mask.
+    - roi - Optional ROIs for calculating on multiple objects in an image.
+    - percentile - cutoff for considering pixels in the average. Expressed as a percent of maximum distance from the object's center (default = 50).
+    - label - Optional label parameter, modifies the variable name of observations recorded. Can be a prefix or list (default = pcv.params.sample_label).
+- **Outputs:**
+    - A list of average values for RGB or gray channels for each object (if ROIs are provided), or for the single object in the image. 
+- **Example use:**
+    - Useful for calculating the intensity of the middle of seeds from an X-ray image.
+    - Also could be useful in determining if there are color differences in the middle of a plant rosette. 
+
+- **Output data stored:** Data ('gray_X%_avg', or 'red_X%_avg', 'green_X%_avg', 'blue_X%_avg') automatically gets stored to
+the [`Outputs` class](outputs.md) when this function is ran. These data can always get accessed during a workflow (example
+below). For more detail about data output see [Summary of Output Observations](output_measurements.md#summary-of-output-observations)
+
+**Multi-ROI object on original image**
+
+![Screenshot](img/documentation_images/analyze_radial/radial_doc1.png)
+
+```python
+
+from plantcv import plantcv as pcv
+
+# Caclulates the average values of pixels that fall within the distance percentile from the center of an object.
+list_of_averages  = pcv.analyze.radial_percentile(img=img, mask=mask, roi=rois, percentile=40)
+
+# Access data stored out from analyze.radial_percentile
+gray_avg_seed1 = pcv.outputs.observations['default_1']['gray_40%_avg']['value']
+
+```
+
+
+**Debug depicting one example of seed cropped to 40% of maximum distance from center**
+
+![Screenshot](img/documentation_images/analyze_radial/radial_doc2.png)
+
+**Source Code:** [Here](https://github.com/danforthcenter/plantcv/blob/main/plantcv/plantcv/analyze/radial.py)
diff --git a/docs/img/documentation_images/analyze_radial/radial_doc1.png b/docs/img/documentation_images/analyze_radial/radial_doc1.png
new file mode 100644
index 000000000..90a367b31
Binary files /dev/null and b/docs/img/documentation_images/analyze_radial/radial_doc1.png differ
diff --git a/docs/img/documentation_images/analyze_radial/radial_doc2.png b/docs/img/documentation_images/analyze_radial/radial_doc2.png
new file mode 100644
index 000000000..3cb96089e
Binary files /dev/null and b/docs/img/documentation_images/analyze_radial/radial_doc2.png differ
diff --git a/docs/updating.md b/docs/updating.md
index ccb9203ff..dff3919c1 100644
--- a/docs/updating.md
+++ b/docs/updating.md
@@ -367,6 +367,10 @@ pages for more details on the input and output variable types.
 * pre v4.0: NA
 * post v4.0: npq, npq_hist = **plantcv.analyze.npq**(*ps_da_light, ps_da_dark, labeled_mask, n_labels=1, auto_fm=False, min_bin=0, max_bin="auto", measurement_labels=None, label=None*)
 
+#### plantcv.analyze.radial_percentile
+* pre v4.9: NA
+* post v9.9: avgs = **plantcv.analyze.radial_percentile**(*img, mask, roi=None, percentile=50, label="default"*)
+
 #### plantcv.analyze.size
 
 * pre v4.0: (see plantcv.analyze_object)
diff --git a/mkdocs.yml b/mkdocs.yml
index db666df16..a88ebd0fb 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -40,6 +40,7 @@ nav:
         - 'Analyze Spectral Index': analyze_spectral_index.md
         - 'Analyze YII': analyze_yii.md
         - 'Analyze NPQ': analyze_npq.md
+        - 'Analyze Radial Average': analyze_radial_percentile.md
       - 'Annotation Tools':
           - 'Points': Points.md
       - 'Apply Mask': apply_mask.md
diff --git a/plantcv/plantcv/analyze/__init__.py b/plantcv/plantcv/analyze/__init__.py
index 414d1e3d9..e001d80ed 100644
--- a/plantcv/plantcv/analyze/__init__.py
+++ b/plantcv/plantcv/analyze/__init__.py
@@ -9,6 +9,7 @@
 from plantcv.plantcv.analyze.yii import yii
 from plantcv.plantcv.analyze.npq import npq
 from plantcv.plantcv.analyze.distribution import distribution
+from plantcv.plantcv.analyze.radial import radial_percentile
 
 __all__ = ["color", "bound_horizontal", "bound_vertical", "grayscale", "size", "thermal", "spectral_reflectance",
-           "spectral_index", "yii", "npq", "distribution"]
+           "spectral_index", "yii", "npq", "distribution", "radial_percentile"]
diff --git a/plantcv/plantcv/analyze/radial.py b/plantcv/plantcv/analyze/radial.py
new file mode 100644
index 000000000..2be125630
--- /dev/null
+++ b/plantcv/plantcv/analyze/radial.py
@@ -0,0 +1,146 @@
+"""Outputs the average pixel values from a percentile radially outward from an object's center."""
+from plantcv.plantcv import params, Objects, apply_mask, auto_crop, outputs
+from plantcv.plantcv import roi as roi_
+from plantcv.plantcv._debug import _debug
+import cv2
+import os
+import numpy as np
+
+
+def _calc_dists(img, mask, percentile, store_debug, ind=None):
+    """Calculates the distances of each pixel,
+    the cutoff based on given percentile,
+    and the average pixel values within that cutoff.
+
+    Parameters
+    ----------
+    img : numpy.ndarray
+        RGB or grayscale image cropped to a focal object
+    mask : numpy.ndarray
+        Segmented mask of the object from img
+    percentile : number
+        Cutoff for inclusion of pixels (percent from center)
+    store_debug : str or None
+        Value for params.debug when function is called
+    ind : int, optional
+        Tracker for which in multi ROI gets a debug, by default None
+
+    Returns
+    -------
+    avgs : float or list
+        average pixel values (gray or RGB) within the distance percentile
+    """
+    # if mask is empty then return NaNs, nothing is in the mask.
+    if np.sum(mask) == 0:
+        return ["nan", "nan", "nan"]
+    # Analyze shape properties
+    m = cv2.moments(mask, binaryImage=True)
+    cmx = m['m10'] / m['m00']
+    cmy = m['m01'] / m['m00']
+    center = (cmx, cmy)
+    # Calculate point distances from center
+    y, x = np.ogrid[:img.shape[0], :img.shape[1]]
+    distances = np.sqrt((x - center[0])**2 + (y - center[1])**2)
+    if len(img.shape) == 3:
+        distances = np.stack([distances for _ in range(3)], axis=2)
+    # Distance cutoff based on percentile
+    cutoff = (np.sqrt((img.shape[0]/2)**2 + (img.shape[1]/2)**2))*(percentile/100)
+    img_cutoff = np.where(distances < cutoff, img, np.nan)
+    # One example debug
+    if ind == 0:
+        example = np.where(distances < cutoff, img, 0)
+        params.debug = store_debug
+        _debug(visual=example, filename=os.path.join(params.debug_outdir, str(params.device) + "_radial_average.png"))
+        params.debug = None
+
+    avgs = float(np.nanmean(img_cutoff))
+
+    if len(img.shape) == 3:
+        avgs = [float(np.nanmean(img_cutoff[:, :, [2]])),
+                float(np.nanmean(img_cutoff[:, :, [1]])),
+                float(np.nanmean(img_cutoff[:, :, [0]]))]
+    return avgs
+
+
+def radial_percentile(img, mask, roi=None, percentile=50, label=None):
+    """_summary_
+
+    Parameters
+    ----------
+    img : numpy.ndarray
+        RGB or grayscale image
+    mask : numpy.ndarray
+        Binary mask with objects of interest segmented
+    roi : plantcv.plantcv.Objects, optional
+        Region of Interest, single or multi, to identify objects, by default None
+    percentile : int, optional
+        Percentile of max distance from center in which to average pixel values, by default 50
+    label : str, optional
+        Optional label for outputs (default = pcv.params.sample_label)
+
+    Returns
+    -------
+    avgs : list
+        average pixel values (gray or RGB) within the distance percentile
+    """
+    if label is None:
+        label = params.sample_label
+
+    store_debug = params.debug
+    params.debug = None
+    if roi:
+        avgs = []
+        for i, _ in enumerate(roi.contours):
+            # Loop through rois (even if there is only 1)
+            roi_ind = Objects(contours=[roi.contours[i]], hierarchy=[roi.hierarchy[i]])
+            # Filter mask by roi and apply it to the image
+            filt = roi_.filter(mask=mask, roi=roi_ind)
+            # Check for empty
+            if len(np.unique(filt)) == 1:
+                noavg = ["nan"]
+                if len(img.shape) == 3:
+                    noavg = ["nan", "nan", "nan"]
+                avgs.append(noavg)
+            else:
+                masked = apply_mask(img=img, mask=filt, mask_color='black')
+                # Crop the image and the mask to the roi
+                crop_img = auto_crop(img=masked, mask=filt, padding_x=1, padding_y=1, color='black')
+                crop_mask = auto_crop(img=filt, mask=filt, padding_x=1, padding_y=1, color='black')
+
+                # Calculate average of each channel
+                avgs.append(_calc_dists(img=crop_img, mask=crop_mask, percentile=percentile, store_debug=store_debug, ind=i))
+
+    else:
+        masked = apply_mask(img=img, mask=mask, mask_color='black')
+        # Crop the image to the mask if the mask is not empty
+        crop_img = masked
+        crop_mask = mask
+        if np.sum(mask) > 0:
+            crop_img = auto_crop(img=masked, mask=mask, padding_x=1, padding_y=1, color='black')
+            crop_mask = auto_crop(img=mask, mask=mask, padding_x=1, padding_y=1, color='black')
+        # Calculate averages of each channel
+        avgs = [_calc_dists(img=crop_img, mask=crop_mask, percentile=percentile, store_debug=store_debug, ind=0)]
+
+    # Outputs
+    for idx, i in enumerate(avgs):
+        if isinstance(i, float):
+            outputs.add_observation(sample=label+"_"+str(idx+1), variable='gray_'+str(percentile)+'%_avg',
+                                    trait='gray_'+str(percentile)+'%_radial_average',
+                                    method='plantcv.plantcv.analyze.radial',
+                                    scale='none', datatype=float, value=i, label='none')
+        elif isinstance(i, list):
+            outputs.add_observation(sample=label+"_"+str(idx+1), variable='red_'+str(percentile)+'%_avg',
+                                    trait='red_'+str(percentile)+'%_radial_average',
+                                    method='plantcv.plantcv.analyze.radial',
+                                    scale='none', datatype=float, value=i[0], label='none')
+            outputs.add_observation(sample=label+"_"+str(idx+1), variable='green_'+str(percentile)+'%_avg',
+                                    trait='green_'+str(percentile)+'%_radial_average',
+                                    method='plantcv.plantcv.analyze.radial',
+                                    scale='none', datatype=float, value=i[1], label='none')
+            outputs.add_observation(sample=label+"_"+str(idx+1), variable='blue_'+str(percentile)+'%_avg',
+                                    trait='blue_'+str(percentile)+'%_radial_average',
+                                    method='plantcv.plantcv.analyze.radial',
+                                    scale='none', datatype=float, value=i[2], label='none')
+
+    params.debug = store_debug
+    return avgs
diff --git a/tests/conftest.py b/tests/conftest.py
index 451d854c3..897d22bd6 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -89,6 +89,11 @@ def __init__(self):
         self.kmeans_classifier_gray_dir = os.path.join(self.datadir, "kmeans_classifier_gray_dir")
         # nd2 file
         self.nd2_img = os.path.join(self.datadir, "test_nd2_img.nd2")
+        # Image, mask, and roi for radial percentile
+        self.rgb_seed = os.path.join(self.datadir, "rgb_seed.jpg")
+        self.rgb_seed_mask = os.path.join(self.datadir, "rgb_seed_mask.jpg")
+        self.small_circle = os.path.join(self.datadir, "small_circle.png")
+        self.empty_mask = os.path.join(self.datadir, "seed_empty_mask.png")
         # flir thermal img
         self.flir_img = os.path.join(self.datadir, "FLIR_test.jpg")
         # parallel results to process
diff --git a/tests/plantcv/analyze/test_radial_percentile.py b/tests/plantcv/analyze/test_radial_percentile.py
new file mode 100644
index 000000000..350f671e7
--- /dev/null
+++ b/tests/plantcv/analyze/test_radial_percentile.py
@@ -0,0 +1,49 @@
+"""Tests for pcv.analyze.radial"""
+import cv2
+import numpy as np
+from plantcv.plantcv.analyze.radial import radial_percentile
+from plantcv.plantcv import Objects
+from plantcv.plantcv._helpers import _cv2_findcontours
+
+
+def test_radial_RGB(test_data):
+    """Test for PlantCV."""
+    # Read in test data
+    img = cv2.imread(test_data.rgb_seed)
+    mask = cv2.imread(test_data.rgb_seed_mask, -1)
+    avgs = radial_percentile(img=img, mask=mask)
+    assert int(avgs[0][0]) == 138
+
+
+def test_radial_gray(test_data):
+    """Test for PlantCV."""
+    # Read in roi
+    circle = cv2.imread(test_data.small_circle, -1)
+    roi_contour, roi_hierarchy = _cv2_findcontours(bin_img=circle)
+    roi = Objects(contours=[roi_contour], hierarchy=[roi_hierarchy])
+    # Read in test data as gray
+    img = cv2.imread(test_data.rgb_seed, 0)
+    mask = cv2.imread(test_data.rgb_seed_mask, -1)
+    avgs = radial_percentile(img=img, mask=mask, roi=roi)
+    assert int(avgs[0]) == 96
+
+
+def test_radial_empty(test_data):
+    """Test for PlantCV."""
+    # Read in test data
+    img = cv2.imread(test_data.rgb_seed)
+    mask = mask = np.zeros((100,100),dtype="uint8")
+    avgs = radial_percentile(img=img, mask=mask)
+    assert avgs[0][0] == "nan"
+
+
+def test_radial_roi_empty(test_data):
+    """Test for PlantCV."""
+    # Read in test data
+    img = cv2.imread(test_data.rgb_seed)
+    mask = mask = np.zeros((100,100),dtype="uint8")
+    circle = cv2.imread(test_data.small_circle, -1)
+    roi_contour, roi_hierarchy = _cv2_findcontours(bin_img=circle)
+    roi = Objects(contours=[roi_contour], hierarchy=[roi_hierarchy])
+    avgs = radial_percentile(img=img, mask=mask, roi=roi)
+    assert avgs[0][0] == "nan"
diff --git a/tests/testdata/gray_seed.jpg b/tests/testdata/gray_seed.jpg
new file mode 100644
index 000000000..435cb1555
Binary files /dev/null and b/tests/testdata/gray_seed.jpg differ
diff --git a/tests/testdata/gray_seed_mask.jpg b/tests/testdata/gray_seed_mask.jpg
new file mode 100644
index 000000000..4e756589a
Binary files /dev/null and b/tests/testdata/gray_seed_mask.jpg differ
diff --git a/tests/testdata/rgb_seed.jpg b/tests/testdata/rgb_seed.jpg
new file mode 100644
index 000000000..73355110f
Binary files /dev/null and b/tests/testdata/rgb_seed.jpg differ
diff --git a/tests/testdata/rgb_seed_mask.jpg b/tests/testdata/rgb_seed_mask.jpg
new file mode 100644
index 000000000..f323aeee0
Binary files /dev/null and b/tests/testdata/rgb_seed_mask.jpg differ
diff --git a/tests/testdata/small_circle.png b/tests/testdata/small_circle.png
new file mode 100644
index 000000000..41d74bb1d
Binary files /dev/null and b/tests/testdata/small_circle.png differ