Learning Point Processes and Convolutional Neural Networks for Object Detection in Satellite Images

Abstract

Convolutional neural networks (CNN) have shown great results for object-detection tasks by learning texture and pattern-extraction filters. However, object-level interactions are harder to grasp without increasing the complexity of the architectures. On the other hand, Point Process models propose to solve the detection of the configuration of objects as a whole, allowing the factoring in of the image data and the objects’ prior interactions. In this paper, we propose combining the information extracted by a CNN with priors on objects within a Markov Marked Point Process framework. We also propose a method to learn the parameters of this Energy-Based Model. We apply this model to the detection of small vehicles in optical satellite imagery, where the image information needs to be complemented with object interaction priors because of noise and small object sizes.

1. Introduction

While object detection has been thoroughly studied for the last 20 years [1], the detection of small objects in optical satellite images remains a challenging task: With limited spatial resolution (around 0.5 m/pixel), objects such as cars can be only a few pixels wide. Moreover, this tiny-object scattering density varies greatly, and when closely packed, instances might be difficult to differentiate. However, as the geometrical configuration of one object is often dependent on its neighbors, we can leverage these priors to improve the detection results.

Methods based on convolutional neural networks (CNN) such as Faster R-CNN [2], YOLO [3], or RetinaNet [4] propose to detect objects in “natural” images. Some of these approaches first extract features through convolutions, then propose a series of boxes (anchors) that are refined by regression afterward [5,6]. The specificities of remote-sensing data (high number of small objects) motivate the use of anchor-free methods, relying on heatmap (probability of object center) inference for oriented vehicle detection [7,8,9], or ship detection [10] (here approximating the heatmap through a vector field).

In remote sensing, the sensor is intrinsically quite distant from the observed objects, implying a limited spatial resolution (usually measured in meters per pixel). This limited resolution, in turn, induces a limited amount of visual information to perform object detection. Given the sensor noise or atmospheric perturbations, remote-sensing detection methods need to be innovative to compensate for the limited signal. Approaches such as [11,12] use the temporal information from image time series to improve the detection of small objects. However, these are not suitable for detecting small static objects. To extract object orientation, Ref. [13] focuses on extracting fine-grained features while ignoring any inter-instance interactions. The interaction between nearby objects is key to analyzing a scene for human operators. However, those are often ignored in CNN-based approaches (at least not beyond a post-processing step such as non-maximum suppression) since CNN models are—by construction—most efficient at extracting texture and local information [14]. Meanwhile, long-range relationships are only caught when drastically increasing the depth of the network or by introducing dense connections. For instance, Ref. [15] models interactions through a cascade of Transformers at a great complexity cost, while [16] uses text-modal descriptors to introduce prior knowledge about the objects and their relations into the model.

Meanwhile, Point Processes [17] have been used to perform object detection while relying on a stochastic geometry model to jointly solve the detection and selection of objects with priors. Point Process approaches model the probability on the whole set of points in the image, thus taking into account the interactions between objects. These approaches rely on decomposing the Point Process density into a data term that measures the correspondence of the points against the image and some prior terms that measure the coherence of the point configuration itself. These approaches have shown good results for detection tasks in images with many small objects, such as biological imagery [18] or remote sensing for road [19,20] and building [21] extraction, or tree detection [22].

Previous methods make use of data terms built upon image contrast measures between the interior and exterior of the object [18,20,23] or image gradient [24], which perform great in images with clear contrast between objects and their background. However, in the case of optical satellite images, background diversity, variable visual aspects of objects, and inconsistent illumination make these contrast measures unreliable.

In this paper, we propose the incorporation of interaction models into object-detection methods while taking advantage of the capabilities of deep convolutional neural networks. Satellite images are inherently imperfect due to atmospheric disturbances, partial occlusions, and limited spatial resolution. To compensate for this lack of visual information, it becomes essential to incorporate prior knowledge about the layout of objects of interest.

Instead of increasing the complexity of the model by adding, for example, Transformers to a deep CNN architecture, we propose in this paper to combine CNN pattern extraction with the Point Process approach. The starting point of this approach is to use the output of a CNN as the data term for a Point Process model (Section 2.3). From the latter, we derive more efficient sampling methods for the Point Process that do not rely on application-specific heuristics (Section 2.4). Then, we propose to bridge the gap in terms of parameter estimation using modern learning techniques inspired by Energy-Based Models (Section 2.5). Ultimately, we introduce a novel scoring function that takes into account object interaction in measuring the confidence in the detections and enables explainable results (Section 2.6). The final part presents quantitative and qualitative results on optical satellite data (Section 3).

2. Materials and Methods

In this paper, we propose an object-detection method leveraging CNN-extracted information within a Point Process framework that models object interactions. First, we introduce the foundations for the Point Process on which our model is built, then detail the energy model, along with its sampling and parameter estimation method. We also propose a novel scoring method to evaluate the confidence in the output while considering the interactions and providing explainable results.

2.1. Materials: Datasets

Our application goal with Airbus Defense and Space (ADS) is the detection of small objects in images from satellites such as Pléiades, (Constellation Pléiades [https://pleiades.cnes.fr/en/PLEIADES/GP_systeme.htm, pleiades.cnes.fr, (accessed on 20 October 2023)]) or CO3D, (Constellation Optique 3D [https://www.eoportal.org/satellite-missions/co3d-constellation, eoportal.org, (accessed on 20 October 2023)]), which have a typical spatial resolution (or ground sampling distance) of $0.5$ m per pixel.

To train the CNN and infer the model parameters $\theta$, we compile ${\mathit{DOTA}}_{0.5}$, a dataset of vehicles in remote-sensing images with $0.5$ m per pixel resolution, from the DOTA dataset [25]. This dataset ${\mathit{DOTA}}_{0.5}$ is built by sub-sampling images from DOTA to the desired spatial resolution and keeping only small-vehicle and large-vehicle classes. The data consists of images from satellite or aerial sources of urban or countryside scenes with variable densities of objects (from isolated cars on country roads to dense parking lots). The images are labeled with sets of oriented rectangles for the objects of interest. We also compile a noisy version of ${\mathit{DOTA}}_{0.5}$: ${\mathit{DOTA}}_{0.5}$+noise, by adding a Gaussian noise on each image ($\sigma =0.3$) to test the resilience of the methods against noise. Moreover, we evaluate models on data provided by ADS. These aerial data are sub-sampled to the desired resolution, emulating the above-mentioned satellite characteristics. This dataset is unlabeled; thus, it will only serve to evaluate the models qualitatively and was not used for training.

2.2. Point Process for Object Detection

Point Processes model configurations of points (a finite non-ordered set $\mathbf{y}=\{y^1,\cdots ,y^n\}$ of elements of $\mathcal{S}$) as realizations of a random variable $\Phi$ in the set of all possible configurations $\mathcal{Y}={\bigcup }_{n=0}^{\infty }{(\mathcal{S}\times \mathcal{M})}^n$ (with an arbitrary amount of points). Space $\mathcal{S}$ corresponds to the image space, and $\mathcal{M}$ to the mark space. The marks can be any random variable that describes the point beyond its location, from the radius of a circle to a discrete categorization of the object. Additionally, we denote $n\left(\mathbf{y}\right)$ the number of points in the configuration $\mathbf{y}$ and ${\mathcal{Y}}^n$ the set of all n points configurations.

As we look to detect vehicles, we model points as rectangles described by the following marks: width $y_a$, length $y_b$, and angle $y_{\alpha }$ as shown in Figure 1 (with the mark space defined as $\mathcal{M}=[a_{min},a_{max}]\times [b_{min},b_{max}]\times [0,\pi ]$).

A Point Process $\Phi$ is described by its density h relative to the uniform Point Process [17]. The model of selection and interaction of points derives from an energy U, through a Gibbs density:

$$h\left(\mathbf{y}|\mathbf{X}\right)=\frac{1}{Z}exp\left(-U\left(\mathbf{y},\mathbf{X},\theta \right)\right),$$

(1)

where $\mathbf{X}$ is the image, $\theta$ the parameters of the model, and Z is a normalizing constant. This constant is intractable due to the non-fixed dimension of $\mathcal{Y}$. However, it can be shown that Z exists (The existence and finiteness of Z are sufficient as the following computations will consider ratios of h; thus Z will cancel out.) and thus h is properly defined as the energy per point is bounded (see Appendix A.2). In short, the energy function U measures the compatibility of the configuration $\mathbf{y}$ with the image $\mathbf{X}$; the lower the energy, the higher the compatibility (see Energy-Based Models [26]). It follows that the most compatible output ${\mathbf{y}}^{*}$ minimizes the energy U (i.e., maximizes the density h).

2.2.1. Markovian Point Process

For most applications [18,20,21,23], the physical constraints of the system that is being modeled imply that the marginal density of a point $y\in \mathbf{y}$ should only depend on a neighborhood around it (e.g., vehicles only align with others within a limited distance). This property is formalized through the Markovianity of the Point Process.

Definition 1.

Point Process $\Phi$ is a Markov Point Process under the symmetric and reflexive relation∼if and only if, for every configuration $\mathbf{y}\in \mathcal{Y}$ such that $h\left(\mathbf{y}\right)>0$ [17]:

• $\forall \mathbf{x}\subset \mathbf{y},h\left(\mathbf{x}\right)>0$
• For every point $u\in \mathcal{S}$, $h\left(\right\{u\}\cup \mathbf{y})/h\left(\mathbf{y}\right)$ only depends on u and its neighbors according to ∼ (i.e., $\{y\in \mathbf{y},y\sim u\}$).

Theorem 1.

For a Point Process derived from the energy:

$$U(\mathbf{y},\mathbf{X},\theta )=\sum _{y\in \mathbf{y}}V(y,{\mathcal{N}}_y^{\mathbf{y}},\mathbf{X},\theta ),$$

(2)

where ${\mathcal{N}}_y^{\mathbf{y}}=\{y^{\prime }\in \mathbf{y},y^{\prime }\sim y,y\ne y^{\prime }\}$ is the set of neighbors of y for relation ∼. Then, the Point Process is Markovian for relation ${\sim }^2$ with $y{\sim }^2{y}^{\prime }\equiv \exists u\in \mathcal{S}\times \mathcal{M},y\sim u\sim y^{\prime }$.

The proof is provided in Appendix A.1.

In short, if the energy $V(y,{\mathcal{N}}_y^{\mathbf{y}},\mathbf{X},\theta )$ of a point y depends on its neighbors within a distance of $d_{\mathit{max}}$, the Point Process is Markovian for a distance $2d_{\mathit{max}}$. The Markovianity will prove useful to simplify the sampling procedure and run it in parallel over the whole image.

2.2.2. Papangelou Conditional Intensity

The reference Poisson Point Process has an intensity $\lambda$ that is either constant or depends on the location (For any compact $A\subseteq \mathcal{S}$, $\mathbb{E}\left[n(\Phi )\right]=\lambda \left|A\right|$ if $\lambda$ is constant.). The density in (1) implies the intensity is now a function of the location and neighborhood of a point:

Definition 2.

The Papangelou conditional intensity $\lambda (·;·)$ [17] associated with a Point Process $\Phi$, is defined as:

$$\lambda (y;\mathbf{y})dy=p\left(n_{\Phi }\left(dy\right)=1|\Phi \cap {\left(dy\right)}^c=\mathbf{y}\cap {\left(dy\right)}^c\right),$$

(3)

i.e., the infinitesimal probability of finding a point in region $dy$around$y\in \mathcal{S}$, given the configuration$\mathbf{y}$outside$dy$ (i.e., ${\left(dy\right)}^c$). When $y\in \mathbf{y}$, the Papangelou conditional intensity can be computed from the energy function U as:

$$\begin{array}{cc}\hfill \lambda (y;\mathbf{y}\setminus \{y\left\}\right)& =\frac{h\left(\mathbf{y}\right|\mathbf{X})}{h(\mathbf{y}\setminus \{y\left\}\right|\mathbf{X})}\hfill \\ \hfill & =exp\left(U(\mathbf{y}\setminus \{y\},\mathbf{X},\theta )-U(\mathbf{y},\mathbf{X},\theta )\right).\hfill \end{array}$$

(4)

2.3. Energy-Based Model

The energy formulation of the Point Process allows the compositing of several simple point behaviors into one model. We illustrate the composition of several simple point distributions through the addition of energy terms in Figure 2: There, we show that the point distribution in (Figure 2d) can be obtained by simply adding the energies used to produce Figure 2b,c.

When performing object detection, the image provides information about the location and properties of objects (e.g., location-based energy as in Figure 2b), while prior knowledge of the objects of interest complements this information (e.g., a repulsion term as in Figure 2c); the composability of the energy model allows the factoring in of both pieces of information into a single model (as in the composition of Figure 2d).

More specifically: For every point $y\in \mathbf{y}$ we compute a set of energy terms $V_e(\mathbf{y},\mathbf{X},{\mathcal{N}}_y^{\mathbf{y}})$, $e\in \xi$ and combine those—per point—as a weighted sum [20,21,23]:

$$\begin{array}{cc}\hfill U\left(\mathbf{y},\mathbf{X},\theta \right)& =\sum _{y\in \mathbf{y}}V\left(y,\mathbf{X},{\mathcal{N}}_y^{\mathbf{y}},\theta \right)\hfill \\ \hfill & =\sum _{y\in \mathbf{y}}{w}_0+\sum _{e\in \xi }{w}_e{V}_e\left(y,\mathbf{X},{\mathcal{N}}_y^{\mathbf{y}},\theta \right)\hfill \end{array}$$

(5)

with $V(y,\mathbf{X},{\mathcal{N}}_y^{\mathbf{y}},\theta )$ is the energy contribution of a single object y. Energy terms can be grouped into two categories:

• Prior terms (of the form $V_e(y,{\mathcal{N}}_y^{\mathbf{y}},\theta )$) that only depend on y and its neighborhood ${\mathcal{N}}_y^{\mathbf{y}}$; they measure the coherence of the configuration itself considering the known properties of the objects of interest (e.g., objects do not overlap).
• Data terms ($V_e(y,\mathbf{X},\theta )$) are a function of y and the image $\mathbf{X}$ (also referred to as observation in remote sensing); these terms measure the correspondence of the points with the image.

The weights w are part of the parameters $\theta$ that will need to be set before inference on any unseen image. We discuss the estimation of parameters $\theta$ in Section 2.5.

In the rest of this section, we define the multiple energy terms $V_e$: First, the data terms, then the prior terms. Contrary to previous Point Process (PP) approaches, we do not rely on contrast measures but rather on data terms pulled from the output of a CNN: this allows for more reliable measures and faster computation. Moreover, this new model formulation opens up new improvements in sampling the energy model and estimating its parameters.

2.3.1. Data Terms from a CNN

Classical likelihood measurements [18,20,23] are based on contrast measures that are crafted to fit each application. Moreover, these measures rely heavily on the high contrast between objects and their background, along with limited background diversity. We illustrate the limitations of those in Section 3.1.

On the other hand, CNN models have shown to be very efficient at extracting features from images for object detection and classification. In the following section, we will show how to interpret a CNN-based object detector output to obtain a scalar energy that measures the fitness of a configuration against an image. It allows us to go past the contrast measure design by utilizing a pre-trained CNN output.

Position Likelihood Term

CNN-based object-detection methods such as [7,8,27] make use of a heatmap to find object centers. This object center probability map is obtained by passing the image $\mathbf{X}$ of shape $H,W,3$ (here with 3 color channels) into a fully convolutional model such as Unet [28]. It outputs a tensor ${\widehat{\mathbf{Z}}}_{\mathit{pos}}\in {\mathbb{R}}^{H\times W}$, from which a probability-like measure of an object center at coordinates $(y_i,y_j)$ is obtained as $p(y_i,y_j|\mathbf{X})=\sigma \left({\widehat{\mathbf{Z}}}_{\mathit{pos}}[y_i,y_j]\right)$ where $\sigma (·)$ is the sigmoid function, and ${\widehat{\mathbf{Z}}}_{\mathit{pos}}[y_i,y_j]$ the value of map ${\widehat{\mathbf{Z}}}_{\mathit{pos}}$ at coordinates $(y_i,y_j)$.

We propose to reinterpret this output as the position energy as follows:

$$V_{\mathit{pos}}\left(y,\mathbf{X},\theta \right)=ln\left(1+exp\left(-{\widehat{\mathbf{Z}}}_{\mathit{pos}}[y_i,y_j]+t_{\mathit{pos}}\right)\right),$$

(6)

with $t_{\mathit{pos}}$ a scalar threshold, allowing the movement of the inflection point of the Softplus function $x\mapsto ln(1+exp(-x\left)\right)$.

Mark Likelihood Term

To compute the energy associated with marks, we discretize each mark $\kappa$ ($\kappa \in \{a,b,\alpha \}$ in the case of a Point Process (PP) of rectangles), into $n_{\kappa }$ classes (or value bins) in the range $[{\kappa }_{min},{\kappa }_{max}]$. We define the integer class of a value v for mark $\kappa$ as:

$$c_{\kappa }\left(v\right)=n_{\kappa }\frac{v-{\kappa }_{min}}{{\kappa }_{max}-{\kappa }_{min}},\forall v\in [{\kappa }_{min},{\kappa }_{max}],$$

(7)

with $\lceil c_{\kappa }\left(v\right)\rceil$ being the corresponding integer class.

Now, supposing we have a model trained to classify the mark of an object at a given position $(i,j)\in \mathcal{S}$, such a model will produce a probability estimate of mark $\kappa$ to be in class $c\in 〚1,n_{\kappa }〛$, from model output ${\widehat{\mathbf{Z}}}_{\kappa }^{y_i,y_j}\in {\mathbb{R}}^{n_{\kappa }}$ as:

$$\widehat{p}\left(c\right|i,j,\mathbf{X})=Softmax{\left({\widehat{\mathbf{Z}}}_{\kappa }^{y_i,y_j}\right)}_c=\frac{exp\left({\widehat{\mathbf{Z}}}_{\kappa }^{y_i,y_j}\left[c\right]\right)}{{\sum }_{c^{\prime }=1}^{n_{\kappa }}exp\left({\widehat{\mathbf{Z}}}_{\kappa }^{y_i,y_j}\left[c^{\prime }\right]\right)}.$$

(8)

As in [29], we can reinterpret the model output into a potential, giving:

$$V_{\kappa }(y,\mathbf{X})=-{\widehat{\mathbf{Z}}}_{\kappa }^{y_i,y_j}\left[\lceil c_{\kappa }\left(y_{\kappa }\right)\rceil \right]+ln\left(\sum _{c=1}^{n_{\kappa }}exp\left({\widehat{\mathbf{Z}}}_{\kappa }^{y_i,y_j}\left[c\right]\right)\right).$$

(9)

As the CNN outputs a tensor ${\widehat{\mathbf{Z}}}_{\kappa }$ of shape $(H,W,n_{\kappa })$, we obtain vector ${\widehat{\mathbf{Z}}}_{\kappa }^{y_i,y_j}$ by sampling tensor ${\widehat{\mathbf{Z}}}_{\kappa }$ at location $(\lceil y_i\rceil ,\lceil y_j\rceil )$, as illustrated in Figure 3a. Vector ${\widehat{\mathbf{Z}}}_{\kappa }^{y_i,y_j}$ gives an estimate of the probability of each class c once passed through a Softmax (as given by (8)); We illustrate this with the histogram, on top of ${\widehat{\mathbf{Z}}}_{\kappa }^{y_i,y_j}$ in Figure 3a. The value of $V_{\kappa }(y,\mathbf{X})$ is then obtained by applying the formula in (9) over vector ${\widehat{\mathbf{Z}}}_{\kappa }^{y_i,y_j}$ (as illustrated on the right-hand side of Figure 3).

Interpolation

The energies described in (6) and (9) are defined over integer coordinates in $\mathcal{S}$ and integer classes in $\mathcal{M}$. However, points y in the configuration $\mathbf{y}$ have real-valued locations $(y_i,y_j)\in {\mathbb{R}}^2$ as is for marks in $\mathcal{M}$. Thus, we propose interpolating values between exact pixel locations. For (9) we now have:

$$\begin{array}{cc}\hfill V_{\kappa }(y,\mathbf{X})& ={\mathbf{V}}_{\kappa }[y_i,y_j,c_{\kappa }\left(y_{\kappa }\right)]\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\mathbf{V}}_{\kappa }[i,j]& =-{\widehat{\mathbf{Z}}}_{\kappa }[i,j]+ln\sum _{c=1}^{n_{\kappa }}exp\left({\widehat{\mathbf{Z}}}_{\kappa }[i,j,c]\right),\forall i,j\in {\mathcal{S}}_{\mathit{d}}.\hfill \end{array}$$

(11)

The above has two benefits: first, (10) defines the energy for continuous values of the marks using bilinear interpolation, thus ensuring Lipschitz continuity of the energy. Second, we can pre-compute (11) once for a given image, making the mark energy computation a simple value lookup and interpolation for each point. In practice, we proceed as illustrated in Figure 3b, where a function is applied to the CNN output ${\widehat{\mathbf{Z}}}_{\kappa }$ to obtain ${\mathbf{V}}_{\kappa }$, which is stored in memory. The interpolation gives $V_{\kappa }$ for any continuous position and mark. This pre-computation is quite fast, as it is defined as an operation over a tensor, which is implicitly parallelized by tensor computation libraries such as PyTorch [30]. The same holds for the position term (6).

2.3.2. Energy Priors

The total energy model encompasses several priors as energies, which allows regularizing configuration against the data terms.

2.3.2.1. Object Priors

These are functions of the current point y only. For $k=\mathit{ratio},\mathit{area}$:

$$V_k(y,\theta )=-exp\left(-0.5{\left(f_k\left(y\right)-{\mu }_k\right)}^2{\sigma }_k^{-2}\right)$$

(12)

with $f_{\mathit{ratio}}\left(y\right)=y_b/y_a$ and $f_{\mathit{area}}\left(y\right)=y_a{y}_b$, and ${\mu }_k$,${\sigma }_k$ being, respectively, the target value and the standard deviation. These two are illustrated in the first line of Figure 4.

For our objects of interest, we find two modes in the distribution of areas and ratios: Cars and trucks. To make sure to only favor these two modes and not objects with the area of a truck and ratio of a car, we introduce a joint area and ratio prior:

$$V_{\mathit{jntAR}}(y,\theta )=-exp\left(-\frac{{\left(\frac{y_a}{y_b}-{\mu }_{\mathit{ratio}}\right)}^2}{2{\sigma }_{\mathit{ratio}}^2}-\frac{{\left(y_a{y}_b-{\mu }_{\mathit{area}}\right)}^2}{2{\sigma }_{\mathit{area}}^2}\right).$$

(13)

2.3.2.2. Interaction Priors

The following priors depend on the neighborhood of the point y. The term in (14) penalizes overlapping objects (with $t_{overlap}$ the overlap threshold), (15) favors aligned objects ($t_{\alpha }=0$ favors parallel objects, while $t_{\alpha }=\pi /2$ prefers perpendicular objects), (16) and (17) are, respectively, repulsive and attractive priors. Finally, (18) allows the adjustment of the energy of neighborless points. Some of these priors are illustrated in Figure 4.

$$\begin{array}{cc}\hfill V_{\mathit{ovrlp}}\left(y,{\mathcal{N}}_y^{\mathbf{y}}\right)& =\underset{y^{\prime }\in {\mathcal{N}}_y^{\mathbf{y}}}{max}\left\{Relu\left(\frac{\mathrm{area}\left(y^{\prime }\cap y\right)}{min\{\mathrm{area}\left(y^{\prime }\right),\mathrm{area}\left(y\right)\}}-t_{\mathit{ovrlp}}\right)\right\}\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill V_{\mathit{align}}\left(y,{\mathcal{N}}_y^{\mathbf{y}}\right)& =\underset{y^{\prime }\in {\mathcal{N}}_y^{\mathbf{y}}}{min}\left\{-cos\left(\left|y_{\alpha }-y_{\alpha }^{\prime }\right|\right)\right\}\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill V_{\mathit{repls}}\left(y,{\mathcal{N}}_y^{\mathbf{y}}\right)& =\underset{y^{\prime }\in {\mathcal{N}}_y^{\mathbf{y}}}{max}\left\{Relu\left(1-\frac{d\left(y,y^{\prime }\right)}{d_{max}}-t_{\mathit{repls}}\right)\right\}\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill V_{\mathit{attrc}}\left(y,{\mathcal{N}}_y^{\mathbf{y}}\right)& =\underset{y^{\prime }\in {\mathcal{N}}_y^{\mathbf{y}}}{min}\left\{Relu\left(\frac{d\left(y,y^{\prime }\right)}{d_{max}}-t_{\mathit{attrc}}\right)\right\}\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill V_{\mathit{zrNbr}}\left(y,{\mathcal{N}}_y^{\mathbf{y}}\right)& ={\mathbb{1}}_{|{\mathcal{N}}_y^{\mathbf{y}}|=0}\hfill \end{array}$$

(18)

where $Relu\left(x\right)=max(0,x)$ for $x\in \mathbb{R}$. Note we introduced a few parameters ${\mu }_{\mathit{ratio}},{\sigma }_{\mathit{ratio}},$$t_{\mathit{ovrlp}},\cdots$, these will have to be either set by trial and error or through the parameter estimation method presented in Section 2.5.

2.3.3. Model Pipeline

We summarize the implementation of the overall energy model in Figure 5. As is, our architecture has the following advantages:

• We simplify the data term computation from complex contrast measures [18,20,23] into a simple sampling and interpolation of a tensor (highlighted in green in Figure 5).
• The backbone CNN architecture remains very simple as it does not need to model any object interactions.
• The CNN inference—which is the most complex operation within the graph in Figure 5 (FCN block)—is limited to once per image, as the energy U can be computed for any configuration $\mathbf{y}\in \mathcal{Y}$ without having to infer ${\mathbf{V}}_{\mathit{pos}}$ and ${\mathbf{V}}_{\kappa }$ a second time (pre-computed maps highlighted in lilac in Figure 5).
• New interactions or priors can be easily added to the model as energies are aggregated through a simple linear combination (see right-hand part of the pipeline in Figure 5).
• The implementation with tensor computation libraries (here PyTorch [30]) allows for implicit parallelization on GPU (or CPU) using the tensors batch dimension.
• The latter also enables automatic differentiation of energy U, both relative to configuration $\mathbf{y}$ and to parameters $\theta$, which will be, respectively, useful in Section 2.4.3 and Section 2.5.4.
• Lastly, as the energy maps are defined over a raster space (thus easy to normalize), those will be of use to design the birth map in Section 2.4.1.

2.4. Sampling Configurations from the Energy Model

For a given energy model U (supposing parameters $\theta$ are set), we aim at sampling configurations that follow the Gibbs density defined in (1). This will allow us to find the most compatible configuration ${\mathbf{y}}^{*}$ for any image $\mathbf{X}$, i.e.,:

$${\mathbf{y}}^{*}=\underset{\mathbf{y}\in \mathcal{Y}}{arg\; min}U(\mathbf{y},\mathbf{X},\theta ).$$

(19)

As the space $\mathcal{Y}$ in which we are trying to minimize the energy is not of fixed density, we must resort to an adapted sampling method such as Reversible Jump Monte Carlo Markov Chain (RJMCMC) [31]. While this method ensures proper sampling of the desired law, it often requires some application-specific adaptations to improve sampling time and efficiency (e.g., [20,23]).

Here, we propose to use the already computed energy maps ${\mathbf{V}}_{\mathit{pos}}$ and ${\mathbf{V}}_{\kappa }$ to add new points in relevant areas and to focus on the parallel sampling of the algorithm. Moreover, the implementation of our energy model allows us to leverage automatic gradient computation to implement diffusion dynamics and explore at a fixed number of points (i.e., within ${\mathcal{Y}}_n$) guided by the whole energy model.

2.4.1. Birth Maps for Efficient Point Proposals

At their core, both RJMCMC [31] (an adaptation of the Metropolis-Hastings algorithm [32] to variable dimension problems) and Jump Diffusion [33] make use of Birth and Death moves to build a Markov chain ${\left({\mathbf{y}}_t\right)}_{t=1,\cdots }$ of stationary density $h{(·|\mathbf{X})}^{1/T_t}$ where $T_t$ is a temperature parameter that decreases towards zero (i.e., Simulated Annealing) (Convergence is proven for a logarithmic temperature decrease, here (and in the literature) we approximate it with a geometric decrease at it is faster [17]). That way, the chain converges towards the global minimum energy. Proposed Birth (addition of points to the current $y_t$) and Deaths (removal of points) are accepted or rejected given an acceptance ratio that depends on the energy change induced by the respective addition or removal of points (for more details see [31,34]).

The simplest birth move is to propose new points uniformly in $\mathcal{S}\times \mathcal{M}$. However, this proves highly inefficient as the density of objects in the image varies a lot. Ideally, the birth move would propose a new point u to the current configuration $\mathbf{y}$, sampled with the marginal density $p\left(u\right|\mathbf{y})$ as is carried out within Gibbs sampling (maximizing the chance for the point addition to be accepted). Knowing that $p\left(u\right|\mathbf{y})\propto h(\mathbf{y}\cup \{u\left\}\right)/h\left(\mathbf{y}\right)$, we could compute the marginal density using energy change induced by adding the point. However, this cannot easily be computed because of the interaction terms. Thus, we propose to approximate it by only considering the data terms (i.e., bypassing the energy change on prior terms):

$$\frac{h(\mathbf{y}\cup \{u\left\}\right)}{h\left(\mathbf{y}\right)}\simeq exp\left(-w_{\mathit{pos}}{V}_{\mathit{pos}}(u,\mathbf{X})-\sum _{\kappa \in \{a,b,\alpha \}}{w}_{\kappa }{V}_{\kappa }(u,\mathbf{X})\right).$$

(20)

From this, we define the density d to sample a new point u:

$$d\left(u\right)=\frac{1}{Z_d}exp\left(-w_{\mathit{pos}}{\mathbf{V}}_{\mathit{pos}}\left[u\right]-\sum _{\kappa \in \{a,b,\alpha \}}{w}_{\kappa }{\mathbf{V}}_{\kappa }\left[u\right])\right),$$

(21)

with ${\mathbf{V}}_{\mathit{pos}}$ and ${\mathbf{V}}_{\kappa }$ the pre-computed position and marks energy maps. Since both are defined as tensors with a finite number of elements, the normalizing constant $Z_d$ in (21) can be computed over this discretized space of pixels and mark classes.

This allows for efficient sampling of new points, without any application-specific heuristics, by taking into account a truncated version of the energy model U that can be normalized easily.

2.4.2. Focused Parallel Sampling

In its canonical implementation, the RJMCMC algorithm operates sequentially over the points in the image—It may only add/remove/transform one point at a time —, which increases simulation time linearly with the number of objects. With parallelization, we aim to take advantage of the spatial Markovianity of the Point Process, as Theorem 1 states that moves further than $2d_{\mathit{max}}$ apart are independent.

As in [35], we split the image into square cells of size $d_c$ that are each assigned one set s (each set corresponding to one color in Figure 6). These sets are referred to as mic-sets in [35] for Mutually Independent Cells. Set size $d_c$ is chosen as $2d_{\mathit{max}}+2{\delta }_{\mathit{max}}$ (${\delta }_{\mathit{max}}$ is the maximum point translation distance, introduced in the next Section), so that moves in cells of one set (or color) have independent acceptance ratios: i.e., those moves can be performed in any sequential order –thus in parallel—without any change in the probability of acceptance. In practice, we have $d_c=48$ ($d_{\mathit{max}}=16$, ${\delta }_{\mathit{max}}=8$).

We aim at leveraging the birth map stemming from the pre-computed energy maps ${\mathbf{V}}_{\mathit{pos}}$ and ${\mathbf{V}}_{\kappa }$ to focus the sampling on parts of the image where the density of objects is high. Contrary to [35], we do not use a quadtree structure of cells, but rather a regular grid of $d_c$ sized cells. In order to achieve an efficient sampling of cells (i.e., avoiding spending time in cells with no evidence of objects), weights are assigned to each cell to focus the sampling on the relevant ones. The sampling procedure goes as follows:

• pick a move type (Birth, Death, or Diffusion),
• pick one set s with probability $p\left(s\right)$,
• keep each cell in set s with probability $p\left(c\right|s)$
• for every cell c in set s run move restricted to c (Birth density becomes $d_c\left(u\right)$)

The process of cell selection (step 3) allows for limiting the number of cells processed at once on big images, hence limiting the computational cost and focusing the sampling on areas with high object density. We pick the following densities:

$$\begin{array}{cc}\hfill p\left(s\right)& =d\left(s\right),\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill p\left(c\right|s)& =min(1,n_p\frac{d\left(c\right)}{d\left(s\right)}),\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill d_c\left(u\right)& =\frac{d\left(u\right)}{d\left(c\right)},\hfill \end{array}$$

(24)

denoting (by extension) $d\left(s\right)$ and $d\left(c\right)$ the density d integrated over all cells of set s and cell c, respectively. In practice, parameter $n_p=1$ ensures new points are sampled over the whole image with density d; higher values of $n_p$ allow for more parallel cells to be sampled each step at the cost of straying away from density d.

2.4.3. Jump Diffusion: Leveraging Gradients

The canonical RJMCMC algorithm also uses local transform moves to explore $\mathcal{Y}$ at a fixed number of points, picking one point to translate, rotate, and scale at random in ${\mathbf{y}}_t$. This is wildly inefficient as it does not consider the energy function U it is trying to minimize. Instead, we propose to use the energy gradient to explore the space more efficiently. This is the idea behind stochastic diffusion (or Langevin) dynamics [36,37]. If ${\mathbf{y}}_t$ and $T_t$ denote, respectively, the configuration and temperature at time t, then the configuration for the next step of the Markov chain is given by ${\mathbf{y}}_{t+1}={\mathbf{y}}_t+d{\mathbf{y}}_t$, with step size $\delta$:

$$d{\mathbf{y}}_t=-\delta \frac{\partial U({\mathbf{y}}_t,\mathbf{X},\theta )}{\partial {\mathbf{y}}_t}+d{w}_t\sqrt{2T_t},d{w}_t\sim \mathcal{N}(0,\delta ).$$

(25)

At high temperatures, the Brownian motion from $d{w}_t$ ensures an exploration of space. At low temperatures, the Brownian motion is negligible, and the diffusion performs as a gradient descent. This allows fine-tuning the configuration at low temperatures while considering the whole energy model (contrary to the truncated energy used for birth maps in Section 2.4.1).

The diffusion alone allows for the exploration of space $\mathcal{Y}$ locally around ${\mathbf{y}}_t$ at a fixed dimension (fixed number of points). To explore (or jump) across dimensions, we make use of the Birth and Death moves as defined previously (see Section 2.4.1).

In practice, as parallelization requires some maximum displacement on points (see Section 2.4.2), the step $d{\mathbf{y}}_t$ from ${\mathbf{y}}_t$ to ${\mathbf{y}}_{t+1}$ is clipped; for each point $y\in {\mathbf{y}}_t$ updated to $y^{\prime }$, we bound the i and j components by ${\delta }_{\mathit{max}}$(i.e., $|y_i-y_i^{\prime }|\le {\delta }_{\mathit{max}}$, and similarly for j).

2.4.4. Resulting Sampling Method

The resulting sampling method is outlined in Algorithm 1. The method requires the following inputs:

• ${\mathbf{x}}_0$: initial configuration;
• $\mathbf{X}$: image;
• $\theta$: energy model parameters;
• $T_0$: initial temperature;
• $n_{\mathit{s}}$: number of samples;
• $\alpha$: temperature decay rate.

Here, we denote by $C_c\left(\mathbf{x}\right)$ the restriction of configuration $\mathbf{x}$ to cells c. Please note that to compute U for iteration of the loop on t, and the energy maps ${\mathbf{V}}_{\mathit{pos}}$, ${\mathbf{V}}_{\kappa }$ are only computed once as stated in Section 2.3.3.

Algorithm 1: Sampling method.

Input  ${\mathbf{x}}_0,\mathbf{X},\theta ,T_0,n_{\mathit{s}},\alpha$ 1: for$t=0,\cdots ,n_{\mathit{s}}-1$do 2:    Pick diffusion with probability $0.8$, else jump 3:    Pick mic-set s with probability $p\left(s\right)$ 4:    Keep each c in s with probability $p\left(c\right|s)$ to make $\tilde{s}$ 5:    ${\mathbf{x}}_{t+1}\leftarrow {\mathbf{x}}_t$             unvisited cells will keep the previous state 6:    for all $c\in \tilde{s}$ 7:      if diffusion then 8:        $dw\sim \mathcal{N}(0,\beta )$ 9:        $\Delta C_c\left({\mathbf{x}}_t\right)=-\beta {\nabla }_{C_c\left({\mathbf{x}}_t\right)}U({\mathbf{x}}_t,\mathbf{X},\theta )+dw\sqrt{2T_t}$ 10:        $C_c\left({\mathbf{x}}_{t+1}\right)\leftarrow C_c\left({\mathbf{x}}_t\right)+\Delta C_c\left({\mathbf{x}}_t\right)$ 11:      else 12:        $Q_c\leftarrow Q_{c,\mathrm{B}}$ with probability $0.5$ else $Q_{c,\mathrm{D}}$        pick birth or death 13:        ${\mathbf{x}}^{\prime }\sim Q_c(\mathbf{x}\to ·)$  remove a point at random or add one using the birth map 14:        $r\leftarrow \frac{Q_c({\mathbf{x}}^{\prime }\to \mathbf{x})}{Q_c(\mathbf{x}\to {\mathbf{x}}^{\prime })}exp\left(-\frac{U({\mathbf{x}}^{\prime },\mathbf{X},\theta )-U(\mathbf{x},\mathbf{X},\theta )}{T_t}\right)$        compute Green ratio 15:        $C_c\left({\mathbf{x}}_{t+1}\right)\leftarrow C_c\left({\mathbf{x}}^{\prime }\right)$ with probability $min(1,r)$ 16:      end if 17:   end for 18:   $T_{t+1}\leftarrow \alpha T_t$        decrease temperature 19: end for Output:  ${\mathbf{x}}_{n_{\mathit{s}}}$

The above implements the sampling improvements described previously, namely:

• sampling new points using the birth map (line 13),
• sampling in parallel over the whole image (line 6), in practice using tensor batch dimensions,
• using the energy gradient to perform Diffusion (line 9).

2.5. Estimating Model Parameters

In this section, we estimate the parameters $\theta$, so that, for each image $\mathbf{X}$, the associated ground truth configuration ${\mathbf{y}}^{\mathit{GT}}$ matches the most compatible configuration, i.e.:

$${\mathbf{y}}^{\mathit{GT}}=\underset{\mathbf{y}\in \mathcal{Y}}{arg\; min}U(\mathbf{y},\mathbf{X},\theta )$$

(26)

Previous PP approaches for object detection often relied on trial-and-error parameter estimation or used linear programming [23,38]. The latter method would generate a set of constraints (These constraints can be seen as a local reformulation of global constraint (26).) $U\left({\mathbf{y}}^{-}\right)\le U\left({\mathbf{y}}^{\mathit{GT}}\right)$ where wrong configurations ${\mathbf{y}}^{-}$ are generated by applying strong perturbations on the ground truth configuration ${\mathbf{y}}^{\mathit{GT}}$. However, this only estimates the weights $w_e$ and is prone to over-constraining when the ground truth is noisy or when considering too many constraints. Moreover, this method requires designing the procedure to generate the configurations ${\mathbf{y}}^{-}$, which we find has a great effect on the estimated parameters.

2.5.1. Maximum Likelihood Learning

To tackle the above limitations, we turn towards the Energy-Based Model (EBM) literature for new parameter estimation methods. The authors in [26] propose to learn EBM parameters by maximizing the likelihood for the data $\mathcal{D}$:

$$p\left({\mathbf{y}}_1^{\mathit{GT}},\cdots ,{\mathbf{y}}_{\left|\mathcal{D}\right|}^{\mathit{GT}}|{\mathbf{X}}_1,\cdots ,{\mathbf{X}}_{\left|\mathcal{D}\right|},\theta \right)=\prod _{k=1}^{\left|\mathcal{D}\right|}p\left({\mathbf{y}}_k^{\mathit{GT}}|{\mathbf{X}}_k,\theta \right).$$

(27)

The negative log-likelihood is then given, for parameters ${\theta }_n$ at step n, as:

$${\mathcal{L}}_{\mathit{NLL}}\left({\theta }_n,\mathcal{D}\right)=\frac{1}{\left|\mathcal{D}\right|}\sum _{k=1}^{\left|\mathcal{D}\right|}\left(U\left({\mathbf{y}}_k^{\mathit{GT}},{\mathbf{X}}_k^{+},{\theta }_n\right)+\frac{1}{\beta }log{\int }_{\mathbf{y}\in \mathcal{Y}}exp\left(-\beta U\left(\mathbf{y},{\mathbf{X}}_k,{\theta }_n\right)\right)\right)$$

(28)

with $\beta$ an inverse temperature parameter (from the Gibbs distribution) that does not affect the position of the minimum [26]. The gradient of ${\mathcal{L}}_{NLL}$ on ${\theta }_n$ is shown to be

$$\frac{\partial {\mathcal{L}}_{\mathit{NLL}}\left({\theta }_n,{\mathbf{y}}_i^{\mathit{GT}},{\mathbf{X}}_i\right)}{\partial {\theta }_n}=\frac{\partial U\left({\mathbf{y}}_i^{\mathit{GT}},{\mathbf{X}}_i,{\theta }_n\right)}{\partial {\theta }_n}-{\int }_{\mathbf{y}\in \mathcal{Y}}\frac{\partial U\left(\mathbf{y},{\mathbf{X}}_i,{\theta }_n\right)}{\partial {\theta }_n}p\left(\mathbf{y}|{\mathbf{X}}_i,{\theta }_n\right)$$

(29)

where $p\left(\mathbf{y}\right|{\mathbf{X}}_i,{\theta }_n)=exp(-\beta U(\mathbf{y},{\mathbf{X}}_i,{\theta }_n))/{\int }_{\tilde{\mathbf{y}}\in \mathcal{Y}}exp(-\beta U(\tilde{\mathbf{y}},{\mathbf{X}}_i,{\theta }_n))$.

While the integral ${\int }_{\mathbf{y}\in \mathcal{Y}}\frac{\partial U(\mathbf{y},{\mathbf{X}}_i,{\theta }_n)}{\partial {\theta }_n}p\left(\mathbf{y}\right|{\mathbf{X}}_i,{\theta }_n)$ remains intractable, it can be approximated through Monte Carlo sampling, where ${\tilde{\mathbf{y}}}^0,\cdots ,{\tilde{\mathbf{y}}}^N$ are drawn from the law defined by $p(·|{\mathbf{X}}_i,{\theta }_n)$, which yields:

$${\int }_{\mathbf{y}\in \mathcal{Y}}\frac{\partial U(\mathbf{y},{\mathbf{X}}_i,{\theta }_n)}{\partial {\theta }_n}p\left(\mathbf{y}|{\mathbf{X}}_i,{\theta }_n\right)\simeq \frac{1}{N}\sum _{k=1}^N\frac{\partial U\left({\tilde{\mathbf{y}}}^k,{\mathbf{X}}_i,{\theta }_n\right)}{\partial {\theta }_n}$$

(30)

2.5.2. Contrastive Divergence

The authors in [39,40] propose to use a single sample in their Contrastive Divergence (CD) method. This method also uses a few simulation steps for the Monte Carlo Markov chain (MCMC) to generate ${\mathbf{y}}^{-}$, starting from the desired answer ${\mathbf{y}}^{\mathit{GT}}$.

The general idea is to generate contrastive samples ${\mathbf{y}}^{-}$ that follow the density derived from $U(·,\mathbf{X},{\theta }_n)$ at step n of the optimization. Then we proceed to update ${\theta }_n$ to ${\theta }_{n+1}$, by gradient descent, to minimize the energy of the valid sample ${\mathbf{y}}^{\mathit{GT}}$, while maximizing the energy of the contrastive sample${\mathbf{y}}^{-}$ (see Figure 7). Alternatively, we can augment the data and use positive samples ${\mathbf{y}}^{+}={\mathbf{y}}^{\mathit{GT}}+\mathcal{N}(0,{\sigma }^{+})$ to replace ${\mathbf{y}}^{\mathit{GT}}$ as in [41]. In the case of imperfect GT, this models the uncertainty over the labels while providing some data augmentation.

The loss to minimize is then:

$$\begin{array}{cc}\hfill \mathcal{L}\left({\theta }_n,{\mathbf{y}}^{+},{\mathbf{y}}^{-},\mathbf{X}\right)& =U({\mathbf{y}}^{+},\mathbf{X},{\theta }_n)-U({\mathbf{y}}^{-},\mathbf{X},{\theta }_n)+\gamma R_V\hfill \\ \hfill R_V& =\sum _{\mathbf{y}\in \{{\mathbf{y}}^{+},{\mathbf{y}}^{-}\}}\frac{1}{n\left(\mathbf{y}\right)}\sum _{y\in \mathbf{y}}\left|V\left(\mathbf{y},\mathbf{X},{\mathcal{N}}_y^{\mathbf{y}},{\theta }_n\right)\right|,\hfill \end{array}$$

(31)

with $\gamma >0$ the weight of regularization term $R_V$. We introduce the regularization term to avoid an explosion of the per-point energy $V(y,\mathbf{X},{\mathcal{N}}_y^{\mathbf{y}},{\theta }_n)$. To ensure a sparse weighting of energies (i.e., minimize the number of non-zero weights $w_e$) we can introduce a new regularization term as an $L^1$ norm on the vectors of weights [42]:

$$R_1=\sum _{e\in \xi }\left|w_e\right|.$$

(32)

The broad strokes of the estimation procedure for parameters $\theta$ are as follows:

• Pick a pair $({\mathbf{y}}^{\mathit{GT}},\mathbf{X})$ from data $\mathcal{D}$.
• Generate positive sample ${\mathbf{y}}^{+}={\mathbf{y}}^{\mathit{GT}}+\mathcal{N}(0,{\sigma }^{+})$.
• Generate negative sample ${\mathbf{y}}^{-}\sim exp(-{\beta }^{-1}U({\mathbf{y}}^{-},\mathbf{X},{\theta }_n))$.
• Compute the loss $\mathcal{L}\left({\theta }_n,{\mathbf{y}}^{+},{\mathbf{y}}^{-},\mathbf{X}\right)$ (see (31)).
• Update ${\theta }_n$ to ${\theta }_{n+1}$ according to the gradient $\nabla \mathcal{L}$, with the Stochastic Gradient Descent (SGD) method [43], as illustrated in Figure 7.
• Loop back to step 1 until convergence.

2.5.3. Replay Buffer

As sampling contrastive samples ${\mathbf{y}}^{-}$ can be time-consuming, we use the replay buffer introduced by [41]. The replay buffer saves Markov chain results at the current optimization step to use for initialization in the next steps, thus saving computing time. This allows for reducing the long simulation time necessary to pass the burn-in period of the Markov Chain [44] and obtain samples ${\mathbf{y}}^{-}$. We use a sample from the law derived from $U(·,\mathbf{X},{\theta }_{n-1})$ to initialize the chain that samples the law derived from $U(·,\mathbf{X},{\theta }_n)$. The use of the replay buffer within one step n of the optimization loop goes as follows:

• With probability $p_{\mathcal{B}}$ set configuration ${\mathbf{x}}_0$ as a value of the replay buffer $\mathcal{B}$ picked at random. With probability $1-p_{\mathcal{B}}$ (or if the buffer is empty) set configuration ${\mathbf{x}}_0$ as a random configuration in $\mathcal{Y}$.
• Run the Markov Chain ${\left({\mathbf{x}}_t\right)}_{t=0,\cdots ,n_{\mathit{s}}}$ to simulate model of energy $U(·,\mathbf{X},{\theta }_n)$.
• Use ${\mathbf{y}}^{-}={\mathbf{x}}_{n_{\mathit{s}}}$ as a negative sample for the loss computation and update of ${\theta }_n$ to ${\theta }_{n+1}$.
• Update buffer $\mathcal{B}\leftarrow \mathcal{B}\cup \left\{{\mathbf{y}}^{-}\right\}$.

This buffer allows running longer chains across epochs (i.e., obtaining higher-quality samples) while maintaining a low computational cost.

2.5.4. Parameters Estimation Algorithm

The final parameter estimation algorithm is as follows:

At line 4 we select a temperature at which to run the Markov chain to sample the negative configuration ${\mathbf{y}}^{-}$ (${\beta }_1=1$, ${\beta }_2={10}^{-2}$, ${\beta }_3={10}^{-4}$); higher temperatures allow sample diversity, while lower ones seek points of minimal energy within the current energy model. Scaling $\varsigma$ computed at line 6 allows rescaling the energy to a unit variance model in order to ensure consistency in sample diversity for a set temperature. It is computed over 1000 samples $\mathbf{y}\sim \mathcal{U}\left(\mathcal{Y}\right)$. Finally, for each temperature $\beta$ and image $\mathbf{X}$, we set a specific buffer ${\mathcal{B}}_{\mathbf{X},\beta }$. The parameter K is the number of steps the Markov chain runs for every gradient descent step. In practice, we set $K=0.02d_c^2$ so that the number of iterations scales with the number of pixels in a cell.

Algorithm 2 we propose above provides a single procedure to estimate all parameters in $\theta$ rather than only the weights $w_e$ as done previously with linear programming in [23,38]. On top of avoiding over-constraining by having no hard constraints, our method also removes the need to design a procedure to sample configurations ${\mathbf{y}}^{-}$ as those are sampled from the current energy model $U(·,\mathbf{X},{\theta }_n)$. Finally, the replay buffer allows for qualitative samples while maintaining a limited computational burden.

Algorithm 2: Contrastive Divergence parameters estimation.

1: $\mathcal{B}\leftarrow \varnothing ,n\leftarrow 0$ 2: while not converged do 3:     for all $({\mathbf{y}}^{\mathit{GT}},\mathbf{X})$ in $\mathcal{D}$ 4:      $\beta \sim \mathcal{U}\left(\{{\beta }_1,{\beta }_2,{\beta }_3\}\right)$             select temperature to sample negative config. 5:      ${\mathbf{x}}_0\sim \mathcal{U}\left({\mathcal{B}}_{\mathbf{X},\beta }\right)$ with probability $p_{\mathcal{B}}$, else $\mathcal{U}\left(\mathcal{Y}\right)$           retrieve relevant buffer 6:      $\varsigma \leftarrow \sqrt{Var\left(U_{T=\infty }(·,\mathbf{X},{\theta }_n)\right)}$                  compute energy scale 7:      ${\mathbf{y}}^{-}\leftarrow Sample\left({\mathbf{x}}_0,\mathbf{X},{\theta }_n,\varsigma \beta ,K,1\right)$                   see Algorithm 1 8:      ${\mathbf{y}}^{+}\leftarrow {\mathbf{y}}^{\mathit{GT}}+\mathcal{N}(0,{\sigma }^{+})$ 9:      $\Delta {\theta }_n\leftarrow {\nabla }_{{\theta }_n}\mathcal{L}\left({\theta }_n,{\mathbf{y}}^{+},{\mathbf{y}}^{-},\mathbf{X}\right)$ 10:      Update ${\theta }_{n+1}$ with $\Delta {\theta }_n$ using SGD 11:      ${\mathcal{B}}_{\mathbf{X},\beta }\leftarrow \left\{{\mathbf{y}}^{-}\right\}\cup {\mathcal{B}}_{\mathbf{X},\beta }$                       update the buffer 12:      $n\leftarrow n+1$ 13:   end for 14: end while

2.6. Papangelou Intensity as a Confidence Score

Classical CNN-based object-detection models for object detection (such as [8,45,46]) yield a confidence score $s\left(y\right)\in \mathbb{R}$ for each proposed object y in the image. This confidence score is often interpreted, for each detection, as proportional to the probability of proposed element y to be a true positive, $s\left(y\right)\propto p\left(y\right|\mathbf{X})$. Applying a score (or confidence) threshold $t_s$ gives a set of detections for which metrics such as precision and recall can be computed by matching the detections with the ground truth. This allows adapting the threshold according to the need for the application; i.e., some applications may require low false positives (high precision) while others require fewer missed detections (high recall). To assess the performance independently of the threshold selection, the Average Precision (AP) metric sums up the performance of a model as the area under the Precision-Recall curve.

Previous Point Process approaches [18,47,48] only compute simple metrics such as precision, recall, or F1 score for the configuration given by the sampling procedure, as no score is associated with each object detection.

2.6.1. Papangelou Intensity as Score

With our PP approach, we propose to introduce a scoring function, first to filter the detections given a confidence threshold and second to be able to compare our method to others using the widely used AP metric. Within the PP framework, the probability of one proposed point being an object of interest depends on the rest of the inferred configuration $\widehat{\mathbf{y}}$; thus the scoring function reflects it: $s\left(y\right|\widehat{\mathbf{y}}\setminus \left\{y\right\})\propto p\left(y\right|\widehat{\mathbf{y}}\setminus \left\{y\right\},\mathbf{X})$. From (3), we have that the Papangelou conditional intensity is proportional to the probability of finding a point $y\in \mathbf{y}$ in a small neighborhood $dy$ knowing the rest of the configuration $\mathbf{y}\setminus \left\{y\right\}$. We propose to use the Papangelou conditional intensity as a score:

$$s\left(y\right|\mathbf{y}\setminus \left\{y\right\})=\lambda (y;\mathbf{y}\setminus \left\{y\right\})$$

(33)

2.6.2. Pruning Sequence

However, the dependency of the score on the current configuration yields a complication while computing the AP: when applying a threshold $t_s$ to prune the configuration $\mathbf{y}$ into ${\mathbf{y}}^{\prime }\subset \mathbf{y}$, for any $y\in {\mathbf{y}}^{\prime }$, the score $s\left(y\right|\mathbf{y}\setminus \left\{y\right\})$ may differ from $s\left(y\right|{\mathbf{y}}^{\prime }\setminus \left\{y\right\})$. With a score of the form $s\left(y\right)$ that only depends on y and the image $\mathbf{X}$—such as those from classical CNN models—the score from one object after pruning is unchanged.

In the PP case, we compute the scores by sequentially removing the lowest scoring point until none is left; i.e., we build a sequence of configurations ${\mathbf{y}}_1\supset {\mathbf{y}}_2\cdots {\mathbf{y}}_{n\left(\widehat{\mathbf{y}}\right)-1}\supset {\mathbf{y}}_{n\left(\widehat{\mathbf{y}}\right)}\supset \varnothing$, with ${\mathbf{y}}_1=\widehat{\mathbf{y}}$, for $n=1,\cdots ,|\widehat{\mathbf{y}}|$:

$$\begin{array}{c}{\mathbf{y}}_{n+1}={\mathbf{y}}_n\setminus \left\{y_n\right\},y_n=\underset{y\in {\mathbf{y}}_n}{\arg \min }\lambda (y;{\mathbf{y}}_n\setminus \left\{y\right\}),\end{array}$$

(34)

$$\begin{array}{c}s\left(y_n\right|{\mathbf{y}}_n\setminus \left\{y_n\right\})=s\left(y_n\right|{\mathbf{y}}_{n+1})=\lambda (y;{\mathbf{y}}_{n+1}).\end{array}$$

(35)

Equation (34) provides a pruning order $y_1,\cdots ,y_{n\left(\widehat{\mathbf{y}}\right)}$ of points in $\widehat{\mathbf{y}}$. This ordering allows plotting the precision and recall curve. Indeed, to trace a Precision-Recall curve, one only requires the sequence of $(\mathrm{Recall}\left(t_s\right),\mathrm{Precision}\left(t_s\right))$ pairs, which are obtained by sequentially pruning the lowest scoring points. Equation (35) provides a score for each point $y_n$.

2.6.3. Contrastive Divergence Loss and Papangelou Intensity

On the one hand, the energy model is trained by minimizing the loss function in (31) derived from the likelihood maximization of the parameters regarding the annotated data. On the other, we evaluate the performance of the inferred configuration with the scoring method in (35) sourced from the Papangelou intensity. Here, we show that while the two are derived differently, the minimization of the loss function leads to good properties on the score function.

Here, we consider a simplified loss with only the two energy terms (as $\gamma \simeq 0$). Denoting the energy change induced by the move from configuration $\mathbf{y}$ to $\mathbf{x}$ as $\Delta U(\mathbf{y}\to \mathbf{x})=U\left(\mathbf{x}\right)-U\left(\mathbf{y}\right)$, we have:

$$\mathcal{L}(\theta ,{\mathbf{y}}^{+},{\mathbf{y}}^{-})=\Delta U({\mathbf{y}}^{-}\to {\mathbf{y}}^{+}).$$

(36)

Similarly, the Papangelou intensity can be rewritten as such:

$$\lambda (u;\mathbf{y})=exp(\Delta U(\mathbf{y}\cup \left\{u\right\}\to \mathbf{y}\left)\right)$$

(37)

Single Point Addition

Thus, for a simple negative sample ${\mathbf{y}}^{-}={\mathbf{y}}^{+}\cup \left\{u\right\}$ in which we add a non-valid point u to ${\mathbf{y}}^{+}$, we have:

$$\mathcal{L}(\theta ,{\mathbf{y}}^{+},{\mathbf{y}}^{-})=log\left(\lambda (u;{\mathbf{y}}^{+})\right).$$

(38)

This leads to the expected behavior: minimizing the loss $\mathcal{L}$ leads to minimizing the score of non-valid point u. The same stands for the removal of a valid point $y\in {\mathbf{y}}^{+}$ and maximizing its score.

Arbitrary Sequence of Moves

This is also valid for the generic case where ${\mathbf{y}}^{-}$ is generated from an arbitrary sequence of additions or removal of points (A translation/rotation/scaling can be viewed as removal and addition.) from ${\mathbf{y}}^{+}$. This defines a sequence ${\left({\mathbf{y}}_k\right)}_{k=0,\cdots ,n}$ of n configurations as:

$$\forall k=1,\cdots ,n,{\mathbf{y}}_k=\left\{\begin{array}{cc}{\mathbf{y}}_{k-1}\setminus \left\{y_k\right\}\hfill & \mathrm{if}{y}_k\in {\mathbf{y}}^{+}\hfill \\ {\mathbf{y}}_{k-1}\cup \left\{y_k\right\}\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(39)

with ${\mathbf{y}}_0={\mathbf{y}}^{+}$, ${\mathbf{y}}^{-}={\mathbf{y}}_n$, and $y_k$ elements of either $\mathcal{S}\times \mathcal{M}$ or ${\mathbf{y}}^{+}$. Without loss of generality, we can reorder the sequence to match the pruning order defined in (34). The energy change for one move is given as:

$$\Delta U({\mathbf{y}}_{k-1}\to {\mathbf{y}}_k)=\left\{\begin{array}{cc}log\left(\lambda (y_k;{\mathbf{y}}_{k-1}\setminus \left\{y_k\right\})\right)\hfill & \mathrm{if}{y}_k\in {\mathbf{y}}^{+}\hfill \\ -log\left(\lambda (y_k;{\mathbf{y}}_{k-1})\right)\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(40)

As we have (by definition) $\Delta U(\mathbf{x}\to {\mathbf{x}}^{\prime \prime })=\Delta U(\mathbf{x}\to {\mathbf{x}}^{\prime })+\Delta U({\mathbf{x}}^{\prime }\to {\mathbf{x}}^{\prime \prime })$, the loss is given as:

$$\mathcal{L}=\sum _{y_k\notin {\mathbf{y}}^{+}}log(\underset{(\mathrm{a})}{\underbrace{\lambda (y_k;{\mathbf{y}}_{k-1})}})-\sum _{y_k\in {\mathbf{y}}^{+}}log(\underset{(\mathrm{b})}{\underbrace{\lambda (y_k;{\mathbf{y}}_{k-1}\setminus \{y_k\})}}).$$

(41)

By ordering the $y_k,{\mathbf{y}}_k$ to match the pruning order in (34) each $\lambda (y_k;\cdots )$ can be matched to their respective score:

• (41)(a) corresponds to non-valid points added to ${\mathbf{y}}^{+}$, their score is minimized as the loss is decreased;
• (41)(b) corresponds to valid points removed from ${\mathbf{y}}^{+}$. Their score is increased as the loss is minimized.

With this, we showed that minimization of the loss at a configuration level leads to the expected results on object scores.

While this score function allows taking into account the interaction between objects that the PP model introduces, the decomposable nature of the PP (see (5)) allows for the interpretability of the results as demonstrated in Section 3.4.

2.7. Method Summary

In this section, we proposed a model for object detection that combines CNN architectures with a Point Process (PP) framework. Before applying our approach to real data, we shall summarize its function here. Considering configurations of objects $cy$, image $\mathbf{X}$, and parameters $\theta$, we look to build an energy model $U(\mathbf{y},\mathbf{X},\theta )$ such that the minimum energy point of it corresponds to the ground truth configuration for that image ${\mathbf{y}}^{\mathit{GT}}$. As such, we build U to be a combination of multiple energy terms (Equation (5)), with prior terms that encode the interaction model of our objects (where high energies repel unwanted configurations and lower ones favor others), and data terms that we build from the features from a CNN (Section 2.3.1). We refer the reader to Figure 5 for an overall view of the model pipeline. Within this framework, we can combine the information from the image with the prior knowledge of the configurations.

To estimate the parameters $\theta$, we propose to use Contrastive Divergence. In short, we initialize the parameters at random and alternate between sampling low-energy configurations with the current parameters—that are nowhere close to the ground truth at the beginning of the procedure—and updating the parameters to increase the energy of these samples while maintaining low energy for the ground truth (Section 2.5). This allows us to estimate any differentiable parameters in the model, while previous methods would often become over-constrained and could only estimate linear parameters.

Now, given parameters $\theta$ are set, we want to estimate the configuration of an unseen image $\mathbf{X}$. This is achieved by sampling the configuration ${\mathbf{y}}^{*}$ that minimizes the energy through a Jump Diffusion mechanism. Alternatively, proposing to add or remove points in the configuration guided by the energy maps the CNN provides (Section 2.4.1) and lowering the energy locally at a fixed number of points with gradient descent (Section 2.4.3). The latter can be applied in parallel over the image while focusing on patches of the image with higher expected object density (Section 2.4.2).

Lastly, we propose a metric based on the Papangelou intensity that allows us to compare our approach to others and takes into account the object interactions in the likelihood of detection (Section 2.6).

3. Results

3.1. Comparison to PP Based on Contrast Measures

To demonstrate the need for learned likelihood measures for our data, we evaluate a few contrast measures found in the literature [20,24] in the task of classification of objects and non-objects (Figure 8).

Table 1. Various contrast measures from literature and Average Precision (AP) computed on the DOTA image (${\mathrm{AP}}_{\mathrm{real}}$) or synthetic image (${\mathrm{AP}}_{\mathrm{synthetic}}$) from Figure 8.

Measure	${\mathit{V}}_{\mathit{c}}$ Formulation	${\mathbf{AP}}_{\mathbf{real}}$	${\mathbf{AP}}_{\mathbf{synthetic}}$
CNN output (ours)	$V_{\mathit{pos}}(y,\mathbf{X}){\sum }_{\kappa }{V}_{\kappa }(y,\mathbf{X})$	0.99	0.88
t-test [20]	$\frac{|{\mu }_y-{\mu }_s|}{\sqrt{\frac{{\sigma }_y^2}{n_y}+\frac{{\sigma }_s^2}{n_s}}}$	0.13	0.99
Gradient [24]	$\int n_y\left(t\right)·\frac{\nabla \mathbf{X}\left(s\right(t\left)\right)}{\sqrt{{|\nabla \mathbf{X}\left(s\left(t\right)\right)|}^2+ϵ}}dt$	0.27	0.99
Constant	$1.0$	0.01	0.01

formulates these measures by denoting ${\mu }_y$, ${\sigma }_y^2$ the empirical mean and variance of pixels inside shape y; ${\mu }_s$, ${\sigma }_s^2$ mean and variance on shape contour s; $n_y$ the normal vector of contour s, and $\nabla \mathbf{X}$ the image gradient.

To assert the viability of such measures for these data, we compute a score $s\left(y\right)=exp(-V_c\left(y\right)$ for proposed detections y (both ground truth and randomly scattered rectangles). The true objects represent $1\%$ of this experiment’s dataset (A higher proportion of true objects would allow high average precision for useless estimators predicting a constant value). Here we consider around 32k object proposals (i.e., 313 true objects, and $31,300$ non-objects). A well-suited measure $s\left(y\right)$ should allow classifying easily between ground truth objects and random, non-valid objects. We plot Precision-Recall curves for several of those contrast measures in Figure 8 both for a sample image of ${\mathit{DOTA}}_{0.5}$ and a synthetic highly contrasted image. The Average Prevision (AP) values are reported in Table 1. We compare with the proposed detection energy given by $V_{\mathit{pos}}+V_a+V_b+V_{\alpha }$.

The proposed measure, based on CNN output, outperforms classical contrast-based measures, as the CNN has learned to differentiate between relevant and irrelevant contrasted objects (Irrelevant objects include, for instance, AC units—a highly contrasted, car-sized object—or some pavement features.).

3.2. Models

In this article, we show results on two CNN-based models and two PP models. A Python implementation of our models is available at https://github.com/Ayana-Inria/PP-EBM (accessed on 31 January 2024).

• CNN-LocalMax.: Naive detection from the CNN backbone (used in the PP models); we find objects through local maxima in the output probability maps (or local minima in the energy maps).
• ${\mathrm{CNN-PP}}_m$: PP model with minimal inferred parameters: Pre-set priors (i.e., with manually set parameters), only the weight vector $w_e$ are estimated with the CD method. The values of manually set parameters are detailed in

  Table 2. Values for manually set parameters. Some energy terms such as $V_a$ have no parameters.

  Energy Term	Parameter	Value
  $V_{\mathit{pos}}$	$t_{\mathit{pos}}$	0
  $V_a$ *	-	-
  $V_b$ *	-	-
  $V_{\alpha }$ *	-	-
  $V_{\mathit{align}}$	-	-
  $V_{\mathit{ovrlp}}$	$t_{\mathit{ovrlp}}$	0
  $V_{\mathit{repls}}$	$t_{\mathit{repls}}$	0
  $V_{\mathit{attrc}}$	$t_{\mathit{attrc}}$	0
  $V_{\mathit{jntAR},\mathit{car}}$	${\mu }_{\mathit{ratio}}$	0.46
  	${\mu }_{\mathit{area}}$	42
  	${\sigma }_{\mathit{ratio}}$	0.1
  	${\sigma }_{\mathit{area}}$	20
  $V_{\mathit{jntAR},\mathit{truck}}$	${\mu }_{\mathit{ratio}}$	0.23
  	${\mu }_{\mathit{area}}$	123
  	${\sigma }_{\mathit{ratio}}$	0.1
  	${\sigma }_{\mathit{area}}$	20

  * Mark energy terms $V_{\kappa }$, $\kappa \in \{a,b,\alpha \}$.

  .
• ${\mathrm{CNN-PP}}_w$: PP model with the whole parameter vector $\widehat{\theta }$ (which encompasses the weights $w_e$ and priors’ parameters) estimated with the CD method.
• BBA-Vec. and YOLOV5-OBB: Lastly, we compare all our models above with BBA-Vec. from [8] and YOLOV5-OBB from [46]. These models are trained on the same data as the above-mentioned models.

3.3. Results on Remote-Sensing Data

The models are evaluated on a test split of the ${\mathit{DOTA}}_{0.5}$ dataset (for quantitative and qualitative results) and on ADS (for qualitative results as these data are not annotated). Metrics are computed by matching the object proposals and ground truth using the IOU (Intersection Over Union) with a $0.25$ threshold. The Average Precision (AP) metrics (which are threshold independent) are shown in

Table 3. Average Precision (AP) values. ${\mathrm{AP}}_0$ and ${\mathrm{AP}}_N$ correspond, respectively, to the AP on the ${\mathit{DOTA}}_{0.5}$ and ${\mathit{DOTA}}_{0.5}$+noise datasets. Precision-recall curves from which the AP is derived are shown in Figure 9.

Method	${\mathbf{AP}}_0$	${\mathbf{AP}}_{\mathit{N}}$
BBA-Vec.	0.82	0.19
YOLOV5-OBB	0.86	0.10
CNN-LocalMax.	0.86	0.55
${\mathbf{CNN-PP}}_m$	0.91	0.58
${\mathbf{CNN-PP}}_w$	0.92	0.62

and are derived from Precision-Recall (PR) curves shown in Figure 9. We show detection results in Figure 10 (lines 1–2 for ${\mathit{DOTA}}_{0.5}$ and line 3 for ${\mathit{DOTA}}_{0.5}$+noise), with a score threshold set for each model to the one maximizing the F1 score. The ADS data inference results in Figure 11 are obtained with the models trained on the ${\mathit{DOTA}}_{0.5}$ training set.

Computational Complexity

Our two PP-based methods—running on an Nvidia Quadro RTX 8000 Graphics Processing Unit—execute in 300 s on average, for 16 k parallel iterations (equivalent to 77 k sequential iterations) for one image. With an efficient implementation, considering a density of $1.9\times {10}^{-3}$ ($95\mathrm{th}$ percentile of the observed object densities (i.e., $95\%$ of observed object densities are below this value)), we estimate the cost of one iteration to be around 5 operations per pixel per iteration, thus around $4\times {10}^5$ operations per pixel in total. We show the derivation of those values in Appendix B. As a point of comparison, Transformer models for image classification range from $5\times {10}^6$ to $1.8\times {10}^7$ operations (see [49]). The complexity could be greatly reduced by reducing the number of iterations of the Point Process sampler at the cost of straying away from the proper convergence properties.

3.4. Results Interpretability

Due to the decomposition of the total energy into energy terms introduced in (5), the object score can be decomposed similarly:

$$\begin{array}{c}\hfill s\left(y\right|\mathbf{y}\setminus \left\{y\right\})=\prod _{e\in \xi }{s}_e\left(y\right|\mathbf{y}\setminus \left\{y\right\})\end{array}$$

(42)

with $s_e\left(y\right|\mathbf{y}\setminus \left\{y\right\})=exp(w_e\Delta V_e(\mathbf{y}\to \mathbf{y}\setminus \left\{y\right\}))$, the Papangelou intensity obtained by considering the single energy term e. This allows viewing the contribution of each component. Moreover, we propose grouping these contributions into the data and prior contributions to, respectively, obtain $s_{\mathit{data}}$ and $s_{\mathit{prior}}$ such that the final score is a product of the two: $s\left(y\right|\mathbf{y}\setminus \left\{y\right\})=s_{\mathit{data}}\left(y\right)s_{\mathit{prior}}\left(y\right|\mathbf{y}\setminus \left\{y\right\})$.

In Figure 12, we illustrate how the two components of the score s can help analyze the results. Here we manually cluster the results into 4 basic categories: green objects correspond to detection with high prior and data scores (high s, top $25\%$), while blue detections have a higher data contribution ($s_{\mathit{prior}}<s_{\mathit{data}}$). The yellow detections correspond to objects with lower data scores ($s_{\mathit{prior}}>s_{\mathit{data}}$), often located in ambiguous locations. Purple objects have a low total score s (i.e., low confidence, bottom $15\%$). These detections are often pruned out when using a score threshold (as used in Figure 10 and Figure 11). This allows for identifying where the model is confident, either due to the data from the image or the interactions with nearby objects. One can access even more detailed information by looking into the individual scores $s_e$ for each energy term $V_e$.

4. Discussion

4.1. Quantitative and Qualitative Results

In Figure 10, we observe that despite noisy ground truth labels—a potential cause for over-constraining issues when using linear programming (see Section 2.5)—our models infer more regular configurations that align better with the physical constraints of the objects of interest. Sample 1 in Figure 10 BBA-Vec. struggles to properly align and detect all the objects. In sample 2, model YOLOV5-OBB shows patches of objects with no overlap regularization. Meanwhile, our model outputs regularized configurations even in these dense areas. Figure 10 (line 3) and Figure 11 showcase the robustness of our PP-based models when tested on noisy data or areas of limited information, such as shaded areas. The inference results on the ADS dataset, obtained with models trained on the ${\mathit{DOTA}}_{0.5}$ training set, highlight the model’s ability to generalize to unseen, albeit similar, data.

Comparison between ${\mathrm{CNN-PP}}_m$ and ${\mathrm{CNN-PP}}_w$ against CNN-LocalMax. in Figure 10 and Table 3 reveals the improvement brought by the PP-based approach, as the performance of the CNN alone (CNN-LocalMax.) falls short of our proposed method especially when confronted with noisy data.

Lastly, the comparable performances of ${\mathrm{CNN-PP}}_m$ and ${\mathrm{CNN-PP}}_w$, both qualitatively and quantitatively, highlight the efficacy of the CD method in inferring model parameters. This approach significantly reduces the need for manual parameter tuning (more specifically, prior parameters), where traditional PP methods, such as linear programming, fall short.

4.2. Limitations

Our method, while promising, faces several limitations that warrant consideration. First, the computational complexity involved in sampling a PP can be significant, particularly with dense configurations of objects. To address this, we could explore approximate sampling methods, although this might come at the expense of result accuracy, for instance, by implementing faster annealing techniques or discretizing the state space [48].

Second, while the model we present is tailored to a specific type of object, the underlying tools are adaptable to various object types. First, the PP framework allows for a wide diversity of objects to be modeled, thanks to the versatility of marks. Moreover, the energy formulation of our model allows for easy integration or removal of interaction models and priors as simple energy functions. For instance, we could incorporate classification as a mark in the Point Process model, utilizing features from the CNN or a prior law based on the geometric features of the objects (e.g., distinguishing trucks from cars based on length).

Lastly, our method relies on annotated data for training. However, energy models offer compositional flexibility; one potential avenue is to train the data terms of the model on annotated images while utilizing synthetic datasets (devoid of visual components) for prior terms training. This approach could help mitigate the reliance on annotated data and improve generalization to diverse scenarios.

5. Conclusions

In this paper, we present a method that integrates interaction models into object-detection algorithms, leveraging deep convolutional neural networks. We address challenges inherent in satellite imagery, such as atmospheric disturbances and limited resolution, by incorporating prior knowledge about object arrangements.

While CNN-based methods excel at pattern extraction, they struggle with learning object-to-object interactions without complex attention mechanisms. Conversely, Point Process methods offer a promising alternative, addressing both object likelihoods and arrangement coherence. However, existing approaches relying on contrast measures have limitations.

Our method combines CNN pattern extraction with the Point Process framework, simplifying contrast measure computations into efficient energy maps. By leveraging pre-computed potential maps, we enhance Point Process sampling methods, improving exploration efficiency within the state space.

We depart from conventional parameter estimation approaches for Point Process methods like linear programming and adopt Contrastive Divergence to estimate energy weights and internal parameters automatically. This improves accuracy and reduces reliance on manual parameter tuning.

To facilitate model comparison, we introduce a novel scoring function based on the Papangelou conditional intensity, providing a more comprehensive evaluation metric that takes into account interactions and allows for the explainability of results thanks to the easy decomposition of the energy model.

Experimental results on diverse datasets demonstrate our method’s effectiveness, particularly in noisy environments. Our model’s ability to generalize to new data showcases its robustness and potential for practical applications.

In conclusion, our method offers a promising approach to object detection, combining the strengths of CNNs with Point Process models. It shows potential for improving accuracy and efficiency in satellite imagery analysis or other applications where the priors are strong: e.g., object tracking (priors on dynamics) or road extraction (geometrical priors).