\section{Introduction}\label{sec:intro}
3D pose, shape and location reconstruction for hundreds of people in a large scene will help with modeling crowd behavior for simulation and security monitoring. However, no existing methods can achieve this with global consistency. In this paper, we aim to reconstruct the 3D poses, shapes and locations of hundreds of people in the global camera space from a single large-scene image, as shown in Fig.~\ref{fig:show}.

Although monocular human pose and shape estimation \cite{HMR, zhang2021pymaf} has been extensively explored over the past years, estimating global space locations together with human poses and shapes for multiple people from a single image is still a difficult problem due to the depth ambiguity. Existing methods \cite{CRMH, BEV} reconstruct 3D poses, shapes and relative positions of the reconstructed human meshes by assuming a constant focal length. But the
\YL{methods are limited to}
small scenes with a common FoV (Field of View). These methods cannot regress the people from a whole large-scene image \cite{PANDA} due to the relatively small and varying human scales in comparison to the image size.
Even with an image cropping strategy, these methods cannot obtain consistent reconstructions in the global camera space due to independent inference from the cropped images.
Besides, existing methods hardly consider the coherence of the reconstructed people with the outdoor scene, especially with the ground, since the ground is a common and significant element of outdoor scenes.
Taking a usual urban scene as an example, these methods may include wrong positions and rotations so that the reconstructed people
\YL{do not appear} to be standing or walking on the ground.

In general, there are three challenges in reconstructing hundreds of people with global consistency from a single large-scene image:
1) there are \YL{a large number} of people with relatively small and highly
\YL{varying}
2D scales;
2) due to the depth ambiguity from a single view, it is difficult to directly estimate absolute 3D positions and 3D poses of people in the large scene;
3) there is no large-scene image datasets with hundreds of people for supervising crowd reconstruction in large scenes.

In this paper, to address these challenges, we propose \emph{Crowd3D}, the first framework for crowd reconstruction from a single large-scene image.
To deal with the large number of people and various human scales, we propose an adaptive human-centric cropping scheme for a consistent scale proportion of people among different cropped images by leveraging the observation of pyramid-like changes in the scales of people in large-scene images.
To ensure the globally consistent spatial locations and coherence with the scene, we propose a progressive ground-guided reconstruction network Crowd3DNet to reconstruct globally consistent human body meshes from the cropped images by pre-estimating a global scene-level camera and a ground plane.
To alleviate the ambiguity brought \YL{in} by directly estimating absolute 3D locations from a single image,
we present a novel concept called \emph{Human-scene Virtual Interaction Point (HVIP)} for effectively converting the 3D crowd spatial localization problem into a progressive 2D pixel localization problem with intermediate supervisions. Benefiting from HVIP, our model can reconstruct the people with various poses including non-standing.

We also construct \emph{LargeCrowd}, a benchmark dataset with over 100K labeled humans (2D bounding boxes, 2D keypoints, 3D ground plane and  HVIPs) in 733 gigapixel images (19200×6480) of 9 different scenes. To our best knowledge, this is the first large-scene crowd dataset, which enables the training and evaluation on large-scene images with hundreds of people.
Experimental results demonstrate that our method achieves globally consistent crowd reconstruction in a large scene.
Fig.~\ref{fig:show} gives an example.
\emph{The code and dataset will be available for research purposes}.

To summarize, our main contributions include:
\begin{itemize}
\item [1)]
We propose Crowd3D,
a multi-person 3D pose, shape and location estimation framework for large-scale scenes with hundreds of people.
We design an adaptive human-centric cropping scheme and a joint local and global strategy to achieve the globally consistent reconstruction.
\item [2)]
We propose a progressive reconstruction network with \YL{the} newly defined HVIP, to alleviate the depth ambiguity and obtain global reconstructions in harmony with the scene.
\item [3)]
We contribute \emph{LargeCrowd}, a benchmark dataset with over 100K labeled crowded people in 733 gigapixel large-scene images (19200×6480),
which are valuable for the training and evaluation of crowd reconstruction and spatial reasoning in large scenes.
\end{itemize}

\section{Related Work}
\label{RelatedWork}
\textbf{Multi-person 3D Pose Estimation.} These methods can be divided into top-down \cite{3DMPPE,PANDANet,Lcr-net++,hmor} or bottom-up \cite{smap,mehta2020xnect,mehta2018single,fabbri2020compressed,kundu2020unsupervised,PAF} paradigms. The top-down methods first detect the people and then estimate the 3D pose of each person separately.
Moon \emph{et al}.~\cite{3DMPPE} estimate root location and root-relative pose separately after detecting the persons. They regard the area of 2D bounding box as a prior and adopt a neural network to learn a correction factor. HMOR \cite{hmor} divides human relations into three levels and formulates pair-wise ordinal relations in each level. Different from the top-down paradigm, the bottom-up methods directly detect all the joints and group them. However, most methods either optimize the translation in a post-processing way \cite{mehta2020xnect} or ignore the root localization. Inspired by monocular depth estimation methods, SMAP \cite{smap} utilizes \YL{a} deep convolutional neural network (CNN) to estimate a normalized root depth map and part relative-depth maps. The final root map is recovered with the given focal length and hence the camera parameters need to be known to obtain the absolute positions.

All the above methods only estimate 3D poses in the form of skeletons while missing shape information that is important for many applications, such as interpenetration reasoning to avoid impossible poses, person re-identification and crowd analysis.

\begin{figure*}[htb]
    \centering
    \includegraphics[width=0.98\linewidth]{materials/figure2_v60.pdf}
    \vspace{-0.2cm}
    \caption{Overview of Crowd3D framework. First, Crowd3D adopts an adaptive human-centric cropping scheme to crop the large-scene image into patches with hierarchical sizes for more appropriate scales of people. Then, Crowd3D estimates the scene-camera intrinsics and ground plane equation with human pose priors. Finally, Crowd3DNet takes the cropped image, the patch location, the estimated camera and ground parameters as inputs and outputs the crowd \YL{reconstruction} with consistent spatial locations in the global camera space.}
    \label{fig:overview}
    \vspace{-0.2cm}
\end{figure*}

\textbf{Multi-person 3D Pose and Shape Estimation.} Parametric human body models, \emph{e.g}., SMPL \cite{smpl}, have been widely adopted to represent the 3D pose and shape of a person.
Single-person 3D pose and shape estimation has been achieved with tremendous progress \cite{KEEPitSMPL,HMR,vibe,spec,kolotouros2019learning,moon2020pose2pose,pavlakos2019texturepose,sun2019human,zeng20203d,zhang2021pymaf}, while multi-person 3D pose and shape estimation still faces many challenges.

Some methods adopt a two-stage framework by utilizing a single-person reconstruction method for each detected person. 3DCrowdNet \cite{3DCrowdNet} leverages 2D poses to distinguish different people and uses a joint-based regressor to estimate human model parameters. This kind of approaches focuses more on the accuracy of pose and shape but ignores 3D spatial locations of the people which are important for holistic understanding of the scene.
To get coherent reconstruction results, Jiang \emph{et al}.~\cite{CRMH} propose CRMH, an R-CNN-based architecture, to detect all the people in an image and estimate their SMPL parameters by using an interpenetration loss and a depth ordering-aware loss in training.
This method calculates human depths based on the assumption that people are consistent in height, which will estimate excessive depths for short individuals like kids.
To solve the inherent body size and depth ambiguity problem, Ugrinovic \emph{et al}.~\cite{bodysize} propose a multi-stage optimization-based method to optimize the 3D translations and scales of body meshes estimated by CRMH \cite{CRMH}. Different from multi-stage methods with computation redundancy, BMP \cite{BMP} is a single-stage solution for multi-person mesh regression, which correlates the depth of a person with the features of different scales.  In ROMP \cite{ROMP}, the mesh and location information can be obtained in combination with the camera map and SMPL map according to the center map. However, this method is based on the assumption of weak perspective projection and can only reason about the 2D locations of people in the image plane. It uses an approximation method to obtain depth ordering.
To address this, BEV \cite{BEV} uses Bird's-Eye-View representation to simultaneously reason about body centers in image and in depth. As mentioned by itself \cite{BEV}, BEV is not trained or designed to deal with large ``crowds'' (\emph{e.g}., \YL{100s} of people) with a constant focal length assumption. In general, all the above methods can only get relative depths rather than absolute 3D positions, and they cannot be applied directly to large scenes.

\textbf{Multi-person Datasets.}
Multi-person datasets can be collected in  indoor controlled environments or in outdoor scenes.
Datasets such as \emph{Panoptic} \cite{Panoptic} build multi-view capture systems to obtain relatively high accurate ground-truths. For outdoor scenes, some datasets enable in-the-wild 3D capture from videos with IMUs \cite{3DPW} or from videos in which humans have to stay still \cite{smply}, while other datasets \cite{MPI} annotate in-the-wild images in 2D only. Besides, datasets such as \emph{Agora} \cite{Agora} and \emph{MuCo-3DHP} \cite{MuCo3DHP} generate synthetic images including 3D people and background images or 3D background scenes.
However, all the above datasets only contain a few people in small scenes.

In this paper, We propose the first work to reconstruct hundreds of people in a large scene with global consistency from a single RGB image.
We also contribute a benchmark dataset, \emph{LargeCrowd}, for the training and evaluation of crowd reconstruction in large scenes.

\section{Method}
\label{method}

Our work aims to recover a globally coherent reconstruction of crowd from a single large-scene image with hundreds of persons.
Fig.~\ref{fig:overview} shows the framework of our method.
The highlight of our method is that we design progressive position transform with our newly defined concept HVIP to establish a mapping between local image points and global spatial positions.
Our method consists of three main steps:
1) we adopt an adaptive human-centric cropping scheme (Sec. \ref{subsec:cropping}) to crop the large-scene image into patches with hierarchical sizes which
ensures that people in different cropped images have appropriate scales;
2) we estimate the camera intrinsics and ground plane equation (Sec. \ref{subsec:groundEstimation}) of the scene with human pose priors for subsequent inference;
3) taking the cropped images, ground plane and camera parameters as inputs, we design the Crowd3DNet (Sec. \ref{subsec:crowd3dnet}) with the progressive position transform based on HVIP (Sec. \ref{subsubsec:ground-guided}) to directly estimate the human meshes in the large-scene camera coordinate system.

\subsection{Adaptive Human-centric Cropping}

\label{subsec:cropping}
Instead of using uniform cropping \cite{ferreira2020s,li2022region} that cannot deal with people of various \YL{image} sizes,
we propose an adaptive human-centric cropping strategy to ensure that the height ratio between people and the corresponding cropped image is as consistent as possible among different cropped images.
It is crucial for accurate and reasonable estimation.
Inspired by the observation that
human heights hierarchically vary like a pyramid in the vertical direction of large-scene image,
the sizes of the cropped images \YL{should also} conform to a similar hierarchical change.
Heuristically,
we use \YL{a} geometric sequence to simulate the hierarchical change, which is simple but effective.
Define the heights of the persons at the top and the bottom of the large-scene image as $h_{\text{t}}$ and $h_{\text{b}}$, respectively. The upper and lower bounds of the image area \YL{to be processed} are defined as $b_u$ and $b_l$.
\YL{Considering} non-overlapping square blocks in \YL{the} vertical direction of image, \YL{we} represent the sizes of blocks from top to bottom \YL{as} $\{c_i\}_{i=1}^{n}$. When we set the height of people in \YL{a} block to be half of the block size and make $\{c_i\}_{i=1}^{n}$ comply with the rule of geometric sequence, we have $c_1=2 \times h_{\text{t}}$, $c_n=c_1 \times q^{n-1}$ and $\sum_{i=1}^{n}c_i=b_l-b_u$, where $q$ is the proportionality coefficient.
This cropping problem is formulated as:
\begin{equation}
    \label{eq:cropping}
     \mathop{\arg\min}\limits_{n, q}|c_n-2\times h_{\text{b}}|.
\end{equation}
For horizontal direction, the block size is consistent.
To ensure that each person can appear completely, we add new overlapped blocks between adjacent cropped images by combining the size of each half of the two cropped images.
The cropping parameters $h_{\text{t}}$, $h_{\text{b}}$, $b_u$, $b_l$ can be set manually or automatically. Details are given in the supplementary document.

\subsection{Camera and Ground Plane Estimation}

\label{subsec:groundEstimation}
We use the ground plane as a guidance for three reasons:
1) it is a common element in large scenes, especially surveillance scenarios;
2) it is the main object interacting with people in the large scenes, reflecting the harmony between people and the scene;
3) it provides the important global information to the local cropped images.

To estimate the ground plane equation and the scene-level camera parameters,
the pose prior of people can be used
\YL{for calibration}.
Note that the estimation of ground plane does not need too many people: more than ten people are enough as shown in the experiment (Sec. \ref{subsubsec:ablation}).
Besides, our method can reconstruct people with various poses, but at the current stage, we only consider the standing or walking people who can be regarded as vertical lines on the ground plane. These people are automatically selected from the 2D keypoints detection obtained from RMPE \cite{fang2017rmpe}.
We use a pinhole camera model with a focal length $f\ (f = f_x = f_y)$ where the principal point $(c_x, c_y)$ of the camera is the image center.
We represent the ground equation as $N^T P_{g} + D = 0$, where $N=(x_n, y_n, z_n)$ is the ground normal with $\Vert N \Vert_2 = 1$, $P_{g} \in \mathbb{R}^3$ is the point on the ground plane and $D$ is \YL{a} constant term.
For these standing people,
we define the midpoints of their two ankle keypoints as $P_a \in \mathbb{R}^3$ and the midpoints of two shoulder keypoints as $P_s \in \mathbb{R}^3$ .
The projections of $P_a$ and $P_s$ are $p_a = (u_a, v_a)$ and $p_s = (u_s, v_s)$, respectively.
\YL{Following} perspective projection,
we have $z_a \times \bar{p}_a=K P_a$, where $\bar{p} = (u, v, 1)^T$ represents the homogeneous \YL{coordinates} of $p = (u, v)$,  $K$ is the intrinsic matrix of the scene-level camera and $z_a$ is the depth of $P_a$.
Similar to \cite{physicalDistance}, we assume that $P_a$ is on the ground plane, and the line from $P_a$ to $P_s$ is parallel to the ground normal. We also set a fixed height prior $h$.
Therefore, we have $N^T P_a + D = 0$ resulting in
\begin{equation}
    \label{eq:ground-transform}

    z_a=- \frac{D}{N^TK^{-1}\bar{p}_a},
\end{equation}
and $P_s$ can be approximated by $P_s' = P_a + h\times N$. Then, the projection $\bar{p}_{s}'$ is computed by
\begin{equation}
    \label{eq:reprojection}
    {z_{s}'}\times\bar{p}_{s}' = {z_{s}'}\times
    \begin{bmatrix}
    u_{s}'\\
    v_{s}'\\
    1
    \end{bmatrix}
    =K (z_a \times K^{-1} \bar{p}_a + h \times N).
\end{equation}
To solve the camera and ground plane parameters $K$, $N$, $D$, we adopt the following optimization loss:
\begin{equation}
\begin{aligned}
    \label{eq:re-projectionLoss}
    L_{\text{param}} = & \lambda_{\text{angle}} L_{\text{cos}}( p_{s}' - p_{a}, p_{s} - p_{a} ) \\
    & + \lambda_{\text{mod}} \frac {\left|\Vert p_{s}' - p_{a} \Vert_2 - \Vert p_{s} - p_{a} \Vert_2\right|}
                        {\Vert p_{s} - p_{a} \Vert_2},
\end{aligned}
\end{equation}
where $L_{\text{cos}}$ is the cosine distance, and $\lambda_{\text{angle}}$ and $\lambda_{\text{mod}}$ are the weights of the corresponding loss terms.
Finally, we translate $0.1$ meters along the ground normal direction to get the ground plane where people stand on.

\subsection{Crowd3DNet}

\label{subsec:crowd3dnet}

As shown in Fig.~\ref{fig:overview}, Crowd3DNet is a one-stage multi-head network based on the body-center-guided representation \cite{ROMP}.
Different from previous methods \cite{ROMP,BMP,BEV}, we define a new concept, \emph{Human-scene Virtual Interaction Point
(HVIP)}, and a progressive position transform (Sec \ref{subsubsec:ground-guided}) to better infer the global 3D positions of people. Crowd3DNet outputs four maps including a body center heatmap, a torso center offset map, a 2D HVIP map and a SMPL parameters map.
The body center heatmap predicts the probability that each location is the center of a human body.
If the body center heatmap gives positive responses, the network samples relevant parameters from other maps at the corresponding center locations to obtain 2D torso center offsets, 2D HVIPs and SMPL parameters of people.
With progressive position transform based on HVIP, Crowd3DNet combines the sampled parameters, the input of ground plane equation and scene-level camera parameters to infer accurate 3D positions of people, \YL{achieving multi-person reconstruction} in the large-scene camera system from the cropped images.

\subsubsection{Progressive Position Transform Based on HVIP}
\label{subsubsec:ground-guided}
We design the progressive position transform based on Human-scene Virtual Interaction Point (HVIP) to help infer the accurate 3D locations of persons in the large-scene camera system.
The core idea is to infer the global 3D position from 2D image pixel points by HVIP and ground transform to avoid the depth ambiguity of estimating from a single view directly.
We define the HVIP which represents the projection point of a person's 3D torso center on the ground plane in the global camera space, marked as $P_{v}=(x_{v}, y_{v}, z_{v})$.
The torso center is a semantic point on human body, \emph{i.e}., the center of two shoulder and two hip joints,
represented as $P_{t}=(x_{t}, y_{t}, z_{t})$.
As show in Fig.~\ref{fig:hvip},
HVIP is a point on the ground plane, which can participate in ground transform directly to establish the mapping from image pixels to 3D points on the ground plane.
HVIP binds a person body's semantic point but it is not on the human body itself. Therefore, different from previous method \cite{bodysize} that forces people's ankle joints \YL{to be} on the ground, \YL{which limits} the posture of people, HVIP \YL{is determined by} the 3D space position of people and can deal with people in various postures.
Because the line from $P_{t}$ to $P_{v}$ is perpendicular to the ground, we have $P_{t}=P_{v}+d\times N$, where $d$ represents the distance from $P_{t}$ to the ground plane.
\begin{figure}[!t]
    \centering
    \includegraphics[width=0.98\linewidth]{materials/fig3.pdf}
    \vspace{-0.2cm}
    \caption{Progressive position transform based on HVIP.}
    \label{fig:hvip}
    \vspace{-0.4cm}
\end{figure}
We represent the projection points of $P_{v}$ and $P_{t}$ with $p_{v}=(u_{v}, v_{v})$ and $p_{t}=(u_{t}, v_{t})$, respectively.
Refer to Eq. \eqref{eq:ground-transform} and Eq. \eqref{eq:reprojection}, we deduce
\vspace{-0.2cm}
\begin{gather}
    \label{eq:tc}
    P_{v}=-\frac{D}{N^TK^{-1}\bar{p}_{v}} \times K^{-1}\bar{p}_{v}, \\
    d=\frac{f\times y_{v} - (v_{t}-c_y)\times z_{v}}{(v_{t}-c_y)\times z_{n} - f \times y_{n}}.
\end{gather}
Therefore, when the network predicts $p_{v}$ and $p_{t}$, $P_{t}$ can be uniquely determined.
Then, the predicted body mesh $M_{cam}$ in the global camera space follows $M_{cam}=M-P_{t-smpl}+P_{t}$, where $M \in \mathbb{R}^{6890 \times 3}$ and $P_{t-smpl} \in \mathbb{R}^{3}$ are the predicted vertices and torso center in SMPL \cite{smpl} space, respectively.
Finally, considering the cropping, our network takes the cropped image as input and predicts the local $p_{v-local}$ and $p_{t-local}$ on the cropped image. We have $p_{v}=p_{v-local}+t_{crop}$ and $p_{t}=p_{t-local}+t_{crop}$, where $t_{crop}$ represents the pixel \YL{coordinates} of the upper left corner of the cropped image.
Our progressive position transform with HVIP builds a mapping from some pixel points on \YL{the} cropped image to the global 3D location, which can simply and effectively predict the precise crowd positions.

\vspace{-0.3cm}
\subsubsection{Representations}
\vspace{-0.1cm}

\noindent\textbf{Input Parameters.}
The network cannot perceive the whole scene information only from the cropped image, hence we take the estimated ground and camera parameters as extra inputs.
We define the camera input \YL{as} $(\frac{f}{W_s}, \frac{\hat{c}_{x} - c_x}{c}, \frac{\hat{c}_{y} - c_y}{c})$ which includes the information of FOV of scene and the principal point shift, where $W_s$, $(\hat{c}_{x}, \hat{c}_{y})$ and $c$ are the width of large-scene image, the image center of \YL{the} cropped image and the size of \YL{the} cropped image, respectively.

\noindent\textbf{Body Center Heatmap.}
The body center heatmap $\boldsymbol{C_m}$ represents the body center likelihood by a Gaussian kernel combining body scales, where $C_m \in \mathbb{R}^{1 \times H \times W}$ and $H=W=64$. We define the body center \YL{the} same as \cite{ROMP}.

\noindent\textbf{Torso Center Offset Map.}
Although we can directly define the body center as the 2D torso center, in practice, the body center heatmap tends to find a person's body salient point, especially when the person is occluded. Therefore, it is necessary to predict the human torso center separately. The torso center offset map $T_m \in \mathbb{R}^{2 \times H \times W}$ contains the offset between 2D torso center $p_{t-local}$ and body center.

\noindent\textbf{2D HVIP Map.}
The goal of 2D HVIP map $H_m \in \mathbb{R}^{1 \times H \times W}$ is to obtain the 2D HVIP projection $p_{v-local}$ on the cropped image.
The line from $P_{v}$ to $P_{t}$ is parallel to the ground normal, following the perspective theory, we have the projection points $p_{v}$, $p_{t}$ and the vanishing point of the ground normal $p_{vp}$ are collinear on image, where $p_{vp}=KN$. Therefore, we only need to estimate the 1D length from $p_t$ to $p_v$ to obtain 2D HVIPs.

\noindent\textbf{SMPL Map.}
The SMPL map $S_m \in \mathbb{R}^{145 \times H \times W}$ includes the parameters of SMPL \cite{smpl} of people and a small 3D offset $\delta t$.
The SMPL parametric model can represent various shape and pose with a \YL{small} number of parameters.
It takes the pose parameters $\theta$ and the shape parameters $\beta$ as inputs and outputs a body mesh $M \in \mathbb{R}^{6890 \times 3}$. We adopt the 6D rotation representation \cite{zhou2019continuity} and drop the last two hand joints. \YL{Considering} the error of the dataset annotations, we predict \YL{an} offset $\delta t$ to further refine the position of people by $P_t=P_v+ d \times N + \delta t$.

\vspace{-0.2cm}
\subsubsection{Loss Function}
\label{subsec:loss functions}

Crowd3DNet is supervised by the weighted sum of multiple loss \YL{terms} as follows:
\vspace{-0.2cm}
\begin{equation}
\begin{aligned}
    L=& \lambda_{\text{center}}L_{\text{center}}+
        \lambda_{\text{mesh}}L_{\text{mesh}}+
        \lambda_{\text{hvip}}L_{\text{hvip}} \\
        &+
        \lambda_{\text{tc}}L_{\text{tc}} +
        \lambda_{\text{root}}L_{\text{root}}+
        \lambda_{\text{gn}}L_{\text{gn}}+
        \lambda_{\text{out}}L_{\text{out}},
\end{aligned}
\end{equation}
\begin{equation}
\begin{aligned}
    L_{\text{mesh}}=& \lambda_{\text{pose}}L_{\text{pose}}+
        \lambda_{\text{shape}}L_{\text{shape}}+
        \lambda_{\text{j2D}}L_{\text{j2D}} \\
        &+
        \lambda_{\text{j3D}}L_{\text{j3D}} +
        \lambda_{\text{paj3D}}L_{\text{paj3D}}+
        \lambda_{\text{gm}}L_{\text{gm}},
\end{aligned}
\end{equation}
where $L_{\text{center}}$ is
the 2D focal loss \cite{lin2017focal}, and $L_{\text{mesh}}$ is the common SMPL related $L_2$ loss including pose parameter loss $L_{\text{pose}}$, shape parameter loss $L_{\text{shape}}$, 2D joint projection loss $L_{\text{j2D}}$, 3D joint loss $L_{\text{j3D}}$ and 3D joint loss after Procrustes alignment $L_{\text{paj3D}}$.
$L_{\text{hvip}}$, $L_{\text{tc}}$ and $L_{\text{root}}$ are all $L_2$ losses, which are used to supervise 2D HVIP projection, 2D torso center and absolute root position, respectively.
$L_{\text{gn}} = L_{\text{cos}}(P_s - P_a, N)$ is a ground normal regularization term to enhance the interaction consistency between people and ground plane, where $P_s - P_a$ is the approximated craniocaudal direction of human.
We also use an out-of-bound loss to prevent people from penetrating the ground.
More concretely, we use $L_1$ loss to punish the point with the most serious penetration into the ground plane, and the out-of-bound loss is defined as
\vspace{-0.2cm}
\begin{equation}
    L_{\text{out}}=\left|\min(\{ \bar{v}_i \cdot G \mid \bar{v}_i \cdot G<0 \})\right|,
\end{equation}
where $v_i \in M_{\text{cam}}$ and $G=[N^T, D]^T$.

\subsection{Scene-specific Optimization and Merging}
\label{subsec:finetune}

To improve the harmony between reconstructed crowd and scene, and the generalization to various camera and ground plane parameters, we add a scene-specific optimization for a new scene at test time.
Please note that the scene-specific optimization is performed only once for a camera-fixed scene, \emph{i.e}., only one image of the scene is \YL{needed}.
Specifically, given a new scene at test time, we optimize a small set of weights in the head layer of Crowd3DNet with the ground normal and 2D poses
\YL{estimated} in the camera and ground plane module. The optimization loss $L_{\text{opt}}$ is
\vspace{-0.1cm}
\begin{equation}
    L_{\text{opt}} =
    \lambda_{\text{j2D}}L_{\text{j2D}} +
    \lambda_{\text{gm}}L_{\text{gm}} +
    \lambda_{\text{gn}} L_{\text{gn}} +
    \lambda_{\text{out}}L_{\text{out}}.
\end{equation}

We finally remove \YL{duplicated} persons in the overlapped adjacent patches by merging. The merging operation retains the people farther away from the boundary of the overlapped region, which tends to keep more complete people to avoid truncation.

\section{Experiments}

\subsection{Large-scene Crowd Dataset}

To train and evaluate crowd reconstruction in a large scene, we contribute \emph{LargeCrowd}, which is a benchmark dataset with over 100K labeled humans in 733 gigapixel images ($19200\times6480$) of 9 different scenes (5 scenes for training and 4 scenes for \YL{testing}).
The images are extracted at a minimum interval of $3s$ from gigapixel streams which are captured by a ZoheTec JMC315 array camera.
We annotate the bounding boxes, 2D poses and 2D HVIPs of all the visible people in the images, with the maximum error less than 5 pixels for $95\%$ labels. We measure 3D landmarks in a world coordinate system and label the corresponding 2D points to solve the camera extrinsic matrix for each scene.
Then, we compute the homography matrix for each ground plane. The homography matrices with labeled ground segmentations and HVIPs provide the true physical positions of the persons.

\begin{table}[!t]

    \caption{Comparison on \emph{LargeCrowd} dataset.}
    \vspace{-0.2cm}
    \centering
    \small
    \begin{tabular}{@{}c|cccc@{}}
    \toprule
    Method  &PPDS$\uparrow$  &PA-PPDS$\uparrow$    & PCOD$\uparrow$      &OKS$\uparrow$   \\
    \midrule
    SMAP\cite{smap}-Large  &60.07     & 61.79                  & 72.18        & 63.16         \\
    CRMH \cite{CRMH}-Large  &55.01      & 64.92                    & 79.42         & 62.64         \\

    BEV \cite{BEV}-Large  & 74.81      & 75.81           & 87.93         & 67.08         \\

    \hline
    Crowd3D w/o HVIP            & 81.02          & 89.63        & \textbf{93.19}           & 64.87        \\
    Crowd3D    &\textbf{81.75}       & \textbf{89.70}                   & 93.03        & \textbf{69.27}        \\

    \bottomrule
    \end{tabular}
\label{tab_large}
\vspace{-0.4cm}
\end{table}

\subsection{Implementation Details}
We use HRNet-32 \cite{cheng2020higherhrnet} as backbone, each head of which is composed of two ResNet \cite{he2016deep} blocks with batch normalization.
We resize input images to $512 \times 512$ with zero padding to keep the same aspect ratio. We also use the collision-aware representation of ROMP \cite{ROMP} to push apart close body centers.
Our training process has two stages: 1) start by training the body center heatmap, torso center offset map and 2D HVIP map for 20 epochs to make sure that the subsequent learning about body mesh has a suitable initial position; 2) train the full model with all losses for 70 \YL{epochs}. We implement Crowd3DNet with PyTorch and adopt the Adam \cite{adam} as optimizer with 5e-5 learning rate. We train our model on \emph{LargeCrowd}, \emph{Agora} \cite{Agora}, \emph{MuCo-3DHP} \cite{Muco} and a single person dataset \emph{Human3.6M} \cite{Human3.6M}.

\subsection{Evaluation Metrics}
We use the torso center as the location of people and evaluate the location distribution of crowd by the distances between people.
We define a metric called pair-wise percentual distance similarity (PPDS) as
\vspace{-0.2cm}
\begin{gather}
    \textit{PPDS}=
    \frac{\sum_{k=1}^{n-1} \sum_{i=k+1}^{n} 1- \min{(d_{ik}, 1)}}{C_n^{2}}, \\
    d_{ik}=\left | \frac{ \Vert E_k - E_i \Vert -\Vert G_k - G_i \Vert }{\Vert G_k - G_i \Vert}\right |,
\end{gather}
where $n$ is the number of people in the image, and $E_i$ and $G_i$ represent the estimated and ground-truth locations of the $i$-th person, respectively.
To evaluate the relative crowd distribution, we also define the procrustes-aligned pair-wise percentual distance similarity (PA-PPDS) which aligns the reconstructed crowd and the ground truth by Procrustes alignment to exclude the influence of scale and rotation.
Due to the lack of 3D pose annotations, we use the object keypoint similarity (OKS) \cite{coco} to evaluate the 2D poses. The percentage of correct ordinal depth (PCOD) \cite{smap} is used to evaluate the ordinal depth relations between all pairs of people in the image.

\begin{figure*}[!h]
    \centering
    \includegraphics[width=0.8\linewidth]{materials/fig4_v11.pdf}
    \vspace{-0.2cm}
    \caption{Qualitative results on \emph{LargeCrowd}. The same color or number corresponds to the same person, and gray indicates that the person is not matched.}
    \label{fig:comparision}
      \vspace{-0.4cm}
\end{figure*}

\vspace{-0.2cm}
\subsection{Comparison}
\vspace{-0.1cm}
Because no existing methods can directly handle large-scene images with hundreds of people, we compare our method with three baselines that are modified \YL{from} the state-of-the-art methods: SMAP \cite{smap}, CRMH \cite{CRMH}, and BEV \cite{BEV}. We denote these baselines as SMAP-Large, CRMH-Large and BEV-Large.
Specifically, we first use our adaptive human-centric cropping to obtain the hierarchical cropped images as their inputs and infer the respective reconstructed results on the cropped images. To obtain the global reconstruction results for these methods, we provide the scene-camera intrinsics estimated by our method to them.
For CRMH \cite{CRMH} which predicts human \YL{bodies} in bounding \YL{boxes} by \YL{a} weak perspective camera model, we use its transform from bounding box position to depth of full image to infer the predicted locations in the global camera space.
Both SMAP \cite{smap} and BEV \cite{BEV} infer the locations through perspective camera models. Following previous method \cite{albiero2021img2pose}, we scale the depths of their results according to the focal length of \YL{the} scene. Please refer to  supplementary material for more details. For fair comparison, we fine-tune all the compared methods on \emph{LargeCrowd}.
Table \ref{tab_large} gives the quantitative results.
Our method outperforms other approaches in terms of all the metrics. Especially, the obvious advantage in PPDS, PA-PPDS and PCOD shows that our method can  predict accurate crowd location distribution, including physical distances and relative arrangements.
Fig.~\ref{fig:comparision} shows qualitative comparison results.
The complete bird's-eye view on the right shows that our predicted crowd distribution is consistent with the input image, while the compared methods are not consistent.
Taking the persons \YL{labeled} with numbers for example, only our method recovers correct relative positions.
The reconstructed people by the existing methods independently inferred from the cropped images are inconsistent in the global large-scene camera space.
Besides, although these methods show reasonable projection results, the wrong global positions mean that their predicted 3D human bodies have wrong scales. We also provide comparison results on public small-scene datasets in supplementary material.

\vspace{-0.1cm}
\subsection{Ablation Study}
\vspace{-0.1cm}
\label{subsubsec:ablation}
\noindent\textbf{Impact of the Number of People on the Estimated Ground and Camera.}
We explore the impact of the number of people on estimating ground and camera parameters by controlling the number of people used in optimization, and the newly added people are randomly selected.
The metrics include a cosine distance for ground normal and a root mean square error for focal length. We calibrate the camera and obtain the ground-truth focal length (about 27000).
As shown in Fig.~\ref{fig:personNum}, more than ten people are enough for estimating ground normal, which is common in real-world large-scale scenes, especially surveillance scenarios with hundreds of people. The focal length is not sensitive to the number of people.

\noindent\textbf{Progressive Position Transform Based on HVIP.}
Our progressive position transform based on HVIP effectively helps the network to predict accurate global 3D positions of people.
To verify this, we compare our full model with a variant of Crowd3DNet, Crowd3D w/o HVIP, which adopts the midpoint of ankle joints to participate in ground transform directly without using HVIP. The result is shown in Table~\ref{tab_large}.
Benefiting from HVIP, which makes use of the ground plane without restricting human posture, Crowd3DNet has obvious advantages on OKS.

\noindent\textbf{Adaptive Human-centric Cropping.}
To verify the adaptive human-centric cropping scheme, we denote a metric called cropping score, which counts the ratio of people with the appropriate scale after cropping.
For a person with appropriate scale, we set the ratio of his height to the corresponding cropped image within [0.3, 0.8], and he is not truncated.
The comparison result on \emph{LargeCrowd} between adaptive human-centric cropping and uniform cropping is 0.923 \emph{vs.} 0.806, which demonstrates the effectiveness of our adaptive human-centric cropping scheme.

\begin{figure}[!t]
    \centering
    \includegraphics[width=0.8\linewidth]{materials/person_num.png}
    \vspace{-0.2cm}
    \caption{The impact of the number of people on the camera and ground plane estimation.}
    \label{fig:personNum}
    \vspace{-0.4cm}
\end{figure}

\vspace{-0.2cm}
\section{Conclusion and Discussion}
\vspace{-0.1cm}
\label{Sec:Conclusion}

\noindent\textbf{Conclusion.}
We propose Crowd3D to reconstruct hundreds of people with global consistency from a single RGB large-scene image. Our method is a joint local and global inference framework which converts the complex crowd localization into pixel localization by our defined HVIP concept and the parameters of pre-estimated scene-level camera and ground plane.
Our adaptive human-centric cropping scheme and progressive position transform based on HVIP solve the challenges of large number of people, large variations in human scale and complex spatial distribution in large scenes.
We also contribute a large-scene dataset called \emph{LargeCrowd} to help train and evaluate crowd reconstruction in large scenes with hundreds of people.
Experimental results demonstrate that our method can achieve globally consistent crowd reconstruction in large scenes.

\noindent\textbf{Limitations and Future Work.}
We focus on outdoor real-world large-scale scenes which contain one or several ground planes.
Our method may be easily extended to multi-ground scenes by using the existing image-based ground plane segmentation methods or manual segmentation, which is taken as our future work.
Although our \YL{Crowd3D} shows effective crowd reconstruction in a global camera space, there are still some cases that we cannot solve well, \emph{e.g}., the people in complex ground conditions and the persons with complicated postures or severe occlusions.
In future work, we will focus on a wider range of large-scale scenes with complex ground and crowd environments.

{\small
\bibliographystyle{ieee_fullname}


\section{Introduction}

With the rapid increase of computational resources available in
computational centers, the numerical representation of flows around
complex configurations of industrial interest is becoming an established
reality. A large effort by the scientific community in the last
decades has allowed to identify strategies for the numerical
investigation of such problems. As a result, new methods studied to
solve problematic aspects of the simulations are getting to maturity. Among
these, the Immersed Boundary Method (IBM)
\cite{Peskin1972,Peskin2002,Mittal_2005,Sotiropoulos_2014,Kim2019_ijhff} is an increasing popular
tool for the representation of moving bodies. In the IBM, the
presence of the immersed body is represented via the inclusion of
source terms in the dynamic equations of the physical model. The
usage of the IBM allows to naturally exclude problematic aspects of
the mesh procedure, such as deformed / stretched mesh elements for
the representation of complex geometries. These aspects can result
in very large numerical errors, significantly affecting the
prediction of the flow features. In addition, IBM can naturally
account for the immersed body motion / deformation with almost zero
computational resources, excluding expensive re-meshing procedures.

While many different IBM proposals can be found in the literature
for purely dynamic effects associated with the immersed body
\cite{Beyer1992,Goldstein_jcp_1993,Fadlun2000,uhlmann_jcp_2005,Cheny2010,Kumar_2020}
(i.e. focusing on the determination of a source forcing term for the
momentum equation) the analysis of heat transfer via IBM is less
investigated and mainly dealing with incompressible flow simulation
\cite{Kim2001,Ren2013,Park2017}.  {Recent works on the subject have
been proposed by Luo et al. \cite{Luo2015,Luo2016} on fixed bodies. Also, simulations of compressible flows in interaction with moving solids are still rare (see for example Qu et al. \cite{Qu_2018} and Khalili et al. \cite{Khalili_2018}) and very few contributions can be found including the coupling with heat exchange effects \cite{Luo2016}.}	 Accounting for
heat transfer effects in IBM applications for compressible flows
will provide significant advances for a number of relevant high
speed applications, such as the flow in aerospace engines. In fact,
many features of the flow such as the behavior of the boundary layer
of coherent structure organization in turbulent wakes is affected by
thermal phenomena occurring between the flow and the surface of the
immersed body. The correct representation of such phenomena is
essential for an accurate prediction of heat dissipation in complex
mechanical systems.

In the present work, the IBM tool recently developed for integration in the opensource platform \textit{OpenFOAM} \cite{Constant_2017,Riahi_2018} is extended to the analysis of heat exchange for compressible flows. To do so, an explicit term is included in the energy equation of the dynamic system. The article is structured as follows. In section \ref{sec:numIBM} numerical details about the dynamic system and the IBM method are provided and the implementation of the IBM method within the numerical solver is detailed. In section \ref{sec::Validation} the research tool is validated with the analysis of academic test cases. In section \ref{sec::Sphere}, a three-dimensional application is analyzed, considering the flow around a heated sphere. In section \ref{sec:FSI}, elements of fluid-structure interaction are discussed via the analysis of an oscillating cylinder. Finally, concluding remarks are drawn in section \ref{sec:conclusion}.

\section{Numerical ingredients and Immersed Boundary Method}\label{sec:numIBM}

The IBM algorithm developed for heat transfer prediction is here
introduced and its analytic derivation is described.
\subsection{Governing equations}
The starting point is represented by the general Navier--Stokes
equations for a compressible fluid:
\begin{eqnarray}
\label{eq:Mass}
\frac{\partial \rho}{\partial t} &+& div(\rho \mathbf{u}) =0 \\
\label{eq:Momentum}
\dfrac{\partial\rho \mathbf{u}}{\partial t} &+& \mathbf{div}(\rho \mathbf{u} \otimes \mathbf{u}) = - \mathbf{grad} p + \mathbf{div} \overline{\overline{\tau}} + \mathbf{f} \\
\label{eq:Energy}
\dfrac{\partial\rho (e +\mathbf{u} \cdot \mathbf{u}/2)}{\partial t} &+& div (\rho \,(e+\mathbf{u} \cdot \mathbf{u}/2) \, \mathbf{u}) = \nonumber \\
 &-& div( p \, \mathbf{u}) + div (\overline{\overline{\tau}} \mathbf{u}) + div (\lambda(\theta) \mathbf{grad} \theta)+ \mathbf{f} \cdot \mathbf{u} + q
\end{eqnarray}
where $\rho$ is the density, $p$ is the pressure, $\mathbf{u}$ is the velocity, $\overline{\overline{\tau}}$ is the tensor of the viscous constraints, $e$ is the internal energy, $\lambda$ is the thermal conductivity, $\theta$ is the temperature and $\mathbf{f}$ is a general force term.
In the case Newtonian fluids are considered, the tensor $\overline{\overline{\tau}}$ can be written as:

\begin{eqnarray}
\overline{\overline{\tau}}=\mu(\theta) \left( (\overline{grad}\mathbf{u}+^{t}\overline{grad}\mathbf{u})- \frac{2}{3} div(\mathbf{u}) \right)
\end{eqnarray}
where $\mu$ is the dynamic viscosity. Its value is calculated using
the Sutherland's law as a function of temperature $\theta$. $q$ is a
general source term for the energy equation. This system is closed using the equation
of state for perfect gas $p=\rho r \theta$ where $r$ is the specific
gas constant.

\subsection{Immersed Boundary Method for compressible flows: numerical formulation}
\label{sec:IBMcompressible} The IBM proposed in this work roots in
recent proposals by Riahi et al. \cite{Riahi_2018} which was
successfully used for the analysis of compressible flows around
immersed bodies with adiabatic walls. The IBM strategy relies on
communication between an Eulerian mesh used for calculation of the
flow and Lagrangian markers representing the discretized shape of
the immersed geometry. The communication between the two frames of reference allows for the calculation of the dynamic effects acting on the body surface in the Lagrangian space and for their representation on the Eulerian mesh as source terms.

The first step is the interpolation of physical quantities calculated on the Eulerian mesh to the Lagrangian markers. Lowcase type is used for information on the Eulerian setting, while capital letters (or via the subscript $L$ for Greek letters) are used to indicate physical quantities on the Lagrangian markers. The interpolation operator $\mathcal{I}$ used in this work is:

\begin{equation}
\label{eq:interp}
\mathcal{I}[\mathbf{p}]_{X_s} = [\mathbf{P}](X_s) = \sum_{j\in D_s} (\mathbf{p})_{j}^n \delta_h(\mathbf{x}_{j} - \mathbf{X}_s) \Delta \mathbf{x}
\end{equation}

where $\mathbf{p}$ and $\mathbf{P}$ are a physical quantity, scalar or vectorial ($\theta$ or $\rho \mathbf{u}$ for example) and $\mathbf{x}$ and $\mathbf{X}$ represent the physical coordinates. For each Lagrangian marker $s$,  $D_s$ represents the set of points of the Eulerian mesh from which information is extracted (computational stencil) and the interpolation kernel $\delta_h$ is the discretized delta function used in Pinelli et al. \cite{Pinelli_JCP_2010}:
 {
\begin{equation}
\delta_h(r) = \begin{cases}
\frac{1}{3} \left( 1 + \sqrt{1 - 3 r^2} \right) \quad 0 \leq r \leq 0.5 \\
\frac{1}{6} \left(5 - 3r - \sqrt{1 - 3 (1 -r)^2} \right) \quad 0.5 \leq r \leq 1.5 \\
0 \quad \text{otherwise}
\end{cases}
\end{equation}
}

 { $\Delta \mathbf{x}$ refers to an Eulerian quadrature, which is $\Delta \mathbf{x}=\Delta x \Delta y \Delta z $ for a Cartesian mesh with uniform distribution.} This step is essential for the calculation of the dynamic effects on the immersed body surface. These effects, which will be detailed in section \ref{sec:forcage_energie1}, are represented as source terms for both the momentum equation ($\mathbf{F}$) and the energy equation ($Q$).

The second step of the IBM procedure is the spreading of the volume force calculated in the first step to the Eulerian mesh elements. We will refer to this term as $\mathbf{R}$, knowing that the numerical procedure for $\mathbf{F}$ and $Q$ is the same. The value of the source term evaluated on the Eulerian mesh, $\mathbf{r}(\mathbf{x}_j)$, is given by:
\begin{equation}
\label{eq:force2}
\mathbf{r}(\mathbf{x}_j) = \sum_{k \in D_j} \mathbf{R}_k \delta_h(\mathbf{x}_{j} - \mathbf{X}_k) \boldsymbol{\epsilon}_k
\end{equation}
 {The $k$-index controls a loop over the Lagrangian markers whose support contains the Eulerian node $j$.} $\boldsymbol{\epsilon}_k$ is the Lagrangian quadrature, which is calculated by solving a linear system to satisfy a partition of unity condition.  {This system can be written in the form}
 {
\begin{equation}
A \mathbf{\epsilon} = \mathbf{1}
\end{equation}
}

 {where $\mathbf{\epsilon} = (\epsilon_1, \cdots, \epsilon_{N_s})^t$ $\mathbf{1} =(1, \cdots, 1)^t$ are vectors whose size is equal to the number of Lagrangian markers $N_s$.
On the other hand, the component $A_{kl}$ of the matrix $A$ is the product between the $k^{th}$ and the $l^{th}$ interpolation kernels:}

 {
\begin{equation}
A_{kl} = \sum_{j \in D_l} \delta_h(x_j-X_k) \, \delta_h(x_j-X_l)
\end{equation}
}

The calculation of the dynamic effects associated with the presence of the immersed body in step one is thus crucial for obtaining a precise representation of the flow. The research team \cite{Riahi_2018} was able to derive an expression for the source term of the momentum equation $\mathbf{F}$ which takes into account the iterative nature of numerical solvers used for the calculation:

\begin{equation}
\label{eq:forceMomentum}
\mathbf{F} = a_s \left( (\rho_L \mathbf{U})^{d} - \overline{(\rho_L \mathbf{U})} \right)
\end{equation}

where $a_s$ is a coefficient derived by the discretization procedure.
The superscript $d$ over the physical quantities represents the target behavior of the flow (i.e. zero velocity at the wall) while the overbar represents the result of the interpolation of the quantities from the Eulerian mesh performed in step one.

The novelty of the present approach is represented by the inclusion
of a volume source term $q$ for the energy equation, which mimics
heat transfer phenomena on the body surface. These effects are
derived from the discretized set of equations via analytic
manipulation. In the present analysis, the heat transfer will be
considered equivalent to a boundary condition applied on the surface
of the immersed body. The analysis will be restricted to the case of
\textit{imposed temperature at the wall}, which is equivalent to
imposing a Dirichlet condition for the temperature field.

\subsubsection{Imposed temperature on the body surface - analytic derivation}
\label{sec:forcage_energie1}
Similarly to previous works
\cite{uhlmann_jcp_2005,Riahi_2018}, the starting point is
represented by the discretized set of equation in the Eulerian frame
of reference for the mesh element of coordinate $x$. Time advancement from the
instant $n$ to $n+1$ is considered. For sake of simplicity, the procedure for the momentum equation and the forcing $\mathbf{F}$ (detailed in \cite{Riahi_2018}) is not reported here, and the development is provided only for the internal energy equation \ref{eq:Energy}. Its discretized form for an iterative solver can be written as:
\begin{equation}
\label{eq:corr3.1}
a_e \, (\rho e)^{n+1} = {\phi}_e^{n+1/2}  + f_e^{n+1/2} + q_e^{n+1/2}
\end{equation}
with:
\begin{itemize}
\item $a_e$ : coefficient derived by the discretization procedure
\item ${\phi}_e$ : discretization term calculated at the intermediate
time $n+\dfrac{1}{2}$ which includes the contribution of the
convection term, the heat flux, the pressure work and the viscous
friction work.
\item $ f_e^{n+1/2} $ : discretized form of the term $\mathbf{f} \cdot \mathbf{u}$
\item $ q_e^{n+1/2} $ : discretization of the volume source term $q$. This term represents the heat exchange
between the flow and the immersed body, calculated at the
intermediate time $n+\dfrac{1}{2}$
\end{itemize}

One must keep in mind that $\mathbf{f}$ represent the discretized
volume source term representing the effect of the body in the
momentum equation, as in classical IBM strategies. Thus, the term
$q_e$ is its correspondent in the energy equation, which describes
the heat transfer between the flow and the immersed body. Thus, if
one targets a behavior of the internal energy $(\rho e)^{d}$ at the
instant $n+1$, then the optimized value of the volume source term
$q_e$ must comply with the equation:
\begin{equation}
\label{eq:corr3.3}
q_e^{n+1/2} =a_e \, (\rho e)^{d} - {\phi}_e^{n+1/2} - f_e^{n+1/2}
\end{equation}

By applying the Eulerian-Lagrangian transformation via the
interpolation operator proposed by Pinelli et
al.\cite{Pinelli_JCP_2010}, the equation  \ref{eq:corr3.3} on the Lagrangian markers representing the body surface in the IBM method is
transformed in:
\begin{equation}
\label{eq:corr3.4}
{Q}_e^{n+1/2}=a_e \, (\rho_L E)^{d} - {\Phi}_e^{n+1/2} - F_e^{n+1/2}
\end{equation}
\begin{itemize}
\item $(\rho_L E)^d$ : target internal energy value on the Lagrangian markers
\item ${\Phi}_e^{n+1/2}$ : total contribution of the convection term, the pressure
work term, the viscous friction work term and heat flux term
interpolated on the Lagrangian marker at the time $n+\dfrac{1}{2}$
\end{itemize}

Let us now consider again equation \ref{eq:corr3.1} in the case of the
presence of the body surface does not introduce direct thermal interactions with the flow. This state will be referred to as neutral behavior. In this case, the time advancement to $n+1$ (a suffix $l$ is here used) reads as:
\begin{equation}
\label{eq:corr3.5}
a_e \, (\rho e)^{l} = {\phi}_e^{n+1/2} +f_e^{n+1/2}
\end{equation}

If one further performs a transformation from the Eulerian system to
the Lagrangian frame of reference, equation \ref{eq:corr3.5}
becomes:
\begin{equation}
\label{eq:corr3.6}
a_e \, (\rho_L E)^{l} = {\Phi}_e^{n+1/2} + F_e^{n+1/2}
\end{equation}
where $(\rho_L E)^{l}$ represents the internal energy at a
Lagrangian point $X_{s}$.

Combining the two equations \ref{eq:corr3.4} et \ref{eq:corr3.6}, an
analytic expression for the volume source term for the energy
equation to be included in the IBM method is derived:
\begin{equation}
\label{eq:corr3.7}
{Q}_e^{n+1/2}=a_e \, ((\rho_L E)^{d} - (\rho_L E)^{l})
\end{equation}
Equation \ref{eq:corr3.7} can be further manipulated if one
considers that $e= c_{v} \theta$ for a perfect gas. Therefore, the
volume source term for the internal energy equation can be written
as a function of the temperature:
 \begin{equation}
\label{eq:corr3.8}
{Q}_e^{n+1/2}=a_e \, c_{v} ((\rho_L \Theta)^{d} - (\rho_L \Theta)^{l})
\end{equation}
with
\begin{itemize}
\item $c_{v}$ : the thermal capacity at constant volume of the interpolated fluid at Lagrangian markers
\item $(\rho_L \Theta)^{d}$ : the desired temperature density at the Lagrangian markers
\item $(\rho_L \Theta)^{l}$ : interpolated temperature density at the Lagrangian markers for
a neutral system.
\end{itemize}

A choice relying on numerical arguments has also been performed for the term $\mathbf{f} \cdot \mathbf{u}$ in the energy equation, which has been systematically set to zero. This term is naturally zero when the body is still, but it shows non-zero values in the case of body movement. However, this term is always very small so interpolation errors are most of the time larger in magnitude, leading to numerical instability and degradation of the accuracy. We therefore decided to remove it from the energy equation.

\subsubsection{IBM implementation in OpenFOAM numerical solvers}\label{sec3:impl}

The implementation of the IBM algorithm has been performed in the
open source library OpenFOAM, owing to previous development of the
research group \cite{Constant_2017,Riahi_2018}. Two solvers for
compressible flows have been considered: \textit{sonicFoam} and
\textit{rhoCentralFoam}. These two solvers are studied to provide
optimal performance for different ranges of $Ma$ numbers. The main
features of the solvers are discussed in \cite{Riahi_2018}. The
newly generated solvers will be referred to in the following as
\textit{IBM-HT} versions of the initial solver modified and are now
presented. One must consider that the versions of the solvers
proposed in \cite{Riahi_2018} are \textit{neutral} with respect to
the heat transfer behavior i.e. the presence of the body surface is
transparent with respect to heat exchange. Details about the algorithmic structure of \textit{IBM-HT-sonicFoam} and \textit{IBM-HT-rhoCentralFoam} are provided in \ref{apendixA} and \ref{apendixB}, respectively.

 {For every test case investigated with these solvers, the numerical discretization used for the gradient and divergence operators in the dynamic equations is linear (i.e. second order centered schemes). A van Leer limiter, which improves the flow prediction in the regions where shocks are observed, is used for the flux interpolation. Also, a grid convergence study is provided in \ref{sec:appC}}

\section{Numerical validation of the \textit{IBM-HT} algorithms} \label{sec::Validation}

Validation of the new solvers is performed on the 2D flow around a
heated circular cylinder. This classical test case has been
extensively investigated in the literature for a large spectrum of
values of $Re$, $Ma$ and temperature ratio
$r_\theta=\frac{\theta_{wall}}{\theta_{inlet}}$, where $\theta_{wall}$ is the wall temperature of the immersed body and $\theta_{inlet}$ is the temperature imposed at the inlet. Numerous
databases are available for comparison.

\subsection{Test case - numerical details}
\label{sec::validationCylinder}
The center of the cylinder is the origin of the domain in $(0,0)$ and its diameter is equal to $D$.
The dimensions of the computational domain are $[-16D,48D] \times
[-16D,16D]$  in the streamwise direction $(x)$ and normal direction
$(y)$, respectively. The mesh elements are uniformly distributed in
the central zone near to cylinder $x \times y \in [-D, \,D] \times
[-D, \, D]$. Their size is equal to $\Delta x = \Delta y = 0.01D$.
Outside this region, the mesh is progressively coarsened moving
through five different levels as shown in Figure
\ref{fig::fig_thermique}. The passage from one level to another
implies a coarsening ratio of the elements of $2$ in both $x$ and
$y$ directions. This value is smoothed in correspondence of the
boundaries of the refinement regions.  {This transition in the mesh is generated via an OpenFOAM native tool, \textit{snappyHexMesh}, and the Finite Volume discretization used in the solvers suitably takens it into account. In fact, the solver imposes conservation of the fluxes through the faces of the mesh elements, granting conservativity of the Navier--Stokes equation even with abrupt changes in the resolution of the mesh. Nonetheless, we have paid attention that these transition regions do not occur near the Eulerian supports for each Lagrangian Marker, for this test case and all the other flow configurations presented in the following Sections. This tasks has been easy to perform as it is connected with the need to impose a refined mesh distribution in the vicinity of the immersed body.} The total number of mesh
elements that result is equal to $1.5 \times 10^5$.  {The positioning of the Lagrangian markers is performed so that each Eulerian mesh element crossed by the body surface includes at least one marker. This strategy has been performed as the original version of the method has shown sensitivity to \textit{holes} in the connection between the Eulerian and the Lagrangian space \cite{Pinelli_JCP_2010}. The IBM library integrates an internal check routine which keeps only the Lagrangian marker closer to the center of each mesh element, if one or more markers are detected. For the mesh previously introduced, $401$ Lagrangian markers have been used.}

 {The boundary conditions have been chosen accordingly to general guidelines in the literature and using the team's experience for each of the $Ma$ number investigated, selecting different options for subsonic and supersonic flows. More precisely, a} constant
velocity inlet condition is imposed upstream (left side), a mass
conserving outlet condition is prescribed downstream and  {slip (subsonic) or non
reflective (supersonic)} conditions are imposed in the normal direction.

For each case analyzed, a number of different coefficients are
compared with available data of the literature. Aerodynamic forces
are determined via the bulk flow coefficients known as drag
coefficient $C_d$ and lift coefficient $C_l$:

\begin{equation}
\label{eq:Bulk}
C_d = \frac{2 F_x}{\rho_{\infty} U_{\infty}^2} \; , \qquad C_l = \frac{2 F_y}{\rho_{\infty} U_{\infty}^2}
\end{equation}

the forces in the streamwise direction $F_x$ and in the normal $F_y$
are directly calculated on the Lagrangian markers. $\rho_{\infty}$
and $U_{\infty}$ indicate asymptotic physical quantities imposed at
the inlet.

\begin{figure}[ht]
\centering
\includegraphics [width=0.9\linewidth, height=0.5\linewidth]{domaine_thermique.png}
\caption{\label{fig::fig_thermique} 2D computational domain used for \textit{IBM-HT} validation.}
\end{figure}

\subsection{Unsteady flow around a heated 2D circular cylinder}
The flow around a heated circular cylinder in unsteady regime is now
investigated to further assess the performance of the source term
for the internal energy equation. The non-dimensional parameters are
set to $Ma=0.01$ using temperature ratio $r_\theta = 1.8$ and three
Reynolds number $Re= \{100 , 140 , 260\}$. For these cases,
variations in the temperature field significantly influence the flow
organization, especially when the temperature ratio $r_\theta$ is
larger than 1.1 \cite{depalma_2006}.

For the different cases, an unsteady behavior characterized by a
periodic von K\'arm\'an wake is observed. Qualitative visualizations
of the instantaneous regimes are shown with the temperature
isocontours in Figure \ref{fig:instationnaire} $(a)$. In addition
the time evolution of the lift coefficient is shown in Figure
\ref{fig:instationnaire} $(b)$.

\begin{figure}
\begin{tabular}{cc}
\includegraphics[width=0.48\linewidth]{t_insta_re100_1_8.pdf} & \includegraphics[width=0.48\linewidth]{cl_re100_teta18.pdf} \\
(a) & (b)
\end{tabular}
\caption{\label{fig:instationnaire} Isocontours of the temperature
(a) and time evolution of the lift coefficient (b) for the flow past
a circular cylinder for $Ma=0.01$, $Re = 100$ and $r_\theta=1.8$.
$t_a = D / U_{\infty}$ is the characteristic advection time.}
\end{figure}

One important parameter describing the flow is the Strouhal number
$S_t=\frac{f D}{U_\infty}$, where $f$ is the shedding frequency
computed using the time evolution of the lift coefficient $C_l$.
Comparison shows a good agreement with results available in the
literature \cite{depalma_2006,deTullio_2007}, see Figure \ref{ta:st_2d}.
\begin{figure}
\center \includegraphics[width=0.8\linewidth]{re_st_depalma.png}
\caption{\label{ta:st_2d} Strouhal number versus Reynolds number for flow past a heated cylinder: comparison between present results and literature data}
\end{figure}
We observe that the vortex shedding frequency, $f$, and thus the Strouhal
number $S_t$, increases for increasing values of $Re$ for a given
$r_\theta$.
Finally, in table \ref{ta:table3} the computed bulk flow
quantities for the flow for all values of $Re$ are reported.  {It can
be seen that the average value of the drag coefficient $C_d$ decreases when $Re$ increases. On the other
hand, the maximum value of the coefficient of lift $C_l^{max}$ exhibits the opposite variation, increasing with $Re$. These results are in agreement with the findings by Homsi et al. \cite{Homsi_2021}.}
\begin{table}[h!]
\centering
\begin{tabular}{c|c|c}
Present study - $\mathbf{r}_\theta=1.8$ & $C_d$  &  $C_l^{max}$  \\ \hline
$Re=100$ & 1.48    & 0.18 \\
$Re=140$ & 1.46    & 0.29 \\
$Re=260$ & 1.45    & 0.53 \\
\end{tabular}
\caption{ $C_d$ and $C_l^{max}$ for unsteady flow around heated
circular cylinder at $Ma=0.01$. } \label{ta:table3}
\end{table}

\section{Forced convection from an heated sphere in a steady regime} \label{sec::Sphere}

The three-dimensional flow around a sphere is now investigated. This test case, which is significantly more complex, allows for a robust validation of the capability of the IBM method to accurately capture heat exchange.

\subsection{Numerical and computational ingredients}

The computational domain is here set to $x \times y \times z =
[-16D, \, 48D] \times [-16D, \, 16D] \times [-16D, \, 16D]$ where
$D$ is the diameter of the sphere. Again, the center of the body is
set in the origin of the system and the mesh is obtained using the tool \textit{snappyHexMesh}. The mesh is made by hexahedral
uniform elements which are progressively refined approaching the
sphere region (see figure \ref{fig::fig7}). The size of the elements
is refined by a factor two in each space direction when crossing the
prescribed interfaces between regions at different resolution. The
central most refined region is defined by the coordinates $x \times
y \times z = [-1.25D, \, 1.25D] \times [-1.25D, \, 1.25D] \times
[-1.25D, \, 1.25D]$. Within this region, the mesh resolution is
$\Delta x=\Delta y=\Delta z = 1/64 D$. This mesh is composed by a
total of $5 \times 10^6$ elements and was already used for analysis of compressible flows \cite{Riahi_2018}.
 {The choice of the number and positioning of the Lagrangian markers follows the same strategy presented for the flow around a circular cylinder discussed in sec. \ref{sec::validationCylinder}. Here $12119$ Lagrangian markers are used.}

\begin{figure}[ht]
\centering
\includegraphics [width=0.9\linewidth, height=0.5\linewidth]{mesh_3d.pdf}
\caption{\label{fig::fig7} Visualization of cutting planes inside the 3D mesh used for the calculation of the flow around a sphere.}
\end{figure}
In order to provide a suitable representation of the physical
features of the flow, the \textit{IBM-HT-rhoCentralFoam} solver is
employed for this study.

\subsection{Physical regimes observed for several $r_\theta$}

This test case has been chosen because of the emergence of different
regimes which exhibit a very high sensitivity to the values of
temperature ratio $r_\theta=1.1, \, 1.5, \, 2$, representing a challenging test case
of validation. We choose to investigate heated flow around sphere at
$Ma=0.8$ and $Re=300$ and the isocontours of the velocity field (represented in terms of $Ma$) and the temperature are shown in Figure
\ref{fig:3d_sphere} for the three cases. One can see that the flow undergoes a transition from an unsteady regime to a steady axisymmetric state. This effect is clearly driven by the heat exchange dynamics. Results are compared with
recent data reported in the literature for body-fitted numerical
simulations using high order discretization schemes \cite{Nagata_thermo} (Table \ref{ta:table6}). The physical quantities compared are the mean drag coefficient $\overline{C_d}$, the length of the recirculation bubble $l_s$ and the Strouhal number $St$. One can see that a reasonably good agreement is obtained, and in particular variations with different $r_{\theta}$ values are very similar. It appears that the IBM tends to provide a slight over prediction (around $10\%$) of the mean drag coefficient $\overline{C_d}$ when compared with body-fitted simulations. This result, which was already observed in previous studies \cite{Riahi_2018}, is probably associated with mesh resolution and the numerical schemes employed for the analysis.

\begin{table}[h!]
\centering
\begin{tabular}{c|c|c|c|c}
&Studies & $\overline{C_d}$  &  $l_s$ & $St$  \\ \hline
\textbf{Ma=0.8} &\textbf{Present ($r_\theta=1.1$)}  & 0.895 & 2.9 & 0.104  \\
\textbf{Re=300}&Nagata \cite{Nagata_thermo} (Num.)  & 0.81       & 2.8 & 0.105  \\\hline
&\textbf{Present ($r_\theta=1.5$)}  & 0.93       &3.18 &   \\
&Nagata \cite{Nagata_thermo} (Num.)  &    0.85    & 3.1 &   \\\hline
&\textbf{Present ($r_\theta=2$)}  & 0.99 & 3.4 &   \\
&Nagata \cite{Nagata_thermo} (Num.) &    0.91    & 3.35 &   \\\hline
\end{tabular}
\caption{Bulk flow quantities for the flow past a sphere obtained
via IBM simulation. Present results are compared with available data
in the literature. $\overline{C_d}$ is the mean, time-averaged drag coefficient,
$l_s$ is the recirculation length and $St$ is the Strouhal number.}
\label{ta:table6}
\end{table}

\begin{figure}
\begin{tabular}{cc}
\includegraphics[width=0.48\linewidth]{Ma08_tr1.pdf} & \includegraphics[width=0.48\linewidth]{T08_tr1.pdf} \\
(a) & (b) \\
\includegraphics[width=0.48\linewidth]{Ma08_tr1_5.pdf} & \includegraphics[width=0.48\linewidth]{T08_tr1_5.pdf} \\
(c) & (d) \\
\includegraphics[width=0.48\linewidth]{Ma08_tr2.pdf} & \includegraphics[width=0.48\linewidth]{T08_tr2.pdf} \\
(e) & (f)
\end{tabular}
\caption{\label{fig:3d_sphere} Iso-contours of $Ma$ number and
temperature for $Ma=0.8$ and $Re=300$: (a) $r_\theta=1.1$;(b)
$r_\theta=1.5$;(c) $r_\theta=2$}
\end{figure}

\subsection{Investigation of the supersonic flow around a sphere}
The supersonic flow for $Ma=2$, $Re=300$ and $r_{\theta}=2$ is now investigated, using the same mesh previously employed. In this case compressibility effects are
very strong and a steady axisymmetric configuration is observed, as qualitatively shown in figure \ref{fig:superHighRe}. The main bulk flow quantities are again compared with results obtained from body-fitted numerical simulations \cite{Nagata_thermo} and reported in table \ref{ta:table7}. It appears that all the physical
features are in agreement with the data in the literature \cite{Nagata_thermo}, assessing the potential of the IBM method developed. A slight over prediction  {(this time around $5 \%$) } of the drag coefficient $C_D$ is again observed, while the recirculation length $l_s$ and the distance of the shock from the wall $\Delta_{shock}$ are very similar.

\begin{table}[h!]
\centering
\begin{tabular}{c|c|c|c|c}
&Studies & $\overline{C_d}$  &  $l_s$ & $\Delta_{shock}$ \\ \hline
\textbf{Ma=2}&\textbf{Present ($r_\theta=2$)}  & 1.48       & 0.85    & 0.23\\
\textbf{Re=300}&Nagata \cite{Nagata_thermo} (Num.) & 1.41       & 0.85  & 0.25 \\\hline
\end{tabular}
\caption{Bulk flow quantities for the flow past a sphere obtained
via IBM simulation. Present results are compared with available data
in the literature. $C_d$ is the time-averaged drag coefficient,
$x_s$ is the recirculation length, $Nu$ Nusselt number and
$D_{shock}$is  the shock distance.} \label{ta:table7}
\end{table}

\begin{figure}
\begin{tabular}{cc}
\includegraphics[width=0.48\linewidth]{Ma_superso.pdf} & \includegraphics[width=0.48\linewidth]{T_superso.pdf} \\
(a) & (b)
\end{tabular}
\caption{\label{fig:superHighRe} Numerical simulation of the
supersonic flow around a sphere for $Ma=2$, $Re=300$ and
$r_\theta=2$ : (a) Isocontours of the $Ma$ number  (b) isocontours
of the Temperature}
\end{figure}

\section{Fluid - Structure interaction : oscillating cylinder with fixed wall temperature}\label{sec:FSI}

In this last section, the flow around an oscillating two-dimensional
cylinder is investigated. The simulation is performed using the $
IBM-HT-rhoCentralFoam $ solver.

The motion of the cylinder is imposed in the streamwise direction $x$, providing an explicit law of movement for the Lagrangian markers following a
sinus law. More precisely, the position in the streamwise direction of the Lagrangian
markers is updated in time following the equation:
\begin{equation}
X_s (t) = X_s (0) + A sin (2 \pi f_c t)
\end{equation}

with $A=0.14 \times D$ and $f_c=2 \times f_0$. A preliminary
calculation is performed on a fixed cylinder to determine the
shedding frequency $f_0$. The values of the parameters have been
chosen in order to allow comparison with a study proposed by Luo et
al. \cite{Luo_2016} with  {$Ma=0.3$}, $Re=100$ and $r_\theta=1$. In figure
\ref{fig:mvt_thermique} the qualitative evolution of the isocontours
of the temperature and of the vorticity are shown after the initial
transient. The comparison of present results with the findings by
Luo et al. \cite{Luo_2016} it terms of the time evolution
of  { the drag coefficient $C_d$} and lift coefficient $C_l$ is very good, as shown in figure \ref{fig::cd_cl_osc_therm_com}. Therefore, the \textit{IBM-HT} model has proven good accuracy even for a more complex case including fluid-structure interaction.

\begin{figure}
\begin{tabular}{cc}
\includegraphics[width=0.48\linewidth]{vorticity_t_1_4.pdf} & \includegraphics[width=0.48\linewidth]{temperature_t_1_4.pdf} \\
(a) & (b)
\end{tabular}
\caption{\label{fig:mvt_thermique} Iso-contours vorticity (a) and
temperature (b) for the flow around a longitudinal oscillating
cylinder $Re=100$, $r_\theta=1$ }
\end{figure}

\begin{figure}
\begin{tabular}{cc}
\includegraphics[width=0.48\linewidth]{cd_thermique_new_comp.png} & \includegraphics[width=0.48\linewidth]{cl_thermique_new_comp.png} \\
(a) & (b)
\end{tabular}
\caption{\label{fig::cd_cl_osc_therm_com}  {Comparison of the (a)
drag coefficient $C_d$ and (b) lift coefficient $C_l$ for the flow around a longitudinal oscillating cylinder, $r_\theta=1$ with results by Luo et al. \cite{Luo_2016}}}
\end{figure}

 {An additional simulation with imposed temperature $r_{\theta}=2$ was carried out in order to study the effects of the increase in temperature on the aerodynamic parameters \ref{fig:mvt_thermique_t2}. The variation in temperature affects the viscosity close to the wall, which modifies the aerodynamic performance.}  {The comparison of the drag coefficient $C_d$ and lift coefficient $C_l$ for the two cases, which are shown in fig. \ref{fig::cd_cl_osc_therm_t_2}, shows important differences. In particular, the maximum value of $C_l$ decreases for the case with $r_\theta=2$, in agreement with the findings by Homsi et al. \cite{Homsi_2021}.}

\begin{figure}
\begin{tabular}{cc}
\includegraphics[width=0.48\linewidth]{vorticity_t_2.pdf} & \includegraphics[width=0.48\linewidth]{temperature_t_2.pdf} \\
(a) & (b)
\end{tabular}
\caption{\label{fig:mvt_thermique_t2}  {Iso-contours vorticity (a) and
temperature (b) for a flow around a longitudinal oscillating
cylinder with $r_{\theta}=2$ }}
\end{figure}

\begin{figure}
\begin{tabular}{cc}
\includegraphics[width=0.48\linewidth]{cd_thermique_t2_new_comp.png} & \includegraphics[width=0.48\linewidth]{cl_thermique_t2_new_comp.png} \\
(a) & (b)
\end{tabular}
\caption{\label{fig::cd_cl_osc_therm_t_2} { (a)  Temporal evolution of the
drag coefficient $C_d$ and (b) of the
lift coefficient $C_l$ for flow around an oscillating circular
cylinder. The cases for $r_{\theta}=1$ and $r_{\theta}=2$ are compared}}
\end{figure}
\section{Conclusion}\label{sec:conclusion}
In this paper, a discrete-forcing immersed boundary method was
proposed for simulating both heat exchange and compressible flows, for
stationary and moving bodies. The proposed methodology introduces an
additional heat source term to the right hand side of the
energy equation to keep the solid boundary at the prescribed
temperature. The compressible Navier-Stokes equations are
discretized on a Cartesian grid and solved
by either the sonicFoam or the rhoCentralFoam solver of the OpenFOAM
platform according to the Mach number of the flow.

To validate the robustness and accuracy of the proposed method, a
series of simulations was  performed on canonical forced convection
flows with comparison with existing results. Both the flow dynamics
and the heat transfer were in good agreement with the data in
literature. The 3D flows around an isothermal sphere in the
transonic and supersonic regimes were investigated. Different wall
temperature conditions were applied leading to the transition
between an unsteady periodical wake and a steady axisymmetric state.
The results obtained with the proposed method were found in close
agreement with recent direct numerical simulations available in the
literature, demonstrating the capability of the numerical strategy
to treat both compressible and heated flows with satisfactory
accuracy. Finally the method was successfully validated on a
benchmark case involving moving boundaries, for which the
computational cost of remeshing techniques would be significantly high for classical body-fitted approaches.
Future works envision an improvement in accuracy for higher Reynolds number configurations, including wall modelling within the formulation of the IBM model.

\section*{Acknowledgement}

{The research work has been developed using computational resources within the framework of the project gen7590-A0012A07590 DARI-GENCI.}

\bibliographystyle{unsrt}