Variations of largest rectangle recognition amidst a bichromatic point set

TL;DR

This paper introduces in-place algorithms for maximum monochromatic rectangle recognition in bichromatic point sets, along with an in-place k-d tree construction supporting efficient range queries, and extends to maximum weight rectangle computation.

Contribution

It presents novel in-place algorithms for maximum monochromatic rectangles and an in-place k-d tree supporting orthogonal range queries with optimal space and time complexities.

Findings

In-place algorithms for maximum monochromatic rectangles in 2D and 3D.

An in-place k-d tree supporting range reporting and counting queries.

Efficient maximum weight rectangle algorithm for weighted point sets.

Abstract

Classical separability problem involving multi-color point sets is an important area of study in computational geometry. In this paper, we study different separability problems for bichromatic point set P=P_r\cup P_b on a plane, where $P_r$ and $P_b$ represent the set of n red points and m blue points respectively, and the objective is to compute a monochromatic object of the desired type and of maximum size. We propose in-place algorithms for computing (i) an arbitrarily oriented monochromatic rectangle of maximum size in R^2, (ii) an axis-parallel monochromatic cuboid of maximum size in R^3. The time complexities of the algorithms for problems (i) and (ii) are O(m(m+n)(m\sqrt{n}+m\log m+n \log n)) and O(m^3\sqrt{n}+m^2n\log n), respectively. As a prerequisite, we propose an in-place construction of the classic data structure the k-d tree, which was originally invented by J. L. Bentley in 1975. Our in-place variant of the $k$-d tree for a set of n points in R^k supports both orthogonal range reporting and counting query using O(1) extra workspace, and these query time complexities are the same as the classical complexities, i.e., O(n^{1-1/k}+\mu) and O(n^{1-1/k}), respectively, where \mu is the output size of the reporting query. The construction time of this data structure is O(n\log n). Both the construction and query algorithms are non-recursive in nature that do not need O(\log n) size recursion stack compared to the previously known construction algorithm for in-place k-d tree and query in it. We believe that this result is of independent interest. We also propose an algorithm for the problem of computing an arbitrarily oriented rectangle of maximum weight among a point set P=P_r \cup P_b, where each point in P_b (resp. P_r) is associated with a negative (resp. positive) real-valued weight that runs in O(m^2(n+m)\log(n+m)) time using O(n) extra space.

Full text

Keywords: Bichromatic point set; obstacle-free rectangle recognition; orthogonal range counting; in-place k-d tree; maximum weight rectangle recognition; space efficient algorithms.

1 Introduction

Given a bichromatic point set $P=P_{r}\cup P_{b}$, where $P_{r}$ is the set of $n$ red points and $P_{b}$ is the set of $m$ blue points, the basic separability problem is to find a separator $S$ such that the points in $P_{r}$ and $P_{b}$ lie in two different sides of $S$ respectively. The motivation for studying this separability problem for a bichromatic point set stems from its various applications in facility location, VLSI layout design, image analysis, data mining, computer graphics and other classification based real life scenarios [19, 22, 23, 25, 26]. The bichromatic separability problem also has its application to detect obstacle free separators. In the literature, different types of separators like hyperplane [30], circle [36], rectangle [18, 25, 39], square [13, 38] has been studied to optimize the objective function of the corresponding problem. In this paper, we focus on designing space-efficient algorithms for the following problems:

• P1

Computing an arbitrarily oriented monochromatic rectangle of maximum size ($LMR$) in $\mathbb{R}^{2}$, where a rectangle $U$ is said to be monochromatic if it contains points of only one color in the proper interior of $U$, and the size of $U$ is the number of points of that color inside or on the boundary of the rectangle $U$.

• P2

Computing an arbitrarily oriented rectangle of maximum weight $(LWR)$ in $\mathbb{R}^{2}$, where each point in the set $P_{b}$ (resp $P_{r}$) is associated with negative (resp. positive) real-valued weight, and the weight of a rectangle $U$ is the sum of weights of all the points inside or on the boundary of $U$.

• P3

Computing a monochromatic axis parallel cuboid222a solid which has six rectangular faces at right angles to each other ($LMC$) of maximum size in $\mathbb{R}^{3}$.

A rectangle of arbitrary orientation in $\mathbb{R}^{2}$ is called red if it does not contain any blue point in its interior333blue points may appear on the boundary. The largest red rectangle ($LRR$) is a red rectangle of maximum size. Similarly, the largest blue rectangle ($LBR$) is defined. The largest monochromatic rectangle ($LMR$) is either $LRR$ or $LBR$ depending on which one is of maximum size. Here, the objective is to compute the $LRR$. In $\mathbb{R}^{3}$, we similarly define the largest axis parallel red cuboid ($LRC$), i.e. a cuboid containing the maximum number of red points and no blue point in its interior1. We use $x(p)$ and $y(p)$ to denote the $x$- and $y$-coordinate of a point $p\in P$ respectively.

Several variations of this problem are well studied in the literature. In the well-known maximum empty rectangle (MER) problem, a set $P$ of $n$ points is given; the goal is to find a rectangle (axis parallel/arbitrary orientation) of maximum area that does not contain any point of $P$ in its interior (see [1, 17, 33, 35] for MER of fixed orientation, and [16, 32] for MER of arbitrary orientation). For fixed orientation version, the best-known algorithm runs in $O(n\log^{2}n)$ time and $O(n)$ space [1]. For arbitrary orientations version, the best-known algorithm runs in $O(n^{3})$ time using $O(n)$ space [16].

For the bichromatic version of the problem, Liu and Nediak [29] designed an algorithm for finding an axis parallel $LRR$ of maximum size in $O(n^{2}\log n+nm+m\log m)$ time using $O(n)$ space. Backer and Keil [3] improved the time complexity to $O((n+m)\log^{3}(n+m))$ using $O(n\log n)$ space adopting the divide-and-conquer approach of Aggarwal and Suri [1]. They also proposed an $O((n+m)\log(m+n))$ time algorithm for finding an axis-parallel red square of maximum size. Recently, Bandyapadhyay and Banik [5] proposed an algorithm for finding the $LRR$ in arbitrary orientation using $O(g(n,m)\log(n+m)+n^{2})$ time and $O(n^{2}+m^{2})$ space, where $g(n,m)\in O(m^{2}(n+m))$ and $g(n,m)\in\Omega(m(n+m))$.

Other variations of the $LRR$ problem, studied in the literature are as follows. For a given bichromatic (red,blue) point set, Armaselu and Daescu [2] considered the problem of finding a rectangle of maximum area containing all red points and minimum number of blue points. In $\mathbb{R}^{2}$, the axis-parallel version of this problem can be solved in $O(m\log m+n)$ time and the arbitrary oriented version requires $O(m^{3}+n\log n)$ time. In $\mathbb{R}^{3}$, the axis-aligned version of the problem can be solved in $O(m^{2}(m+n))$ time. Eckstein et al. [25] considered the axis-parallel version of the $LRR$ problem in higher ($d\geq 3$) dimensions. They showed that, if the dimension $d$ is not fixed, the problem is $NP$-hard. They presented an $O(n^{2d+1})$ time algorithm for any fixed dimension $d\geq 3$. Later, Backer and Keil [4] improved the time bound of the problem to $O(n^{d}\log^{d-2}n)$. Cortés et al. [18] considered the problem of removing as few points as possible from the given bichromatic point set such that the remaining points can be enclosed by two axis-parallel rectangles $A_{R}$ and $A_{B}$ (may or may not be disjoint), where $A_{R}$ (resp. $A_{B}$) contains all the remaining red (resp. blue) points. They solved this problem in $O(n^{2}\log n)$ time using $O(n)$ space. The problem of separating bichromatic point sets by two disjoint axis-parallel rectangles such that each of the rectangles is monochromatic, is solved in $O(n\log n)$ time by Moslehi and Bagheri [31] (if such a solution exists). If these two rectangles are of arbitrary orientation then they solved the problem in $O(n^{2}\log n)$ time. Bitner et al. [10] studied the problem of computing the minimum separating circle, which is the smallest circle containing all the points of red color and as few points as possible of blue color in its interior. The proposed algorithm runs in $O(nm\log m+n\log n)$ time using $O(n+m)$ space. They also presented an algorithm for finding the largest separating circle in $O(nm\log m+k(n+m)\log(n+m))$ time using $O(n+m)$ space, where $k$ is the number of separating circles containing the smallest possible number of points from blue point set. The problem of covering a bichromatic point set with two disjoint monochromatic disks has been studied by Cabello et al. [14], where the goal is to enclose as much points as possible in each of the monochromatic disks. They solved the problem in $O(n^{\frac{11}{3}}\polylog~{}n)$ time. If the covering objects are unit disks or unit squares, then it can be solved in $O(n^{\frac{8}{3}}\log^{2}n)$ and $O(n\log n)$ time respectively [13].

The weighted bichromatic problems are also studied in the literature. The smallest maximum-weight circle for weighted points in the plane has been addressed by Bereg et al. [8]. For $m$ negative weight points and $n$ positive weight points they solved the problem in $O(n(n+m)\log(n+m))$ time using linear space. For a weighted point set Barbay et al. [6] provided an $O(n^{2})$ time algorithms to find the maximum weight axis-parallel square.

Our Contribution

Given a bichromatic point set $P=P_{r}\cup P_{b}$ in a rectangular region $\cal A\subseteq$ $\mathbb{R}^{2}$, where $P_{r}$ and $P_{b}$ are set of $n$ red points and $m$ blue points respectively, we design an in-place algorithm444An in-place algorithm is an algorithm where the input is given in an array, the execution of the algorithm is performed using only $O(1)$ extra workspace, and after the execution of the algorithm all the input elements are present in the array. for finding an arbitrarily oriented $LMR$ of maximum size in $O(m(n+m)(m\sqrt{n}+m\log m+n\log n))$ time, using $O(1)$ extra workspace. We also show that the axis-parallel version of the $LMR$ problem in $\mathbb{R}^{3}$ (called the $LMC$ problem) can be solved in an in-place manner in $O(m^{3}\sqrt{n}+m^{2}n\log n)$ time using $O(1)$ extra workspace. As a prerequisite of the above problems, we propose an algorithm for constructing a k-d tree with a set of $n$ points in $\mathbb{R}^{k}$ given in an array of size $n$ in an in-place manner such that the orthogonal range counting query can be performed using $O(1)$ extra workspace. The construction and query time of this data structure is $O(n\log n)$ and $O(n^{1-1/k})$, respectively. Finally, we show that if the points in $P_{r}$ (resp. $P_{b}$) have positive (resp. negative) real-valued weight, then a rectangle of arbitrary orientation with maximum weight (called $LWR$) can be computed in $O(m^{2}(n+m)\log(n+m))$ time using $O(n)$ space.

2 In-place k-d tree

To perform the orthogonal range reporting query, Bentley [7] invented k-d tree in 1975. It is a binary tree in which every node is a k-dimensional point. Every non-leaf node can be thought of being associated with one of the k-dimensions of the corresponding point, with a hyperplane perpendicular to that dimension’s axis, and implicitly this hyperplane splits the space into two half-spaces. Points to the negative side of this splitting hyperplane are represented by the left subtree of that node and points in the positive side of the hyperplane are represented by the right subtree555For a hyperplane $x_{i}=c$, its negative (resp. positive) side is the half-space $x_{i}<c$ (resp. $x_{i}>c$), where $x_{i}$ is the $i$-th coordinate of a $k$ dimensional point $(x_{1},x_{2},\ldots,x_{k})$. Depending on the level of a node going down the tree, the splitting dimension is chosen one after another in a cyclic manner. Each node $v$ of the tree is associated implicitly with a rectangular region of space, called cell($v$). The cell corresponding to the root of the tree is the entire $\mathbb{R}^{k}$. A child’s cell is contained within its parent’s cell, and it is determined by the splitting hyperplane stored at the predecessor nodes.

Brönnimann et al. [12] mentioned an in-place version of the k-d tree. We note that their approach for both constructing the data structure and querying in the data structure are recursive, and need to remember the subarray and the cell in which the recursive call is invoked. As a result, there is a hidden $O(\log n)$ space requirement for system stack. We present an alternate variant of in-place $k$-d tree data structure that supports both reporting and counting query for orthogonal query range with same query time as the classical one. The advantage of this in-place variant is that both construction and query algorithms are non-recursive, and it takes only $O(1)$ extra workspace during the execution apart from the array containing the input points. The in-place organization of this data structure is similar to the in-place min-max priority search tree proposed by De et al. [20].

2.1 Construction of in-place k-d tree

Let us consider that a set $P$ of $n$ points in $\mathbb{R}^{k}$ is given in an array $P[1,\ldots,n]$. We propose an in-place algorithm to construct the k-d tree $\cal T$ in the array $P$. Here, $\cal T$ is a binary tree of height $h=\lfloor\log n\rfloor$, such that the levels $0,1,\ldots,h-1$ are full and level $h$ consists of $\varkappa=n-(2^{h}-1)$ nodes which are aligned as far as possible to the left. At the end of the construction, the tree $\cal T$ is stored implicitly in the given array $P$. In other words, we store the root of the tree in $P[1]$, its left and right children in $P[2]$ and $P[3]$, etc. This allows us to navigate $\mathord{\emph{parent}}(P[i])$ which is at $P[\lfloor\frac{i}{2}\rfloor]$, and $\mathord{\emph{left-child}}(P[i])$ and $\mathord{\emph{ right-child}}(P[i])$, if they exist, which are at $P[2i]$ and $P[2i+1]$, respectively.

Note that there are $2^{i}$ nodes in the level $i\neq h$ of the tree $\cal T$. As the number of leaf nodes in a full tree of height $h-i$ is $2^{h-i}$, so there are $k_{i}=\lfloor{\frac{n-\varkappa}{2^{h-i}}}\rfloor$ nodes at level $i$ ($0<i<h$) that are roots of subtrees, each of size $K_{1}^{i}=2^{h+1-i}-1$. If $k_{i}=2^{i}$, then all the subtrees at level $i$ are full, and the number of nodes in each subtree is $2^{h+1-i}-1$. Otherwise, we have $k_{i}<2^{i}$, and level $i$ ($0<i<h$) consists of, from left to right,

• •

$k_{i}$ nodes which are roots of subtrees, each of size $K_{1}^{i}=2^{h+1-i}-1$,

• •

one node that is the root of a subtree of size $K_{2}^{i}=2^{h-i}-1+\varkappa-k_{i}\cdot 2^{h-i}$, and

• •

$2^{i}-k_{i}-1$ nodes which are roots of subtrees, each of size $K_{3}^{i}=2^{h-i}-1$.

See Figure 1 for an illustration.

Here, we introduce the notion of block and block median. Assume that $0<i<h$. We refer to the portion of the array $P[(2^{i}+(j-1)K_{1}^{i}),\ldots,(2^{i}+jK_{1}^{i}-1)]$ as block $B_{j}^{i}$, for $j\leq k_{i}$. The portion of the array $P[(2^{i}+k_{i}K_{1}^{i}),\ldots,(2^{i}+k_{i}K_{1}^{i}+K_{2}^{i}-1)]$ is referred to as block $B_{k_{i}+1}^{i}$, and $P[(2^{i}+k_{i}K_{1}^{i}+K_{2}^{i}+(j-1)K_{3}^{i}),\ldots(2^{i}+k_{i}K_{1}^{i}+K_{2}^{i}+jK_{3}^{i}-1)]$ are referred to as blocks $B_{j}^{i}$, for all $j>k_{i}+1$. For $i=0$, we refer to the whole array $P[1,\ldots,n]$ as $B_{1}^{0}$. For $i=h$, we refer to the array element $P[2^{h}+j]$ as block $B_{j}^{h}$, where $1\leq j\leq\varkappa$. For a block $B_{j}^{i}$ ($0<i<h$) of size $K_{1}^{i}$ (resp. $K_{3}^{i}$), we denote block median $m_{j}^{i}$ as a point in $B_{j}^{i}$ whose $(i\mod k)$-th coordinate value is $\lceil\frac{K_{1}^{i}}{2}\rceil$-th (resp. $\lceil\frac{K_{3}^{i}}{2}\rceil$-th) smallest among all the points in $B_{j}^{i}$. If the size of $B_{j}^{i}$ is $K_{2}^{i}$, then depending on whether $K_{2}^{i}-(2^{h-i}-1)<2^{h-i-1}-1$ or $K_{2}^{i}-(2^{h-i}-1)\geq 2^{h-i-1}-1$, we refer the block median $m_{j}^{i}$ as a point in $B_{j}^{i}$ whose $(i\mod k)$-th coordinate value is $K_{2}^{i}-(2^{h-i-1}-1)$-th or $2^{h-i}$-th smallest among all the points in $B_{j}^{i}$. For block $B_{1}^{0}$, depending on whether $n-(2^{h}-1)<2^{h-1}-1$ or $n-(2^{h}-1)\geq 2^{h-1}-1$, we refer the block median $m_{1}^{0}$ as a point in $B_{1}^{0}$ whose $1$st coordinate value is $n-(2^{h-1}-1)$-th or $2^{h}$-th smallest among all the points in $B_{1}^{0}$.

Our algorithm constructs the tree level by level. After constructing the $(i-1)$-th level of the tree, it maintains the following invariants:

Invariant 1**.**

• (i)

The subarray $P[1,\ldots,2^{i}-1]$ stores levels $0,1,\ldots,i-1$ of the tree.

• (ii)

Block $B_{j}^{i}$ contains all the elements of the $j$-th leftmost subtree of level $i$, for $j\in\{1,\ldots,2^{i}\}$ $(j\in\{1,\ldots,\varkappa\}$ when $i=h)$.

At the first iteration of the algorithm, we find the block median $m_{1}^{0}$ using the linear time in-place median finding algorithm of Carlsson and Sundström [15], and swap it with $P[1]$. Next, we arrange the subarray $P[2,\ldots,n]$ such that all the elements whose first coordinate value is greater than $m_{1}^{0}$ appear before all the elements whose first coordinate value is less than $m_{1}^{0}$. We can do this arrangement in linear time using $O(1)$ extra space.

Note that after the first iteration, both the invariants are maintained.

Assuming that the tree is constructed up to level $(i-1)$, now, we construct the tree up to level $i$ by doing the following:

Find block median $m_{j}^{i}$ from each block $B_{j}^{i}$ and swap it with the first location of block $B_{j}^{i}$. Using the median finding algorithm of [15], this needs a total of $O(n)$ time for all the blocks in this ($i$-th) level. 2. 2.

Now depending on the median value $m_{j}^{i}$ we arrange the elements of each block $B_{j}^{i}$ such that all the elements having $(i\mod k)$-th coordinate value greater than $m_{j}^{i}$ appears before all the elements having $(i\mod k)$-th coordinate value less than $m_{j}^{i}$. Thus the block $B^{i}_{j}$ splits into two parts, named first half-block and second half-block. This step again needs time proportional to the size of each block, and hence $O(n)$ time in total. 3. 3.

Now, we need to move all $m_{j}^{i}$ stored at the first position of each block to the correct position of level $i$ of the tree. To do this we do the following. First, we move the last block median $m_{2^{i}}^{i}$ next to $m_{2^{i}-1}^{i}$ by two swaps; (i) swap $m_{2^{i}}^{i}$ with the first element of the second half-block of $B_{2^{i}-1}^{i}$, and (ii) swap $m_{2^{i}}^{i}$ with the first element of the first half-block of $B_{2^{i}-1}^{i}$. Thus, after this swapping step all the elements in the block $B_{2^{i}-1}^{i}$ that are less than $m_{2^{i}-1}^{i}$ will stay before the elements greater than $m_{2^{i}-1}^{i}$. Now, we will move both the pair ($m_{2^{i}-1}^{i}$, $m_{2^{i}}^{i}$) just after $m_{2^{i}-2}^{i}$. It can be shown that, for the move of each element of this pair, we need a pair of swaps as explained above. Next, we will move $m_{2^{i}-2}^{i}$, $m_{2^{i}-1}^{i}$ and $m_{2^{i}}^{i}$ by swapping (as mentioned above) next to $m_{2^{i}-3}^{i}$. In this way, we will continue until all the block medians $\{m_{j}^{i}|j\in 2^{i}\}$ will become consecutively placed. Using $O(1)$ space, this can be done in linear time666The reason is that, during this step of execution each element is moved backward from its present position in the array at most once..

Step 3 ensures that both the invariants are maintained after this iteration.

At the end of $h$-th iteration, we have the tree $\cal T$ stored implicitly in the array $P$. The correctness of this algorithm follows by observing that the invariants are correctly maintained. As there are $O(\log n)$ iterations and each iteration takes $O(n)$ time, in total the algorithm takes $O(n\log n)$ time.

Lemma 1**.**

Given a set of $n$ points in $\mathbb{R}^{k}$ in an array $P$, the in-place construction of $k\text{-}d\text{-}tree$ takes $O(n\log n)$ time and $O(1)$ extra workspace.

2.2 Orthogonal range counting query in the in-place k-d tree

For the simplicity of explanation, we illustrate the range counting query for points in $\mathbb{R}^{2}$. We can easily generalize it for points in $\mathbb{R}^{k}$, for any fixed $k$. Given a rectangular range $Q=[\alpha,\beta]\times[\gamma,\delta]$ as a query, here, the objective is to return a count of the number of elements in $P$ that lie in the rectangular range $Q$.

In the traditional model, to answer counting query in $O(\sqrt{n})$ time each node in pre-proceesed k-d tree stores the subtree size. For our case, we cannot afford to store the subtree size along with each node of the in-place k-d tree. However, if we have the information of the level $\ell$ of a node $P[t]$, then we can on-the-fly compute the subtree size as follows. Note that $P[t]$ is $r=t-(2^{\ell}-1)$-th left most node at $\ell$-th level of the tree $\cal T$. Depending on whether $r\leq k_{\ell}$, $r=k_{\ell}+1$ or $r\geq k_{\ell}+2$, the subtree size of the node corresponding to $P[t]$ is $K_{1}^{\ell}$, $K_{2}^{\ell}$ or $K_{3}^{\ell}$. We want to remind the reader that $k_{\ell}=\lfloor{\frac{n-\varkappa}{2^{h-\ell}}}\rfloor$, $K_{1}^{\ell}=2^{h+1-\ell}-1$, $K_{2}^{\ell}=2^{h-\ell}-1+\varkappa-k_{\ell}\cdot 2^{h-\ell}$ and $K_{3}^{\ell}=2^{h-\ell}-1$, where $\varkappa=n-(2^{h}-1)$.

On the other hand, the traditional query algorithm is a recursive algorithm that starts from the root of the tree. At a node $v$, (i) if $Q\cap cell(v)=\emptyset$, then it returns 0; (ii) else if $cell(v)\subseteq Q$, then it returns the subtree size of $v$; (iii) otherwise, it recursively counts in the two children of $v$ and returns by adding these counts, accordingly. The main issue in implementing this algorithm in the in-place model is that it needs $O(\log n)$ space for system stack to have the knowledge of the corresponding $cell$ of a node. To tackle this situation, we have a new geometric observation which leads to a non-recursive algorithm in the in-place model.

At a node $v$, we can test whether the cells corresponding to both the children are intersecting the query region $Q$ or not, by checking whether the splitting plane stored at $\mathord{\emph{parent}}(v)$ is intersecting the query region $Q$ or not. If the splitting plane does not intersect, then the one of the child’s cell that has non-empty intersection with $Q$, can be decided by checking in which side of the hyperplane the region $Q$ lies. This simple trick works because when we are at a node $v$, we know that the cell corresponding to its parent has non-empty intersection with $Q$. The following observation plays a crucial role here.

Observation 1**.**

If the left (resp. right, bottom, and top) boundary of $cell(v)$ intersects the query region $Q$, then the left (resp. right, bottom, and top) boundary of $cell(v^{\prime})$ corresponding to the left (resp. right, left, right) child ($v^{\prime}$) of node $v$ also intersects the query region $Q$.

To decide whether $cell(v)\subseteq Q$, we do the following. Throughout the query algorithm, we keep a four-tuple $(L,R,B,U)$ each being able to store one coordinate value of the given input points. Initially, all of them are set to $NULL$. Throughout the query algorithm, this four-tuple maintains the following invariant:

Invariant 2**.**

When we are at a node $\mathord{\it Current}$, the non-NULL or NULL value stored at $L$ (resp. $R$, $B$, and $U$) implies that the left (resp. right, bottom, and top) boundary of the $cell(\mathord{\it Current})$ is intersecting or not intersecting the query region $Q$. More specifically, if the value stored at $L$ (resp. $R$, $B$, and $U$) is non-NULL777split-value of some node of the ancestor of $\mathord{\it Current}$, then it represents the left (resp. right, bottom, and top) boundary of the cell corresponding to the lowest level ancestor $v$ of the node $\mathord{\it Current}$, such that left (resp. right, bottom, and top) boundary of $cell(v)$ intersects the query region $Q$.

At a node $v$, if all the entries in the four-tuple is non-NULL, then the $cell(v)\subseteq Q$. We present our algorithm as a pseudocode in Algorithm 1. This is similar to the algorithm Explore in [20]. It uses two variables $\mathord{\it Current}$ and $\mathord{\it state}$ that satisfies the following:

• •

$\mathord{\it Current}$ is a node in $\cal T$.

• •

$\mathord{\it state}\in\{0,1,2\}$.

• •

If $\mathord{\it state}=0$, then either $cell(\mathord{\it Current})\subseteq Q$ or both the children of $\mathord{\it Current}$ need to be processed to compute $cell(\mathord{\it Current})\cap Q$.

• •

If $\mathord{\it state}=1$, then all elements of the set $Q\cap\left(\{\mathord{\it Current}\}\cup{\cal T}_{{\mathord{\emph{left-child}}(\mathord{\it Current})}}\right)$ have been counted, where ${\cal T}_{{\mathord{\emph{left-child}}(\mathord{\it Current})}}$ is the left subtree of $\mathord{\it Current}$ in the tree ${\cal T}$.

• •

If $\mathord{\it state}=2$, then all elements of the set $Q\cap{\cal T}_{\mathord{\it Current}}$ have been counted.

Update of the four-tuple $(L,R,T,B)$ is done as follows. While searching with the query rectangle $Q$ and with $\mathord{\it state}=0$, when $Q$ is split by the split-line of the node and the search proceeds towards one subtree of that node, we store the split-value (corresponding to the split-line) in the corresponding variable of the four-tuple provided it is not set earlier (contains $NULL$ value). During the backtracking, i.e, when $\mathord{\it state}=2$, if the split-value of the current node matches with the corresponding variable in the four-tuple, then the corresponding entity of the four-tuple is set to $NULL$. Now, if backtracking reaches from left, we set $\mathord{\it state}=1$. Since the right child of the current node needs to be processed, we set the corresponding entity of four-tuple with the split-value stored at that node.

The correctness of the algorithm follows from maintaining the invariants and Observation 1. In the worst case, we might have visited all the nodes whose corresponding cells overlap on the orthogonal query rectangle $Q$. As the number of cells stabbed by $Q$ can be shown to be $O(\sqrt{n})$ [9], we have the following result.

Lemma 2**.**

Given the in-place $2$-d tree maintained in the array $P$ of size $n$, the rectangular range counting query can be performed in $O(\sqrt{n})$ time using $O(1)$ extra workspace.

We can generalize, the above algorithm for points in $\mathbb{R}^{k}$. The only difference is that we need $2k$-tuple instead of four-tuple. Assuming $k$ is a fixed constant, we have the following:

Lemma 3**.**

Given the in-place k-d tree maintained in the array $P$ of size $n$, the orthogonal range counting query can be performed in $O(n^{1-1/k})$ time using $O(1)$ extra workspace.

3 $LMR$ problem in arbitrary orientation

In this section, we describe the method of identifying an arbitrarily oriented red rectangle of largest size for a given bichromatic point set $P=P_{r}\cup P_{b}$ in $\mathbb{R}^{2}$. The $LRR$ problem was solved by Bandyapadhyay and Banik [5], considering the blue points as obstacles, using the following observation:

Observation 2**.**

[5]** At least one side of a $LRR$ must contain two points $p,q$ such that $p\in P_{b}$ and $q\in P_{r}\cup P_{b}$, and other three sides either contain at least one point of $P_{b}$, or is open (unbounded) (see Figure 2).

For the sake of formulation of our problem, let us have a general position assumption that no three points are collinear. We will use $\cal A$ to denote the convex hull of the point set $P$.

Definition 1**.**

A pair of points $(p,q)$ is said to be a candidate pair if $p\in P_{b}$ and $q\in P_{r}\cup P_{b}$.

Definition 2**.**

A rectangle with one of its boundaries defined by a candidate pair, and each of the other three boundaries containing at least one point in $P_{b}$ is referred to as a candidate $LRR$, or $cLRR$ in short.

We consider each candidate pair $(p,q)$, and define a line $\ell_{pq}$ passing through $p$ and $q$. We process each side of $\ell_{pq}$ separately to compute all the $cLRR$s’ with $(p,q)$ on one of its boundaries by sweeping a line parallel to $\ell_{pq}$ among the points in $P$ in that side of $\ell_{pq}$, as stated below. After considering all the candidate pairs in $P$, we report the $LRR$. We describe the method of processing the points in $P$ above888A point $(\alpha,\beta)$ is said to be above the line $ax+by+c=0$ if $a\alpha+b\beta+c>0$; otherwise the point $(\alpha,\beta)$ is below the said line. $\ell_{pq}$. A similar method works for processing the points in $P$ below $\ell_{pq}$.

3.1 Processing a candidate pair $(p,q)$

Without loss of generality, we consider $\ell_{pq}$ as the $x$-axis, and $x(p)<x(q)$. Let $P^{\prime}$ be the array containing the subset of $P$ lying above the $x$-axis. Let $P_{b}^{\prime}$ and $P_{r}^{\prime}$ denote the blue and red point set respectively in $P^{\prime}$, $m^{\prime}=|P_{b}^{\prime}|$ and $n^{\prime}=|P_{r}^{\prime}|$. We sort the points of $P_{b}^{\prime}$ with respect to their $y$-coordinates, and construct a range tree $\cal T$ with the red points in $P_{r}^{\prime}$ considering $\ell_{pq}$ as the $x$-axis.

Observe that each $cLRR$ above $\ell_{pq}$ with $(p,q)$ on its one side corresponds to a maximal empty rectangle ($MER$) [21] among the points in $P_{b}^{\prime}$ whose bottom side is aligned with the $x$-axis and containing $(p,q)$. We sweep a horizontal line $H$ in a bottom-up manner to identify all these $cLRR$s’.

During the sweep, we maintain an interval ${\cal I}=[\alpha,\beta]$. $\cal I$ is initialized by $[x_{min},x_{max}]$ at the beginning of the sweep, where $x_{min}$ and $x_{max}$ are the points of intersection of the line $\ell_{pq}$ (the $x$-axis) with the boundary of $\cal A$. For each point $\theta\in P_{b}^{\prime}$ encountered by the sweep line, if $x(\theta)\not\in{\cal I}$, sweep proceeds to process the next point. Otherwise, we have a $cLRR$ with horizontal span $[\alpha,\beta]$, and the top boundary containing $\theta$999Needless to say, its bottom boundary contains the points $(p,q)$.. Its size is determined in $O(\log n^{\prime})$ time by performing a rectangular range counting query in $\cal T$. Now,

• •

if $x(\theta)\in[x(p),x(q)]$ then the sweep stops.

• •

otherwise,

• –

if $\alpha\leq x(\theta)\leq x(p)$ then $\alpha=x(\theta)$ is set,

• –

if $x(q)\leq x(\theta)\leq\beta$ then $\beta=x(\theta)$ is set,

and the sweep continues. Finally, after considering all the points in $P_{b}$, the sweep stops. For a detailed description of our proposed method, see Algorithm 2. A similar method is adopted for the points below $\ell_{pq}$.

Lemma 4**.**

The above algorithm computes the $LRR$ in $O(m(m+n)(m\log n+m\log m+n\log n))$ time using $O(n\log n)$ extra space.

Proof.

The space complexity follows from the space needed for maintaining the range tree $\cal T$ . We now analyze the time complexity. For each candidate pair $(p,q)$, (i) the preprocessing steps sorting of the points in $P_{b}^{\prime}$, and constructing $\cal T$ with the points in $P_{r}^{\prime}$) need $O(n^{\prime}\log n^{\prime}+m^{\prime}\log m^{\prime})$ time, and (ii) during the sweep, reporting the size of each $cLRR$ needs $O(\log n^{\prime})$ time 101010the time for the counting query for a rectangle in a range tree using fractional cascading.. Since, $O(m^{\prime})$ $cLRR$ may be reported for the candidate pair $(p,q)$, the total processing time for $(p,q)$ is $O(m^{\prime}\log n^{\prime}+m^{\prime}\log m^{\prime}+n^{\prime}\log n^{\prime})$ in the worst case. The result follows from the fact that we have considered $O(m(n+m))$ candidate pairs, $m^{\prime}=O(m)$ and $n^{\prime}=O(n)$ in the worst case. ∎

The same method is followed to compute the $LBR$. Finally, $LMR$ is reported by comparing the size of $LRR$ and $LBR$. Lemma 4 says that both the time and space complexities of our proposed algorithm for computing the $LMR$ are an improvement over those of the algorithm of Bandyapadhyay and Banik [5] for the same problem.

It needs to be mentioned that, we can implement the algorithm for the $LRR$ problem in an in-place manner by replacing range tree with the in-place implementation of 2-d tree as described in Section 2 for the range counting. Thus, the preprocessed data structure (the sorted array of $P_{b}$ and the 2-d tree for $P_{r}$) can be stored in the input array $P$ without any extra space. Using the results in Lemmata 1 and 2, we have the following result.

Theorem 1**.**

In the in-place setup, one can compute an $LMR$ in $O(m(m+n)(m\sqrt{n}+m\log m+n\log n))$ time using $O(1)$ extra space.

4 $LWR$ problem in arbitrary orientation

In this section, we consider a weighted variation of P1. Here each point in $P_{r}$ is associated with a non-zero positive weight and each point in $P_{b}$ is associated with a non-zero negative weight. Our goal is to report a rectangle $LWR$ of arbitrary orientation such that the sum of weights of the points inside that rectangle (including its boundary) is maximum among all possible rectangles in that region. Unlike problem P1, here the optimum rectangle may contain points of both the colors.

Observation 3**.**

At least one side of the $LWR$ must contain two points $p,q\in P_{r}$, and other three sides either contain a point of $P_{r}$ or is open. A point $p\in P_{r}$ may appear at a corner of the solution rectangle $LWR$. In that case, $p$ is considered to be present in both the adjacent sides of $LWR$.

We will consider all possible pairs of points $(p,q)\in P_{r}$ and define a line $\ell_{pq}$ joining $p,q$. We process each side of $\ell_{pq}$ separately to compute all the candidate $LWR$, denoted as $cLWR$, among the points in $P$ lying in that side of $\ell_{pq}$. After considering all possible pairs of points, we report $LWR$. We now describe the processing of the set of points $P^{\prime}\in P$ that lies above $\ell_{pq}$.

4.1 Processing a point-pair $(p,q)$

As earlier, assume $\ell_{pq}$ to be the $x$-axis. Consider a rectangle $R$ whose bottom side aligned with $\ell_{pq}$ (see Figure 3); the top side passing through $p_{\theta}$, left and right sides are passing through $p_{b}$ and $p_{c}$ respectively. We can measure the weight of the rectangle $R$ as follows:

Let $U=\{u_{i},i=1,2,\ldots,n\}$ be the projection of all the points on $\ell_{pq}$ having $y$-coordinate (distance from $\ell_{pq}$) less than or equal to that of $p_{\theta}$. Each member $u_{i}$ is assigned an weight equal to the weight of its corresponding point $p_{i}$. Now, compute the cumulative sum of weights $W(u_{i})$ at each projected point of $U$ from left to right. Observe that the weight of the rectangle $R$ is equal to $W(c)-W(\alpha)$, where $u_{\alpha}$ is the rightmost point in $U$ to the left of $p_{b}$.

Thus, in order to get a maximum weight rectangle with its top boundary passing through the point $p_{\theta}$ and having $(p,q)$ on its bottom boundary, we need to search for an element in $u_{\alpha}\in U$ having $x$-coordinate less that $\min(x(p_{\theta}),x(p))$ having minimum weight, and an element $u_{\beta}\in U$ having $x$-coordinate greater than $\max(x(p_{\theta}),x(q))$ having maximum weight. The weight of the rectangle with $(p,q),p_{\alpha},p_{\theta},p_{\beta}$ on its bottom, left, top and right boundaries will be $W(u_{\beta})-W(u_{\alpha})$.

We sweep a horizontal line (see Figure 3) among the points in $P^{\prime}$. During the sweep, we create a projection $u_{i}$ of each point $p_{i}\in P^{\prime}$ and assign its weight $w(u_{i})=\displaystyle w(p_{i})$, and store them in a dynamically maintained weight balanced leaf search binary tree $\cal T$ [37]. Its leaves correspond to the projections of all points that are faced by the sweep line (see Figure 4). Each internal node $u$ in $\cal T$ maintains three pieces of information, namely $EXCESS$, $MAX$ and $MIN$. $MAX$ and $MIN$ store the maximum and minimum of $W(u_{i})$ values stored in the subtree rooted at the node $u$ of $\cal T$. The $EXCESS$ field is initialized with “zero”. Each projected point $u_{j}$ at the leaf also stores the cumulative sum of weights $W(u_{j})$. During the sweep, when a new point $p_{i}\in P^{\prime}$ is faced by the sweep line, $u_{i}$ is inserted in $\cal T$. Now, for all $u_{j}$ with $x(u_{j})>x(u_{i})$, the cumulative sum of weights needs to be updated as $\hat{W}(u_{j})=W(u_{j})+w(u_{i})$. We use $EXCESS$ field to defer this update as follows.

While tracing the search path to insert $u_{i}$ (= $x(p_{i})$) in $\cal T$, if the search goes from a node $v$ to its left child, then we add $w(u_{i})$ with the $EXCESS$ field of the right child $z$ of $v$. This is in anticipation that while processing another point $u_{j}\in P^{\prime}$ if the search goes through $z$, then the $EXCESS$ field of $z$ will be propagated to the $EXCESS$ field of its two children (setting the $EXCESS$ field of $z$ to [math]).

After the insertion of $u_{i}$ in $\cal T$, we trace back up to the root of $\cal T$ and update the $MAX$ and $MIN$ fields (if necessary) of each node on the search path. If the (weight-)balance condition at any node is violated, a time linear in the size of the subtree rooted at that node is spent to rebuild that subtree in the (weight-)balanced manner.

Now, if $p_{i}\in P_{r}$, then we find the $cLWR$ of maximum weight with $(p,q)$ on its bottom boundary and $p_{i}$ on its top boundary by identifying (i) a element $u_{\alpha}\in{\cal T}$ with $W(u_{\alpha})=min\{W(u)|x(u)<min(x(p),x(p_{i}))\}$ using the $MIN$ fields of the nodes on the search path, and (ii) a point $u_{\beta}\in{\cal T}$ with $W(u_{\beta})=max\{W(u)|x(u)>max(x(q),x(p_{i}))\}$ using the $MAX$ fields of the nodes on the search path. As mentioned earlier, the weight of the rectangle on $\ell_{pq}$ with $p_{\alpha},p_{i},p_{\beta}\in P_{r}$ on its left, top, and right sides respectively, is $W(u_{\beta})-W(u_{\alpha})$. The iteration continues until all the points of $P^{\prime}$ are considered by the sweep line.

Lemma 5**.**

The $cLWR$ of maximum weight with $(p,q)$ on its one side can be computed in $O((n+m)\log(n+m))$ time.

Proof.

Follows from the fact that the amortized insertion time of a point in $\cal T$ is $O(\log n)$ [37]. While rebuilding, due to the violation of balance condition, the setting of $EXCESS$, $MIN$ and $MAX$ fields of each node can also be done in $O(|{\cal T}|)$ time, and rebuilding of $\cal T$ is needed after at least $O(\log n)$ updates [37]. ∎

The algorithm proposed above is not in-place. It uses a preprocessed data structure implemented in an $O(n)$ extra space. Lemma 5 and the fact that we need to consider $\ell_{pq}$ for each pair $p,q\in P_{r}$ suggest the following result:

Theorem 2**.**

An $LWR$ of arbitrary orientation for a set of weighted points can be computed in $O(m^{2}(n+m)\log(n+m))$ time using $O(n)$ workspace.

5 Computing largest axis-parallel monochromatic cuboid

$\mathbb{R}^{3}$

We now propose an in-place algorithm for computing a monochromatic axis-parallel cuboid with the maximum number of points. Here, the input is a set of bi-chromatic points $P=P_{r}\cup P_{b}$ inside a $3D$ axis-parallel region $\cal A$ bounded by six axis-parallel planes, where $P_{r}$ is the set of $n$ red points and $P_{b}$ is the set of $m$ blue points. The input points are given in an array, also called $P$. The $x,y,z$ coordinates of a point $p_{i}\in P$ are denoted by $x(p_{i})$, $y(p_{i})$ and $z(p_{i})$ respectively, along with its color information $c(p_{i})$ = red/blue. A cuboid is said to be a candidate for $LRC$ if its every face either coincides with a face of $\cal A$ or passes through a blue point, and its interior does not contain any blue point. Such a cuboid will be referred to as $cLRC$. The objective is to identify an $LRC$, which is a $cLRC$ containing the maximum number of red points. Similarly, a blue cuboid containing the maximum number of blue points ($LBC$) can be defined. The $LMC$ is either $LRC$ or $LBC$ depending on whose number of points is more.

We compute all possible maximal empty cuboid [34] among the $m$ blue points. Each one will be a $cLRC$; we perform a range query to count the number of red points it contains. In our algorithm, three types of $cLRC$s’ inside $\cal A$ are considered separately.

type-1:

the $cLRC$ with both top and bottom faces aligned with the top and bottom faces of $\cal A$,

type-2:

the $cLRC$ whose top face is aligned with the top face of $\cal A$, but bottom face passes through a blue point in $P_{b}$, and

type-3:

the $cLRC$ whose top face passes through some blue point in $P_{b}$. The bottom face may pass through another blue point in $P_{b}$ or may coincide with the bottom face of $\cal A$.

As a preprocessing, we first split the array $P$ into two parts, namely $P_{r}$ and $P_{b}$, such that $P_{r}=P[1,\ldots,n]$ and $P_{b}=P[n+1,\ldots,n+m]$. We construct an in-place 2-d tree $\cal T$ with the points in $P_{r}$ considering their $(x,y)$ coordinates, which will be used for the range-counting query for the $cLRC$s’. We also sort the points in $P_{b}$ in decreasing order of their $z$-coordinates. Thus, the preprocessing needs $O(m\log m+n\log n)$ time.

In [24, 27], it is proved that the number of maximal empty hyper-rectangles among a set of $n$ points in $\mathbb{R}^{d}$ is $O(n^{d})$. In the following subsections, we will analyze the processing of these three types of $cLRC$s’ in an in-place manner. The largest among the type-$i$ $cLRC$ will be referred to as type-$i$ $LRC$, for $i=1,2,3$.

5.1 Computation of type-1 $LRC$

As both the top and bottom faces of the type-1 $cLRC$s’ are aligned with the top and bottom faces of $\cal A$, if we consider the projections of the points in $P_{b}$ on the top face of $\cal A$, then each maximal empty axis-parallel rectangle ($MER$) on the top face of $\cal A$ will correspond to a type-1 $cLRC$. Thus, the problem reduces to the problem of computing all the $MER$s’ using the array $P_{b}$ in an in-place manner, and for each $MER$, count the number of points of $P_{r}$ inside the corresponding type-1 $cLRC$ using the 2-d tree $\cal T$ with the projection of points in $P_{r}$ on the top face of $\cal A$.

Lemma 6**.**

The number of type-1 $cLRC$ is $O(m^{2})$ in the worst case and the one of maximum size can be computed in $O(m^{2}\sqrt{n}+n\log n)$ time.

Proof.

The first part of the result i.e the number of type-1 $cLRC$ follows from [33]. (i) We can generate all the $MER$s with bottom boundary passing through a point $b_{i}$ on the top face of $\cal A$ using the method described in Section 3 in $O(m)$ time, and (ii) for each $MER$, the number of projected red points inside that $MER$ can be obtained in $O(\sqrt{n})$ time using the 2-d tree $\cal T$. The second part in the result follows from the fact that $\cal T$ can be generated in $O(n\log n)$ time (see Section 2.1). ∎

5.2 Computation of type-2 $LRC$

Now we describe the in-place method of computing the largest type-2 $cLRC$ whose top face is aligned with the top face of $\cal A$, but bottom face passes through a point in $P_{b}$.

We will use $p_{1},p_{2},\ldots,p_{m}$ to denote the points in $P_{b}$ in decreasing order of their $z$-coordinates. We consider each point $p_{i}\in P_{b}$ in order and compute $LRC(p_{i})$, the largest type-2 red cuboid whose bottom face passes through $p_{i}$. Let $B_{i}=\{b_{1},b_{2},\ldots,b_{i-1}\}$, $i<m$ be the set containing the projection of all the blue points having $z$-coordinate larger than $z(p_{i})$ on the plane $H(p_{i})$. Similarly, $R_{i}=\{r_{1},r_{2},\ldots\}$ are the projection of all the red points having $z$-coordinate larger than $z(p_{i})$ on the plane $H(p_{i})$. Thus, $LRC(p_{i})$ corresponds to a rectangle on the plane $H(p_{i})$ that contains $p_{i}$, but no point of $B_{i}$ in its interior, and has the maximum number of points of $R_{i}$.

As in the earlier section, we can partition the array $P$ into two contiguous blocks $P_{b}$ and $P_{r}$. The block $P_{b}$ contains all the blue points in decreasing order of their $z$-coordinates. The block $P_{r}$ contains all the red points. The blue points are processed in top-to-bottom order. Global counters $LRC$ and $MAX_{r}$ are maintained to store the $LRC$ detected so far, and its size. While processing each point $p_{i}\in P_{b}$, let $B_{i}$ denote the blue points with their $z$-coordinates greater than $z(p_{i})$. We split $P_{r}$ into two parts. The left part contains an in-place 2-d tree ${\cal T}_{i}$ with all the red points having $z$-coordinates greater than $z(p_{i})$. The right part of $P_{r}$ contains the red points with $z$-coordinates less than $z(p_{i})$. We can compute all the $MER$s using the set $B_{i}$ as in Section 5.1. For each generated $MER$ if it contains $p_{i}$ in its interior, then we perform the range counting query in ${\cal T}_{i}$ to compute the number of red points inside it. $LRC$ and $MAX_{b}$ are updated, if necessary. Thus, we have the following result.

Lemma 7**.**

The type-2 $LRC$ can be computed in $O(m^{3}\sqrt{n}+mn\log n)$ time.

Proof.

The time complexity of processing each point $p_{i}\in P_{b}$ follows from Lemma 6. Since $m$ blue points are processed, the result follows. ∎

5.3 Computation of type-3 $LRC$

Here also we use $p_{1},\ldots,p_{m}$ to denote the points in $P_{b}$ in decreasing order of their $z$-coordinates, and the algorithm processes the members in $P_{b}$ in this order. We now describe the phase of processing of a point $p_{i}\in P_{b}$. It involves generating all the type-3 $cLRC$s whose top face passes through $p_{i}$; their bottom face may pass through another blue point $p_{j}\in P_{b}$ or may coincide with the bottom face of $\cal A$. Consider the horizontal plane $H(p_{i})$ passing through $p_{i}\in P_{b}$ and sweep it downwards until it hits the bottom face of ${\cal A}$. During this phase when the sweeping plane touches $H(p_{j})$ (i.e. hits a point $p_{j}\in P_{b}$), the points inside these two horizontal planes $H(p_{i})$ and $H(p_{j})$ will participate in computing the $cLRC$s with top and bottom faces passing through $p_{i}$ and $p_{j}$, respectively. Let, $B_{ij}=\{b_{i},\ldots,b_{j}\}$ be the projections of these blue points $p_{i},\ldots,p_{j}$ $(1\leq i<j\leq n)$ on the plane $H(p_{i})$. Similarly, consider the projections $R_{ij}$ of the red points on the plane $H(p_{i})$ those lie in between the planes $H(p_{i})$ and $H(p_{j})$. Our objective is to determine a $cLRC$ corresponding to an $MER$ on the plane $H(p_{i})$ with the points in $B_{ij}$ as obstacles that contains the maximum number of points in $R_{ij}$.

In the phase of processing $p_{i}\in P_{b}$, the points of $P$ above $H(p_{i})$ does not participate in this processing. Those points of $P_{b}$ (resp. $P_{r}$) are separately stored at the beginning of the array $P_{b}$ (resp. $P_{r}$). From now onwards, by $P_{b}$ (resp. $P_{r}$) we will mean the blue (resp. red) points below $H(p_{i})$.

We consider two mutually orthogonal axis-parallel lines $x=x(p_{i})$ and $y=y(p_{i})$ on the plane $H(p_{i})$ that partition $H(p_{i})$ into four quadrants. The blue points that belong to the $\theta$-th quadrant, are denoted by $P_{b}^{\theta}$, and are stored consecutively in the array $P_{b}[i+1,\ldots,m]$. We use $m_{\theta}=|P_{b}^{\theta}|$. While processing the point $p_{j}\in P_{b}$ during the sweep in this phase, we use $B_{ij}^{\theta}$ to denote the projections of the subset of points in $P_{b}^{\theta}$ that lie between the planes $H(p_{i})$ and $H(p_{j})$, $\theta=1,2,3,4$. The members in $B_{ij}^{\theta}$ are stored in the consecutive locations of the array $P_{b}^{\theta}$ in decreasing order of their $z$-coordinates. We maintain four index variables $\chi_{\theta}$, $\theta=1,2,3,4$, where $\chi_{\theta}$ indicates the last point hit by the sweeping plane in the $\theta$-th quadrant. Thus, $p_{j}\in P_{b}\setminus(\cup_{\theta=1}^{4}B_{ij}^{\theta})$, and is obtained by comparing the $z$-coordinates of the points $\{P_{b}[\chi_{\theta}+1],\theta=1,2,3,4\}$. We will use $R_{ij}$ to denote the projection of the points in $P_{r}$ lying between $H(p_{i})$ and $H(p_{j})$. These are stored at the beginning of the array $P_{r}$.

In each quadrant $\theta$, we define the unique maximal closest stair $\textit{STAIR}_{\theta}$ around $p_{i}$ with a subset of points of $B_{ij}^{\theta}$ as in [21, 34]. The projection points of $B_{ij}^{\theta}$, that determine $\textit{STAIR}_{\theta}$, are stored at the beginning of the sub-array $P_{b}^{\theta}$ in order of their $y$- coordinates111111The remaining elements ($B_{ij}^{\theta}\setminus\textit{STAIR}_{\theta}$ are stored just after $\textit{STAIR}_{\theta}$ in a contiguous manner in $P_{b}^{\theta}$ so that the first unprocessed element in the quadrant $\theta$ is obtained at $P_{b}[\chi_{\theta}+1]$.. Thus, $\bigcup_{\theta=1}^{4}\textit{STAIR}_{\theta}$ forms an empty ortho-convex polygon $OP$ on $H(p_{i})$ (see Figure 5(a)). As a consequence, the problem of finding a type-3 $LRC$, with top and bottom faces passing through $p_{i}$ and $p_{j}$ respectively, maps to finding an $MER$ inside this ortho-convex polygon that contains $b_{j}$ and maximum number of points in the set $R_{ij}$.

Thus we need to: (i) construct the in-place 2-d tree ${\cal T}_{ij}$ with the points in $R_{ij}$, (ii) compute all maximal empty rectangles in $OP$ that contains both $b_{i}$ and $b_{j}$ (see ), (iii) for each generated maximal empty rectangle ($MER$) perform the rectangular range counting query in ${\cal T}_{ij}$, and (iv) update $OP$ by inserting $b_{j}$ in the corresponding STAIR for processing the next blue point $p_{j+1}\in P_{b}$ during this phase. The tasks (i) and (iii) performed as mentioned in Sections 2.1 and 2.2 respectively. Task (ii) is explained in Section 5.3.1 (also see Algorithm 3). Task (iv) is explained in Section 5.3.2 (also see Algorithm 5).

5.3.1 Computation of $MER(p_{i},p_{j})$

Without loss of generality, assume that $b_{j}$ (projection of $p_{j}$ on the plane $H(p_{i})$) is in the first quadrant. If $b_{j}$ is in some other quadrant, then the situation is similarly tackled.

If there exist any point in the $STAIR_{1}$ which dominates $b_{j}$, i.e., if there exist any blue point $p$ in $STAIR_{1}$ such that $x(p)<x(b_{j})$ and $y(p)<y(b_{j})$, then no axis-parallel $cLRC$ is possible whose top boundary passes through $p_{i}$ and bottom boundary passes through $p_{j}$. Therefore we assume that $b_{j}$ is not dominated by any point in $STAIR_{1}$. We now determine the subset of points in $\textit{STAIR}_{1}\cup\textit{STAIR}_{2}$ that can appear in the north boundary of an $MER$ containing both $b_{i}$ and $b_{j}$.

Let $\textit{STAIR}_{1}=\{b_{k},k=1,2\ldots,\nu_{1}\}\subseteq B_{i}^{1}$. Let $b_{\alpha}\in STAIR_{1}$ be such that $y(b_{\alpha})=\max\{b_{k}\in STAIR_{1}|y(b_{k})<y(b_{j})\}$ (i.e., the $y$-coordinate of $b_{\alpha}$ is maximum among all the points in $STAIR_{1}$ whose $y$-coordinate is lesser than the $y$-coordinate of $b_{j}$). Similarly, let $b_{\beta}\in STAIR_{1}$ be such that $y(b_{\beta})=\min\{b_{k}\in STAIR_{1}|x(b_{k})<x(b_{j})\}$ (i.e., $y$-coordinate of $b_{\beta}$ is minimum among all the points in $STAIR_{1}$ whose $x$-coordinate is lesser than the $x$-coordinate of $p_{j}$). We define $Q=\{b_{\alpha+1},b_{\alpha+2},\ldots,b_{\beta-1}\}$ $=\{b_{k}\in\textit{STAIR}_{1}|x(b_{k})>x(b_{j})~{}\text{and}~{}y(b_{k})>y(b_{j})\}$ (see Figure 5(b)).

All the axis-parallel $MER$s in $OP$ with north boundary passing through $b_{k}$, $k\in\{\alpha+1,\alpha+2,\ldots,\beta\}$ and containing $p_{i}$ in its proper interior will contain $b_{j}$ also. We draw the projections of $b_{j}$ and $b_{\beta}$ on $\textit{STAIR}_{2}$. Let these two points be $\mu$ and $\nu$, respectively. If $x(\mu)=x(\nu)$, then no point on $\textit{STAIR}_{2}$ can appear on the north boundary of a desired axis-parallel $MER$. But if $x(\mu)<x(\nu)$, then all the points $p\in\textit{STAIR}_{2}$ satisfying $x(\mu)<x(p)<x(\nu)$ can appear on the north boundary of a desired axis-parallel $MER$. In Figure 5(b), the set of points that can appear on the north boundary of an $MER$ are marked with empty dots. The method of computing an axis-parallel $MER$ with a point $p\in\textit{STAIR}_{1}\cup\textit{STAIR}_{2}$ on its north boundary is given in Algorithm 4.

5.3.2 Updating $OP$

After computing the set of axis-parallel $MER$s in $OP$ containing both the projected points $b_{i}$ and $b_{j}$ in its interior, instead of recomputing the whole ortho-convex polygon again to process the next point $p_{j+1}\in P_{b}$, we update $OP$ by inserting $b_{j}$ in the respective $STAIR$ (see Figure 6(a)).

Without loss of generality, assume that $b_{j}$ lies in the first quadrant. After inserting $b_{j}$ in $STAIR_{1}$, none of the points in $Q\in STAIR_{1}$ will participate in forming $MER$ while processing points $p_{k}\in P_{b}$ with $z(p_{k})<z(p_{j})$. So, we need to remove the members in $Q$ from $\textit{STAIR}_{1}$. This can be done by using the algorithm for stable sorting [28], where the elements in $Q$ will assume the value 1 of the given (0, 1)-valued selection function $f$, and will stably move to the end of $\textit{STAIR}_{1}$. A simple procedure for this task is given in [11] in the context of stably selecting a sorted subset. We tailored that procedure for our purpose as follows:

We maintain two index variables $\alpha$ and $\beta$; $\alpha+1$ and $\beta-1$ indicates the starting and ending positions of the $Q$, respectively. Now, two cases may arise depending on whether $|Q|=0$ or not.

$|Q|\neq 0$

: See Case 1 of Figure 6(b). Here, we need to remove $Q$ from $STAIR_{1}$ and appropriately insert $p_{j}$ into the stair. We do this by the following way:

• •

First, by swapping $b_{j}$ and $b_{\alpha+1}$, we insert $b_{j}$ in the proper position.

• •

Now, we need to move out $b_{\alpha+1},\ldots,b_{\beta-1}$ from the $STAIR_{1}$. This can be done by a sequence of swap operations: swap($P[r],P[r-(\beta-\alpha-2)]$, starting from $r=\beta$ until $r=\nu_{1}$, where $\nu_{1}$ denotes the end of $STAIR_{1}$.

• •

Finally, we set $\nu_{1}$ as $\nu_{1}-(\beta-\alpha-2)$.

$|Q|=0$

: See Case 2 of Figure 6(b). Here, we need only to insert $b_{j}$ into the stair. We do this by first swapping $(P[\nu_{1}+1],P[j])$ and then a sequence of swapping ($P[r],P[r+1]$) starting from $r=\nu_{1}$ until $r=\beta$. Finally, we set $\nu_{1}$ as $\nu_{1}+1$.

Clearly, this updating $OP$ needs $O(|P_{i}^{1}|)$ time in the worst case.

After computing the largest type-3 axis-parallel $LRC$ with $p_{i}$ on its top boundary, we need to sort the points again with respect to their $z$-coordinates for the processing of $p_{i+1}$.

Thus we have the following result:

Lemma 8**.**

The time required for processing $p_{i}$ is $O(m^{2}+C_{i}^{\prime}\sqrt{n}+mn\log n)$ in the worst case, where $C_{i}^{\prime}$ is the number of type-3 axis-parallel $LRC$s with $p_{i}$ on its top boundary.

Proof.

The worst case time required for computing $MER(p_{i},p_{j})$ is $O(|P_{ij}|+C_{ij})$, where $P_{ij}$ denotes the number of points inside the horizontal slab bounded by $H(p_{i})$ and $H(p_{j})$, and $C_{ij}$ denotes the number of axis-parallel $MER$s containing both $b_{i}$ and $b_{j}$ inside $OP$ with the projection of points $B_{ij}$ on $H(p_{i})$. In order to compute the largest type-3 axis-parallel $LRC$ with $p_{i}$ on its top boundary, we need to compute $MER(b_{i},b_{j})$ for all $j>i$, $C_{i}^{\prime}=\displaystyle\sum_{j=i+1}^{n}C_{ij}$, and $\displaystyle\sum_{j=i+1}^{n}|P_{ij}|=O((m-i)^{2})$. For each of the cuboid $C_{i}^{\prime}$, the in-place counting query in the corresponding ${\cal T}_{ij}$ requires $\sqrt{n}$ time (using Lemma 2). The last part of the time complexity follows due to the fact that, for every point $p_{j},j=i+1,\ldots,m$ we need to construct the in-place $2d$-tree. Also after the processing of each $p_{i}\in P_{b}$, the sorting step takes $O(n\log n)$ time. ∎

Lemma 6, 7 and 8 lead to the following result.

Theorem 3**.**

The worst case time complexity of our in-place algorithm for computing the axis-parallel largest monochromatic cuboid ($LMC$) is $O(m^{3}\sqrt{n}+m^{2}n\log n)$, and it takes $O(1)$ extra space.

Acknowledgment:

The authors acknowledge the valuable constructive suggestions given by the reviewer regarding the presentation of the paper.