How to get all spherical angles $\phi$ and $\theta$, which describe a circle on a sphere?

I have a spherical rendering, where the spherical coordinates $\phi$ and $\theta$ are represented by the x and y axis of the image (similar to how world maps work):

Now given a point on the image with the pixel coordinates p($\phi$, $\theta$) I want to calculate the average color of all the neighboring pixels on the sphere within a radius (radius could maybe more easily be defined as an angle, since we are thinking of a unit sphere). For this I need to sum up all of the pixels within the circle, however the circle on the sphere will not be a circle on my map. So the question is, how can I get all the pixel coordinates (= my $\phi$ and $\theta$ angles) within the circle?

Edit: Forgot to describe, how the rendering was done. With $\phi = [0, 2\pi]$ we divide $2\pi$ by the horizontal image resolution. $\theta = [0, \pi]$, so we divide $\pi$ by the vertical image resolution. Using a horizontal resolution that is twice the vertical resolution one pixel represents area on the sphere that is $(\frac{2 * \pi}{res_h}\times\frac{2 * \pi}{res_h})$.

I will stick to your notations $\phi$ for longitude, $\theta$ for latitude, although one finds rather often the inverse convention.

Given: a center $M_0$ (defined by spherical coordinates $(\phi_0,\theta_0)$) and a radius $R$ ($0<R<\pi/2$, measured as an arc on the unit sphere).

The set of points $M$ of the sphere which are interior to the spherical disk with center $M_0$ and radius $R$ is given by the following dot product constraint:

$$\vec{OM}.\vec{OM_0}>\cos(R)$$

which is equivalent, using classical spherical coordinates, to:

$$\begin{pmatrix}\cos(\phi)\cos(\theta)\\ \sin(\phi)\cos(\theta)\\ \sin(\theta)\end{pmatrix} .\begin{pmatrix}\cos(\phi_0)\cos(\theta_0)\\ \sin(\phi_0)\cos(\theta_0)\\ \sin(\theta_0)\end{pmatrix}>\cos(R)\tag{1}$$

This constraint can be written under the form:

$$\cos(\theta_0)\cos(\theta)\cos(\phi-\phi_0)+\sin(\theta_0)\sin(\theta)>\cos(R)\tag{2}$$

which is an implicit equation in $(\phi,\theta)$ depending upon three parameters $(\phi_0,\theta_0,R)$ that can be visualized with this Geogebra animation (play with the sliders !):

https://www.geogebra.org/calculator/eanc7njp

The top sliders $f$ and $t$ refer to the coordinates $\phi_0$ and $\theta_0$ resp. of center $M_0$.

Here are two examples:

Fig. 1: An "ordinary" circle centered in (0,0) with radius $\pi/4$ rendered as a kind of ellipse.

Fig. 2: A "limit case" image of a circle belonging to northern hemisphere, tangent to the equator, (almost) passing throughout North Pole, explaining the almost linear segment ranging approximately from $-\pi/2$ to $\pi/2$.

Remark 1: constraint (1) has been given in the same form by @blamocur.

Remark 2: formula (2) could have been obtained directly by using the spherical law of cosines..

Remark 3: Explicit equations of the circle can be found here.

Let $a$ be the axis of the circle passing through its center, and let $\theta_0$ be the angle of the radius (i.e. $\text{radius} = \sin \theta_0$ )

The $a$ axis has coordinates:

$a = ( \sin \theta_a \cos \phi_a, \sin \theta_a \sin \phi_a, \cos \theta_a)$

you have to find two orthogonal vectors that are orthogonal to $a$, and these can be chosen as follows:

$u_1 = ( \cos \theta_a \cos \phi_a, \cos \theta_a \sin \phi_a , - \sin \theta_a )$

$u_2 = (- \sin \phi_a, \cos \phi_a , 0 )$

The rectangular coordinates of points at the circular region boundary are given by (with respect to the basis $[u_1, u_2, a]$, utilizing the respective spherical coordinates) are:

$(x', y', z') = ( \sin \theta_0 \cos \phi, \sin \theta_0 \sin \phi, \cos \theta_0 )$

where $0 \le \phi \le 2 \pi $

now the world rectangular coordinates of the circular region boundary is

$(x, y, z) = [u_1, u_2, a] (x', y', z') = x' u_1 + y' u_2 + z' a$

once you've calculated $(x, y, z)$ , you can calculate the corresponding $\theta, \phi$ in the original world coordinates, as follows

$ \theta = \cos^{-1} (z)$ and $ \phi = \text{atan2} (x, y) $

By changing $\phi$ over $[0, 2 \pi]$ you get a closed curve in your image, that encloses all the pixels you want.