Resizing Operators

Introduction

In ODL, resizing of a discretized function is understood as the operation of shrinking or enlarging its domain in such a way that the size of the partition cells do not change. This “constant cell size” restriction is intentional since it ensures that the underlying operation can be implemented as array resizing without resampling, thus keeping those two functionalities separate (see Resampling).

Basic setting

Let now \mathbb{R}^n with n \in \mathbb{N} be the space of one-dimensional real vectors encoding values of a function defined on an interval [a, b] \subset \mathbb{R} (see Discretizations for details). Since values are not manipulated, the generalization to complex-valued functions is straightforward.

Restriction operator

We consider the space \mathbb{R}^m for an m < n \in \mathbb{N} and define the restriction operator

(1)R : \mathbb{R}^n \to \mathbb{R}^m, \quad R(x) := (x_p, \dots, x_{p+m-1})

with a given index 0 \leq p \leq n - m - 1. Its adjoint with respect to the standard inner product is easy to determine:

\langle R(x), y \rangle_{\mathbb{R}^m}
&= \sum_{j=0}^{m-1} R(x)_j\, y_j
= \sum_{j=0}^{m-1} x_{p+j}\, y_j
= \sum_{j=p}^{p+m-1} x_j\, y_{j-p} \\
&= \sum_{i=0}^{n-1} x_i\, R^*(y)_i

with the zero-padding operator

(2)R^*(y)_i :=
\begin{cases}
y_{i-p} & \text{if } p \leq i \leq p + m - 1, \\
0       & \text{else.}
\end{cases}

In practice, this means that a new zero vector of size n is created, and the values y are filled in from index p onwards. It is also clear that the operator R is right-invertible by R^*, i.e. R R^* = \mathrm{Id}_{\mathbb{R}^m}. In fact, any operator of the form R^* + P, where P is linear and P(x)_i = 0 for i \not \in \{p, \dots, p+m-1\} acts as a right inverse for R. On the other hand, R has no left inverse since it has a non-trivial kernel (null space) \mathrm{ker} R = \{x \in \mathbb{R}^n\,|\,x_i = 0 \text{ for } i = p, \dots, p+m-1\}.

Extension operators

Now we study the opposite case of resizing, namely extending a vector. We thus choose m > n and consider different cases of enlarging a given vector x \in \mathbb{R}^n to a vector in \mathbb{R}^m. The start index is again denoted by p and needs to fulfill 0 \leq p \leq m - n - 1, such that a vector of length n “fits into” a vector of length m when starting at index p.

It should be noted that all extension operators mentioned here are of the above form R^* + P with P acting on the “outer” indices only. Hence they all act as a right inverses for the restriction operator. This property can also be read as the fact that all extension operators are left-inverted by the restriction operator R.

Moreover, the “mixed” case, i.e. the combination of restriction and extension which would occur e.g. for a constant index shift x \mapsto (0, \dots, 0, x_0, \dots, x_{n-p-1}), is not considered here. It can be represented by a combination of the two “pure” operations.

Zero padding

In this most basic padding variant, one fills the missing values in the target vector with zeros, yielding the operator

(3)E_{\mathrm{z}} : \mathbb{R}^n \to \mathbb{R}^m, \quad E_{\mathrm{z}}(x)_j :=
\begin{cases}
x_{j-p}, & \text{if } p \leq j \leq p + n - 1, \\
0      , & \text{else}.
\end{cases}

Note that this is the adjoint of the restriction operator R defined in (1). Hence, its adjoint is given by the restriction, E_{\mathrm{z}}^* = R.

Constant padding

In constant padding with constant c, the extra zeros in (3) are replaced with c. Hence, the operator performing constant padding can be written as E_{\mathrm{c}} = E_{\mathrm{z}} + P_{\mathrm{c}}, where the second summand is given by

P_{\mathrm{c}}(x) =
\begin{cases}
0      , & \text{if } p \leq j \leq p + n - 1, \\
c      , & \text{else}.
\end{cases}

Note that this operator is not linear, and its derivative is the zero operator, hence the derivative of the constant padding operator is E_{\mathrm{c}}' = E_{\mathrm{z}}.

Periodic padding

This padding mode continues the original vector x periodically in both directions. For reasons of practicability, at most one whole copy is allowed on both sides, which means that the numbers n, m and p need to fulfill p \leq n (“left” padding amount) and m - (p + n) \leq n (“right” padding amount). The periodic padding operator is then defined as

(4)E_{\mathrm{p}}(x)_j :=
\begin{cases}
x_{j-p + n}, & \text{if } 0 \leq j \leq p - 1,     \\
x_{j-p},     & \text{if } p \leq j \leq p + n - 1, \\
x_{j-p - n}, & \text{if } p + n \leq j \leq m - 1.
\end{cases}

Hence, one can at most get 3 full periods with m = 3n and p = n. Again, this operator can be written as E_{\mathrm{p}} = E_{\mathrm{z}} + P_{\mathrm{p}} with an operator

P_{\mathrm{p}}(x)_j :=
\begin{cases}
x_{j-p + n}, & \text{if } 0 \leq j \leq p - 1,     \\
0,           & \text{if } p \leq j \leq p + n - 1, \\
x_{j-p - n}, & \text{if } p + n \leq j \leq m - 1.
\end{cases}

For the adjoint of P_{\mathrm{p}}, we calculate

\langle P_{\mathrm{p}}(x), y \rangle_{\mathbb{R}^m}
&= \sum_{j=0}^{p-1} x_{j-p+n}\, y_j + \sum_{j=p+n}^{m-1} x_{j-p-n}\, y_j \\
&= \sum_{i=n-p}^{n-1} x_i\, y_{i+p-n} + \sum_{i=0}^{m-n-p-1} x_i\, y_{i+p+n} \\
&= \sum_{i=0}^{n-1} x_i\, \big( P_{\mathrm{p},1}^*(y) + P_{\mathrm{p},2}^*(y) \big)

with

P_{\mathrm{p},1}^*(y)_i :=
\begin{cases}
y_{i+p-n}, & \text{if } n - p \leq i \leq n - 1, \\
0,         & \text{else},
\end{cases}

and

P_{\mathrm{p},2}^*(y)_i :=
\begin{cases}
y_{i+p+n}, & \text{if } 0 \leq i \leq m - n - p - 1, \\
0,         & \text{else}.
\end{cases}

In practice, this means that that besides copying the values from the indices p, \dots, p+n-1 of a vector y \in \mathbb{R}^m to a new vector x \in \mathbb{R}^n, the values corresponding to the other indices are added to the vector x as follows. The first m - n - p - 1 entries of y (negative means 0) are added to the last m - n - p - 1 entries of x, in the same ascending order. The last p entries of y are added to the first p entries of x, again keeping the order. This procedure can be interpreted as “folding back” the periodized structure of y into a single period x by adding the values from the two side periods.

Symmetric padding

In symmetric padding mode, a given vector is extended by mirroring at the outmost nodes to the desired extent. By convention, the outmost values are not repeated, and as in periodic mode, the input vector is re-used at most once on both sides. Since the outmost values are not doubled, the numbers n, m and p need to fulfill the relations p \leq n - 1 (“left” padding amount) and m - (p + n) \leq n - 1 (“right” padding amount). Now the symmetric padding operator is defined as

(5)E_{\mathrm{s}}(x)_j :=
\begin{cases}
x_{p-j},      & \text{if } 0 \leq j \leq p - 1,      \\
x_{j-p},      & \text{if } p \leq j \leq p + n - 1,  \\
x_{2n-2+p-j}, & \text{if } p + n \leq j \leq m - 1.
\end{cases}

This operator is the sum of the zero-padding operator E_{\mathrm{z}} and

P_{\mathrm{s}}(x)_j :=
\begin{cases}
x_{p-j},      & \text{if } 0 \leq j \leq p - 1,      \\
0,            & \text{if } p \leq j \leq p + n - 1,  \\
x_{2n-2+p-j}, & \text{if } p + n \leq j \leq m - 1.
\end{cases}

For its adjoint, we compute

\langle P_{\mathrm{s}}(x), y \rangle_{\mathbb{R}^m}
&= \sum_{j=0}^{p-1} x_{p-j}\, y_j + \sum_{j=p+n}^{m-1} x_{2n-2+p-j}\, y_j \\
&= \sum_{i=1}^p x_i\, y_{p-i} + \sum_{i=2n-1+p-m}^{n-2} x_i\, y_{2n-2+p-i} \\
&= \sum_{i=0}^{n-1} x_i\, \big( P_{\mathrm{s},1}^*(y) + P_{\mathrm{s},2}^*(y) \big)

with

P_{\mathrm{s},1}^*(y)_i :=
\begin{cases}
y_{p-i},   & \text{if } 1 \leq i \leq p, \\
0,         & \text{else},
\end{cases}

and

P_{\mathrm{s},2}^*(y)_i :=
\begin{cases}
y_{2n-2+p-i}, & \text{if } 2n - 1 + p - m \leq i \leq n - 2, \\
0,            & \text{else}.
\end{cases}

Note that the index condition m - (p + n) \leq n - 1 is equivalent to 2n - 1 + p - m \geq 0, hence the index range in the definition of P_{\mathrm{s},2}^* is well-defined.

Practically, the evaluation of E_{\mathrm{s}}^* consists in copying the “main” part of y \in \mathbb{R}^m corresponding to the indices p, \dots, p + n - 1 to x \in \mathbb{R}^n and updating the vector additively as follows. The values at indices 1 to p are updated with the values of y mirrored at the index position p, i.e. in reversed order. The values at the indices 2n - 1 + p - m to n - 2 are updated with the values of y mirrored at the position 2n + 2 - p, again in reversed order. This procedure can be interpreted as “mirroring back” the outer two parts of the vector y at the indices p and 2n + 2 - p, adding those parts to the “main” vector.

Order 0 padding

Padding with order 0 consistency means continuing the vector constantly beyond its boundaries, i.e.

(6)E_{\mathrm{o0}}(x)_j :=
\begin{cases}
x_0,     & \text{if } 0 \leq j \leq p - 1,      \\
x_{j-p}, & \text{if } p \leq j \leq p + n - 1,  \\
x_{n-1}, & \text{if } p + n \leq j \leq m - 1.
\end{cases}

This operator is the sum of the zero-padding operator and

P_{\mathrm{o0}}(x)_j :=
\begin{cases}
x_0,     & \text{if } 0 \leq j \leq p - 1,      \\
0,       & \text{if } p \leq j \leq p + n - 1,  \\
x_{n-1}, & \text{if } p + n \leq j \leq m - 1.
\end{cases}

We calculate the adjoint of P_{\mathrm{o0}}:

\langle P_{\mathrm{o0}}(x), y \rangle_{\mathbb{R}^m}
&= \sum_{j=0}^{p-1} x_0\, y_j + \sum_{j=p+n}^{m-1} x_{n-1}\, y_j \\
&= x_0 \sum_{j=0}^{p-1} y_j + x_{n-1} \sum_{j=p+n}^{m-1} y_j \\
&= x_0 M_{\mathrm{l},0}(y) + x_{n-1} M_{\mathrm{r},0}(y)

with the zero’th order moments

M_{\mathrm{l},0}(y) := \sum_{j=0}^{p-1} y_j, \quad M_{\mathrm{r},0}(y) := \sum_{j=p+n}^{m-1} y_j.

Hence, we get

P_{\mathrm{o0}}^*(y)_i :=
\begin{cases}
M_{\mathrm{l},0}(y), & \text{if } i = 0,     \\
M_{\mathrm{r},0}(y), & \text{if } i = n - 1, \\
0,                   & \text{else},
\end{cases}

with the convention that the sum of the two values is taken in the case that $n = 1$, i.e. both first cases are the same. Hence, after constructing the restriction x \in \mathbb{R}^n of a vector y \in \mathbb{R}^m to the main part p, \dots, p + n - 1, the sum of the entries to the left are added to x_0, and the sum of the entries to the right are added to x_{n-1}.

Order 1 padding

In this padding mode, a given vector is continued with constant slope instead of constant value, i.e.

(7)E_{\mathrm{o1}}(x)_j :=
\begin{cases}
 x_0 + (j - p)(x_1 - x_0),                     & \text{if } 0 \leq j \leq p - 1,      \\
 x_{j-p},                                      & \text{if } p \leq j \leq p + n - 1,  \\
 x_{n-1} + (j - p - n + 1)(x_{n-1} - x_{n-2}), & \text{if } p + n \leq j \leq m - 1.
\end{cases}

We can write this operator as E_{\mathrm{o1}} = E_{\mathrm{o0}} + S_{\mathrm{o1}} with the order-1 specific part

S_{\mathrm{o1}}(x)_j :=
\begin{cases}
(j - p)(x_1 - x_0),                 & \text{if } 0 \leq j \leq p - 1,      \\
0,                                  & \text{if } p \leq j \leq p + n - 1,  \\
(j - p - n + 1)(x_{n-1} - x_{n-2}), & \text{if } p + n \leq j \leq m - 1.
\end{cases}

For its adjoint, we get

\langle S_{\mathrm{o1}}(x), y \rangle_{\mathbb{R}^m}
&= \sum_{j=0}^{p-1} (j - p)(x_1 - x_0)\, y_j +
   \sum_{j=p+n}^{m-1} (j - p - n + 1)(x_{n-1} - x_{n-2})\, y_j \\
&= x_0 (-M_{\mathrm{l}}(y)) + x_1 M_{\mathrm{l}}(y) +
   x_{n-2}(-M_{\mathrm{r}}(y)) + x_{n-1} M_{\mathrm{r}}(y)

with the first order moments

M_{\mathrm{l},1}(y) := \sum_{j=0}^{p-1} (j - p)\, y_j, \quad
M_{\mathrm{r},1}(y) := \sum_{j=p+n}^{m-1} (j - p - n + 1)\, y_j.

Hence, the order-1 specific operator has the adjoint

S_{\mathrm{o1}}^*(y)_i :=
\begin{cases}
-M_{\mathrm{l},1}(y), & \text{if } i = 0,     \\
M_{\mathrm{l},1}(y),  & \text{if } i = 1,     \\
-M_{\mathrm{r},1}(y), & \text{if } i = n - 2, \\
M_{\mathrm{r},1}(y),  & \text{if } i = n - 1, \\
0,                  & \text{else},
\end{cases}

with the convention of summing values for overlapping cases, i.e. if i \in \{1, 2\}. In practice, the adjoint for the order 1 padding case is applied by computing the zero’th and first order moments of y and adding them to the two outmost entries of x according to the above rule.

Generalization to arbitrary dimension

Fortunately, all operations are completely separable with respect to (coordinate) axes, i.e. resizing in higher-dimensional spaces can be written as a series of one-dimensional resizing operations. One particular issue should be mentioned with the extension operators and their adjoints, though. When extending a small, e.g., two-dimensional array to a larger size, there is an ambiguity in how the corner blocks should be handled. One possibility would be use the small array size for the extension in both axes, which would leave the corner blocks untouched (initialized to 0 usually):

../_images/resize_small.svg

However, this is not the behavior one would often want in practice. Instead, it is much more reasonable to also fill the corners in the same way the “inner” parts have been extended:

../_images/resize_large.svg

This latter behavior is implemented in the resizing operators in ODL.

The adjoint operators of these “corner-filling” resizing operator are given by reversing the unfolding pattern, i.e. by “folding in” the large array axis by axis according to the adjoint formula for the given padding mode. This way, the corners also contribute to the final result, which leads to the correct adjoint of the 2D resizing operator. Of course, the same principle can easily be generalized to arbitrary dimension.