A method and related product for obtaining an optimal transport mapping within a two-dimensional region

Abstract

The present disclosure relates to a method and related products for obtaining an optimal transport map within a two-dimensional region. Among them, the method is executed by a computing device and includes: obtaining a source density of the two-dimensional region and a target density of a target region in the two-dimensional region; determining a discrete Monge-Ampère equation associated with the optimal transport map based on the source density and the target density, wherein the discrete Monge-Ampère equation includes a Brenier potential function; constructing a discrete Poisson equation according to the discrete Monge-Ampère equation; and solving the discrete Poisson equation by using a fast Fourier transform to determine the gradient of the Brenier potential function, so as to obtain the optimal transport map from the source density to the target density. By using the solution of the present disclosure, the optimal transport map can be obtained quickly and accurately.

Description

Technical Field

[0001] The present disclosure generally relates to the field of optimal transport technology. More specifically, the present disclosure relates to a method, an apparatus, and a computer-readable storage medium for obtaining an optimal transport map within a two-dimensional region. Background Art

[0002] In recent years, optimal transport has developed rapidly and is widely applied in fields such as vision, deep learning, and medical images. It has been found through research that the optimal transport (OT) map is the most economical way to transform one probability measure to another, and the total transport cost is regarded as the Wasserstein distance between the two probability measures. Thus, the OT map can be used to measure the difference between probability distributions and for virtual magnifying glasses, image registration, etc.

[0003] In practical applications, the OT map with the cost function being the square of the Euclidean distance can be reduced to solving the Monge-Ampère equation. However, due to the strong non-linearity of the Monge-Ampère equation and the high computational complexity, how to improve the computational efficiency of the OT map has become one of the core challenges. Currently, researchers have developed various algorithms to solve this problem, such as the Sinkhorn algorithm, the convex geometric variational algorithm, and the hydrodynamic algorithm. Among them, the Sinkhorn algorithm can greatly improve the optimization speed of the Kantorovich potential by adding an entropy regularization term, but it sacrifices accuracy. The convex geometric variational algorithm has a high solution accuracy, but has a complex dynamic geometric data structure and an adaptive algorithm accuracy. The hydrodynamic algorithm finds the optimal flow in space-time, which not only increases the dimension but also reduces the computational efficiency. Thus, it can be seen that the existing algorithms are difficult to effectively solve the computational efficiency of large-scale OT problems. Summary of the Invention

[0004] To at least partially solve the technical problems mentioned in the background art, the solution of the present disclosure provides a solution for obtaining an optimal transport map within a two-dimensional region. Using the solution of the present disclosure, the optimal transport map can be obtained efficiently and accurately so as to achieve the transport from the source region to the target region with the minimum transport cost. To this end, the present disclosure provides solutions in the following aspects.

[0005] In a first aspect, the present disclosure provides a method for obtaining an optimal transport map within a two-dimensional region, the method being executed by a computing device and comprising: obtaining a source density of the two-dimensional region and a target density of a target region within the two-dimensional region; determining a discrete Monge-Ampère equation associated with the optimal transport map based on the source density and the target density, wherein the discrete Monge-Ampère equation includes a Brenier potential function; constructing a discrete Poisson equation according to the discrete Monge-Ampère equation; and solving the discrete Poisson equation using a fast Fourier transform to determine the gradient of the Brenier potential function to obtain the optimal transport map from the source density to the target density.

[0006] In one embodiment, constructing the discrete Poisson equation according to the discrete Monge-Ampère equation includes: constructing an auxiliary function according to the discrete Monge-Ampère equation; and constructing the discrete Poisson equation using the auxiliary function.

[0007] In another embodiment, constructing the discrete Poisson equation using the auxiliary function includes: discretizing the target region into a Cartesian grid and initializing the auxiliary function; calculating the auxiliary function of the current iteration of the current grid point in the Cartesian grid using the initialized auxiliary function; copying the auxiliary function of the current iteration to a virtual cell region transformed via the current grid point to obtain an extended auxiliary function; and constructing the discrete Poisson equation on the Cartesian grid based on the extended auxiliary function.

[0008] In yet another embodiment, constructing the discrete Poisson equation on the Cartesian grid based on the extended auxiliary function includes: calculating a plurality of first-order partial derivatives and a plurality of second-order partial derivatives of the extended auxiliary function respectively using a finite difference algorithm on the Cartesian grid; determining the current gradient of the extended auxiliary function based on the plurality of first-order partial derivatives using a scan conversion algorithm; determining a loading function of the source density and the target density based on the plurality of second-order partial derivatives, the current gradient, the source density, and the target density; and constructing the discrete Poisson equation based on the loading function.

[0009] In yet another embodiment, determining the loading function of the source density and the target density based on the plurality of second-order partial derivatives, the current gradient, the source density, and the target density includes: determining the loading function of the source density and the target density based on the extended auxiliary function according to the following formula:

[0010]

[0011] where ρ (n)A loading function of an auxiliary function based on the extension for representing the source density and the target density Denotes the second partial derivative with respect to x Denotes the second partial derivative with respect to y Denotes the second mixed partial derivative with respect to xy, f represents the source density, and g represents the target density Denotes the current gradient, and Id denotes the identity transformation

[0012] In yet another embodiment, constructing the discrete Poisson equation based on the loading function includes: constructing the discrete Poisson equation based on the following formula

[0013]

[0014] where Δ represents the Laplace operator Denotes the auxiliary function of the next iteration, ρ (n) Denotes the loading function

[0015] In yet another embodiment, solving the discrete Poisson equation using the fast Fourier transform to determine the gradient of the Brenier potential function to obtain the optimal transport map from the two-dimensional region to the target region includes: solving the discrete Poisson equation using the fast Fourier transform to obtain the corresponding equation solution; updating the auxiliary function of the next iteration according to the corresponding equation solution to obtain the auxiliary function of the next iteration; comparing the L 2 distance between the auxiliary function of the next iteration and the current auxiliary function with a preset error value; and in response to the L 2 distance between the auxiliary function of the next iteration and the current auxiliary function being less than the preset error value, determining the sum of the gradient of the returned auxiliary function of the next iteration and the identity transformation as the gradient of the Brenier potential function to obtain the optimal transport map from the two-dimensional region to the target region

[0016] In yet another embodiment, solving the discrete Poisson equation using the fast Fourier transform to obtain the corresponding equation solution includes: adding an offset to the loading function in the discrete Poisson equation so that the integral of the loading function is zero; performing a discrete cosine transform on the loading function with the offset added using the fast Fourier transform to obtain the corresponding frequency-domain representation; performing a transformation operation on the corresponding frequency-domain representation; and using the fast Fourier transform to perform an inverse discrete cosine transform on the transformed corresponding frequency-domain representation to obtain the equation solution corresponding to the discrete Poisson equation satisfying the Neumann boundary condition

[0017] In a second aspect, the present disclosure further provides a device for obtaining an optimal transport map within a two-dimensional region, including: a processor; and a memory storing program instructions for obtaining an optimal transport map within a two-dimensional region, which, when executed by the processor, cause the device to implement the foregoing multiple embodiments.

[0018] In a third aspect, the present disclosure further provides a computer-readable storage medium storing computer-readable instructions for obtaining an optimal transport map within a two-dimensional region, which, when executed by one or more processors, implement the foregoing multiple embodiments.

[0019] Through the solution of the present disclosure, a discrete Poisson equation is constructed by means of a discrete Monge-Ampère equation determined based on a source density and a target density, and the gradient of a Brenier potential function is determined by solving the discrete Poisson equation using a fast Fourier transform, thereby obtaining an optimal transport map quickly and accurately. Further, the present disclosure improves the solution efficiency by transforming the solution of the Monge-Ampère equation into a fixed-point problem and solving, through iteration and using a fast Fourier transform, the solution of an equation of a discrete Poisson equation satisfying Neumann boundary conditions. Based on this, by using the solution of the present disclosure, technologies such as virtual magnifying glasses and image registration can be implemented, thereby making the solution of the present disclosure have characteristics such as strong versatility, strong scalability, high operation efficiency, and high operation accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] By reading the following detailed description with reference to the accompanying drawings, the above and other objects, features, and advantages of the exemplary embodiments of the present disclosure will become readily understood. In the drawings, several embodiments of the present disclosure are shown in an exemplary rather than restrictive manner, and the same or corresponding reference numerals represent the same or corresponding parts, wherein:

[0021] Figure 1 is an exemplary flowchart showing a method for obtaining an optimal transport map within a two-dimensional region according to an embodiment of the present disclosure;

[0022] Figure 2 is an exemplary overall flowchart showing a method for constructing a discrete Poisson equation by using an auxiliary function according to an embodiment of the present disclosure;

[0023] Figure 3 is an exemplary schematic diagram showing a Cartesian grid and a virtual cell region according to an embodiment of the present disclosure;

[0024] Figure 4 is an exemplary flowchart showing a method for constructing a discrete Poisson equation based on an extended auxiliary function on a Cartesian grid according to an embodiment of the present disclosure;

[0025] Figure 5is an exemplary flow chart showing the overall method of using the fast Fourier transform to solve the discrete Poisson equation to determine the gradient of the Brenier potential function according to an embodiment of the present disclosure;

[0026] Figure 6 is an exemplary flow chart showing the overall method for obtaining the optimal transport map within a two-dimensional region according to an embodiment of the present disclosure;

[0027] Figures 7 - 12 is an exemplary result graph showing the use of the embodiments of the present disclosure and other methods to process different images; and

[0028] Figure 13 is a block diagram showing a device for obtaining the optimal transport map within a two-dimensional region according to an embodiment of the present disclosure. Detailed implementation manners

[0029] The technical solutions in the embodiments of the present disclosure will be clearly and completely described below with reference to the accompanying drawings. It should be understood that the embodiments described in this specification are only some embodiments provided by the present disclosure for the convenience of clear understanding of the solutions and compliance with legal requirements, and not all embodiments that can implement the present disclosure. All other embodiments obtained by those skilled in the art based on the embodiments disclosed in this specification without creative efforts belong to the scope of protection of the present disclosure.

[0030] Figure 1 is an exemplary flow chart showing method 100 for obtaining the optimal transport map within a two-dimensional region according to an embodiment of the present disclosure. Based on the following description, those skilled in the art will understand that method 100 here can be executed by a computing device. According to different embodiments, the computing device can be a general computing device including a processor and a memory, or an artificial intelligence device including a dedicated processor (such as an artificial intelligence processor) and a memory.

[0031] As Figure 1 shown, at step S102, the source density of the two-dimensional region and the target density of the target region in the two-dimensional region are obtained. In one embodiment, the aforementioned two-dimensional region can be, for example, a two-dimensional image region, which includes but is not limited to a medical image region. In an implementation scenario, according to the given two-dimensional region, the source density of the two-dimensional region and the target density of the target region in the two-dimensional region can be obtained. Based on the aforementioned obtained source density and target density, at step S104, the discrete Monge-Ampère equation associated with the optimal transport map is determined based on the source density and the target density. In one implementation scenario, the discrete Monge-Ampère equation includes a Brenier potential function. The following will describe in detail how to determine the discrete Monge-Ampère equation.

[0032] As an example, assume that the above source density and target density are \(d\mu = f(x)dx\) and \(d\nu = g(y)dy\) respectively, and their total masses are equal, i.e., \(\int\) Ω \(f(x)dx=\int\) Ω* \(g(y)dy\). Next, assume that there exists a mapping \(T:\Omega\rightarrow\Omega\) * which is \(C\) 1 smooth and satisfies the Jacobi equation: Then it can be said that the mapping \(T\) is measure-preserving, and it is denoted as \(T\) # \(\mu=\nu\). Among them, the aforementioned \(\Omega\) represents the two-dimensional source region, and \(\Omega\) * represents the target region, and \(DT\) represents the Jacobi matrix in the aforementioned Jacobi equation. Further, assume that the cost function and \(c(x,y)\) represents the cost required to transport a unit mass from \(x\) to \(y\), then the total transport cost of this mapping is \(C(T):=\int\) Ω \(c(x,T(x))f(x)dx\). Thus, the optimal transport problem can be expressed as That is, among all measure-preserving mappings, the mapping that minimizes the total transport cost is the optimal transport mapping (i.e., the OT mapping).

[0033] Further, based on Brenier's theorem, assume Then there exists a convex function (i.e., the Brenier potential function of the embodiment of the present disclosure) such that the optimal transport mapping is equal to the gradient of the Brenier potential function: By substituting the Brenier potential function into the above Jacobi equation, the discrete Monge-Ampère equation can be obtained, which can be specifically expressed as the following formula:

[0034]

[0035] Among them, \(\det(*)\) represents the determinant of the matrix \(*\), \(D\) 2 \(u(x)\) represents the second derivative of \(u(x)\), \(f(x)\) represents the source density, \(g\) represents the target density, represents the gradient of the Brenier potential function, and \(\circ\) represents the composition of functions. For formula (1), it has the second type of boundary condition: Thus, assume that the boundary of the support set and is \(C\) 1 smooth, then the optimal mapping satisfies the tilt condition, that is, the outer normal vector at the boundary point of the region satisfies \(\langle n(x),n(T(x))\rangle>0\), and That is, the inner product of the outer normal vector \(n(x)\) at the boundary point of the two-dimensional region and the outer normal vector \(n(T(x))\) at the boundary point of the target region is non-negative.

[0036] Based on the determined discrete Monge-Ampère equation above, at step S106, a discrete Poisson equation is constructed according to the discrete Monge-Ampère equation. In one embodiment, first, an auxiliary function can be constructed according to the discrete Monge-Ampère equation, and then the auxiliary function is used to construct the discrete Poisson equation. In the implementation scenario, considering the above two-dimensional region Ω and the target region Ω * are the same planar rectangular region. Under the satisfaction of the second boundary condition and the inclination condition, for any boundary point that is not a corner point its image point is also a boundary point that is not a corner point, and satisfies n(x) = n(T(x)). This means that both x and T(x) are on the same side of the rectangle. Therefore, the image of the corner point is equal to itself. Based on this, the above Brenier potential function can be rewritten into the following formula to construct the auxiliary function:

[0037]

[0038] The auxiliary function constructed according to formula (2) Furthermore, the auxiliary function can be used to construct the discrete Poisson equation. For example, the target region can be first discretized into a Cartesian grid, and then the discrete Poisson equation is constructed based on the auxiliary function within the Cartesian grid. As an example, the foregoing can be the fixed point of the operator That is, by transforming the discrete Monge-Ampère equation into a fixed-point problem, and then obtaining the solution of the discrete Poisson equation through iteration to obtain the Brenier potential function, so as to obtain the optimal transport mapping.

[0039] Specifically, first, the target region can be discretized into a Cartesian grid and the auxiliary function is initialized. Then, the auxiliary function of the current iteration of the current grid point in the Cartesian grid is calculated using the initialized auxiliary function. Further, the auxiliary function of the current iteration is copied to the virtual cell region transformed by the current grid point to obtain an extended auxiliary function. Then, the discrete Poisson equation is constructed based on the extended auxiliary function on the Cartesian grid. More specifically, on the Cartesian grid, the finite difference algorithm is used to calculate multiple first-order partial derivatives and multiple second-order partial derivatives of the extended auxiliary function respectively, and the scan conversion algorithm is used to determine the current gradient of the extended auxiliary function based on multiple first-order partial derivatives. Then, based on multiple second-order partial derivatives, the current gradient, the source density, and the target density, the loading function of the source density and the target density is determined, and the discrete Poisson equation is constructed based on the loading function. The discrete Poisson equation includes the Laplace operator, the auxiliary function of the next iteration, and the loading function. The construction of the discrete Poisson equation will be described in detail later in combination with Figures 2 - 4

[0040] ​After obtaining the above discrete Poisson equation, at step S108, the discrete Poisson equation is solved using the fast Fourier transform to determine the gradient of the Brenier potential function, so as to obtain the optimal transport map from the source density to the target density. That is, the obtained gradient of the Brenier potential function is the optimal transport map. Specifically, first, the discrete Poisson equation can be solved using the fast Fourier transform to obtain the corresponding equation solution, and then the auxiliary function for the next iteration is updated according to the corresponding equation solution to obtain the auxiliary function for the next iteration. Further, the L 2 distance between the auxiliary function for the next iteration and the current auxiliary function is compared with a preset error value, so as to determine the gradient of the Brenier potential function based on the comparison result, so as to obtain the optimal transport map from the source density defined on a two-dimensional region to the target density defined on a target region. In one embodiment, the following operations can be performed to solve the discrete Poisson equation using the fast Fourier transform to obtain the corresponding equation solution: that is, an offset is added to the loading function in the discrete Poisson equation to make the integral of the loading function zero, and then the discrete cosine transform is performed on the loading function after adding the offset to obtain the corresponding frequency domain representation. Then, a transformation operation is performed on the corresponding frequency domain representation, and the inverse discrete cosine transform is performed on the transformed corresponding frequency domain representation to obtain the equation solution corresponding to the discrete Poisson equation that satisfies the Neumann boundary condition. How to solve the discrete Poisson equation using the fast Fourier transform will be described in detail later in conjunction with Figure 5 how to solve the discrete Poisson equation using the fast Fourier transform will be described in detail later in conjunction with

[0041] Combined with the above description, it can be seen that the present disclosure transforms the discrete Monge-Ampère problem between the source density and the target density into a fixed-point problem, and solves the discrete Poisson equation through iteration and the fast Fourier transform, thereby greatly improving the calculation efficiency and calculation accuracy. In practical applications, for the densities of any two two-dimensional regions that satisfy the boundary conditions, the optimal transport map can be quickly and accurately obtained using the solution of the present disclosure, making the solution of the present disclosure highly versatile and scalable. Based on the solution of the present disclosure, for example, a virtual magnifying glass (such as magnifying an image) and image registration can be realized.

[0042] Figure 2 is an exemplary flowchart showing the overall method 200 of constructing a discrete Poisson equation using an auxiliary function according to an embodiment of the present disclosure. It should be understood that Figure 2 the method 200 in Figure 1 is a specific embodiment of step S106 of the method 100 in the above Figure 1 above, so the above description about Figure 2 also applies to

[0043] As Figure 2As shown, at step S202, the target region is discretized into a Cartesian grid and the auxiliary function is initialized. In one embodiment, by treating each pixel block in the target region as a point, the target region can be discretized into a Cartesian grid composed of multiple grid points. In this scenario, assuming the size of the Cartesian grid is M*N, the step sizes h x = 2 / M and h y = 2 / N can be obtained, and then the coordinates of each grid point (i, j) can be obtained In addition, at the beginning of the iteration, the auxiliary function can be initialized, that is, is set to 0. Then, at step S204, the auxiliary function of the current iteration of the current grid point in the Cartesian grid is calculated using the initialized auxiliary function. In one embodiment, the coordinates of the current grid point can be substituted into the initialized auxiliary function, and the auxiliary function of the current iteration of the current grid point can be obtained using the above formula (2). Based on the obtained auxiliary function of the current iteration, at step S206, the auxiliary function of the current iteration is copied to the virtual cell region transformed via the current grid point to obtain an extended auxiliary function. In one implementation scenario, the aforementioned virtual cell region refers to the region formed by adding one row above, below, to the left, and to the right of the current grid point in the Cartesian grid (for example Figure 3 as shown). In this scenario, the auxiliary function of the current iteration is copied into the virtual cell region to form an extended auxiliary function, and by adding the aforementioned virtual cell region, it is beneficial to calculate the higher-order derivatives of the extended auxiliary function under Neumann conditions. After obtaining the extended auxiliary function, at step S208, a discrete Poisson equation is constructed on the Cartesian grid based on the extended auxiliary function. How to construct a discrete Poisson equation on the Cartesian grid based on the extended auxiliary function will be described in detail later in conjunction with Figure 4

[0044] Figure 3 is an exemplary schematic diagram showing a Cartesian grid and a virtual cell region according to an embodiment of the present disclosure. As Figure 3 ​As shown, the multiple square blocks shown in the figure are pixel blocks of the target area. As described above, by treating each pixel block as a point, a Cartesian grid composed of multiple grid points can be formed. For example, the multiple hollow dots connected by solid lines shown in the figure form a Cartesian grid, and the hollow dots are the multiple grid points of the Cartesian grid, that is, the current grid points of the Cartesian grid. The figure further shows that one row of grid points is added above, below, to the left, and to the right of the Cartesian grid, thereby forming a virtual cell area. For example, the area formed by the multiple hollow dots connected by dashed lines shown in the figure is the virtual cell area. Further, by copying the auxiliary function of the current iteration to this virtual cell area, an extended auxiliary function is obtained, and then a discrete Poisson equation can be constructed based on the extended auxiliary function.

[0045] Figure 4 FIG. is an exemplary flowchart showing a method 400 for constructing a discrete Poisson equation based on an extended auxiliary function on a Cartesian grid according to an embodiment of the present disclosure. It should be understood that Figure 4 The method 400 in Figure 2 is a specific embodiment of step S208 of the method 200 in the above Figure 2 Therefore, the description made above regarding Figure 4 also applies to

[0046] As Figure 4 shown, at step S402, on the Cartesian grid, the finite difference algorithm is used to calculate the multiple first-order partial derivatives and the multiple second-order partial derivatives of the extended auxiliary function respectively. In one embodiment, assume that the extended auxiliary function is denoted as Then the aforementioned multiple first-order partial derivatives may include the first-order partial derivative of the extended auxiliary function with respect to x and the first-order partial derivative of the extended auxiliary function with respect to y. As an example, denote the first-order partial derivative of the extended auxiliary function with respect to x as Denote the first-order partial derivative of the extended auxiliary function with respect to y as where n represents the current iteration step. Similarly, the aforementioned multiple second-order partial derivatives may include the second-order partial derivative of the extended auxiliary function with respect to x, the second-order partial derivative with respect to y, and the second-order mixed partial derivative with respect to xy, and they are denoted as and In an implementation scenario, the aforementioned and can be calculated by the following formulas:

[0047]

[0048]

[0049]

[0050] According to the foregoing, hx and h y represent the step sizes of the Cartesian grid in the horizontal and vertical directions respectively, where \(i = 0,\cdots,M - 1\), \(j = 0,\cdots,N - 1\), and \(M\) and \(N\) represent the number of rows and columns of the Cartesian grid respectively. Based on the above formulas (3)-(5), a difference equation can also be obtained, which is specifically expressed as follows:

[0051]

[0052] According to the multiple first-order partial derivatives obtained above (the first-order partial derivative with respect to \(x\) and the first-order partial derivative with respect to \(y\) ), at step S404, a scan conversion algorithm is used to determine the current gradient of the extended auxiliary function based on the multiple first-order partial derivatives. As is known to those skilled in the art, the scan conversion algorithm is a mathematical algorithm in computer graphics to convert two-dimensional or three-dimensional graphics into the raster form of a computer monitor. In an implementation scenario, the current gradient of the extended auxiliary function can be obtained through the scan conversion algorithm, for example, denoted as Furthermore, at step S406, a loading function of the source density and the target density is determined based on the multiple second-order partial derivatives, the current gradient, the source density, and the target density. As mentioned above, can be a fixed point of the operator , and thus this operator satisfies the following formula:

[0053]

[0054] Based on this, a loading function \(\rho\) of the source density and the target density based on the extended auxiliary function can be determined (n) , and its specific formula can be expressed as follows:

[0055]

[0056] where \(\rho\) (n) represents the loading function of the source density and the target density based on the extended auxiliary function, represents the second-order partial derivative with respect to \(x\), represents the second-order partial derivative with respect to \(y\), represents the second-order mixed partial derivative with respect to \(xy\), \(f\) represents the source density, \(g\) represents the target density, represents the current gradient, and \(Id\) represents the identity transformation.

[0057] Furthermore, at step S408, a discrete Poisson equation is constructed based on the loading function. In one embodiment, the discrete Poisson equation can be expressed by the following formula:

[0058]

[0059] Among them, Δ represents the Laplace operator, represents the auxiliary function for the next iteration, and ρ (n) represents the loading function. Then, by solving the discrete Poisson equation using the fast Fourier transform during the iteration process, the gradient of the Brenier potential function is determined, and finally the optimal transport mapping from the two-dimensional region to the target region is obtained.

[0060] Figure 5 is an exemplary flowchart showing the overall method 500 for determining the gradient of the Brenier potential function by solving the discrete Poisson equation using the fast Fourier transform according to an embodiment of the present disclosure. It should be understood that, Figure 5 the method 500 in Figure 1 is a specific embodiment of step S108 of the method 100 in the above Figure 1 , so the above description about Figure 5 also applies to

[0061] As Figure 5 shown, at step S502, the discrete Poisson equation is solved using the fast Fourier transform to obtain the corresponding equation solution. In one embodiment, first, an offset can be added to the loading function in the discrete Poisson equation to make the integral of the loading function zero, for example, as shown in the following formula:

[0062]

[0063] where ρ represents the loading function, represents the auxiliary function, and Δ represents the Laplace operator. In the implementation scenario, assuming the offset is c, then ∫ Ω* ρ - c = 0.

[0064] Next, the discrete cosine transform (Discrete Cosine Transform, "DCT") is performed on the loading function with the added offset using the fast Fourier transform to obtain the corresponding frequency-domain representation. Specifically, the DCT can be performed on the loading function with the added offset using the fast Fourier transform through the following formula to obtain the frequency-domain representation of the loading function with the added offset:

[0065]

[0066]

[0067] where is the frequency-domain representation of the loading function with the added offset, m, i = 0,..., M - 1, n, j = 0,..., N - 1, and M and N respectively represent the number of rows and columns of the Cartesian grid.

[0068] Based on the frequency-domain representation of the loading function with the added offset obtained above, a transformation operation can be performed on the corresponding frequency-domain representation. In one embodiment, the transformation operation can be performed on the corresponding frequency-domain representation by the following formula:

[0069]

[0070] where, represents the frequency-domain representation of the loading function with the added offset, is the transformed frequency-domain representation, m, i = 0,..., M - 1, n, j = 0,..., N - 1, and M and N respectively represent the number of rows and columns of the Cartesian grid.

[0071] Further, the inverse discrete cosine transform is performed on the corresponding transformed frequency-domain representation using the fast Fourier transform to obtain the solution of the equation corresponding to the discrete Poisson equation satisfying the Neumann boundary condition. Specifically, the inverse DCT can be performed on the corresponding transformed frequency-domain representation by the following formula and using the fast Fourier transform:

[0072]

[0073]

[0074] where, represents the solution of the equation corresponding to the discrete Poisson equation and satisfies the Neumann boundary condition, represents the corresponding transformed frequency-domain representation, m, i = 0,..., M - 1, n, j = 0,..., N - 1, and M and N respectively represent the number of rows and columns of the Cartesian grid.

[0075] After obtaining the above equation solution, at step S504, the auxiliary function for the next iteration is updated according to the corresponding equation solution to obtain the auxiliary function for the next iteration. That is, the corresponding equation solution is substituted into the above formula (9) for updating to obtain the auxiliary function for the next iteration. At step S506, the L 2 distance between the auxiliary function for the next iteration and the current auxiliary function is compared with a preset error value. For example, assume that the auxiliary function for the next iteration is denoted as the current auxiliary function is denoted as the preset error value is denoted as ε, then the L 2 distance between the auxiliary function for the next iteration and the current auxiliary function is compared to see if it is less than the preset error value, that is, to determine whether the inequality holds. At step S508, in response to the L 2If the distance is less than a preset error value, the sum of the gradient of the auxiliary function for the next iteration returned and the identity transformation is determined as the gradient of the Brenier potential function to obtain the optimal transport map from the source density to the target density. That is, when the L 2 distance between the auxiliary function for the next iteration and the current auxiliary function is less than the preset error value, that is holds, the sum of the gradient of the auxiliary function for the next iteration and the identity transformation is returned. For example, return where Id is the identity transformation, and is the gradient of the auxiliary function for the next iteration. In this scenario, the returned 2 is the gradient of the Brenier potential function, that is, the optimal transport map from the source density to the target density. On the contrary, when the L distance between the auxiliary function for the next iteration and the current auxiliary function is not less than the preset error value, that is Figures 2 - 5 does not hold, then the method described above is repeated based on the auxiliary function for the next iteration for iteration until

[0076] Figure 6

[0077] Figure 6 is an exemplary flowchart showing the overall process for obtaining the optimal transport map within a two-dimensional region according to an embodiment of the present disclosure. As

[0077] shown in Set to 0. According to the initialized auxiliary function, the auxiliary function of the current iteration can be obtained. Then, at step S606, the auxiliary function of the current iteration is copied to the virtual cell region to obtain an extended auxiliary function. Based on the obtained extended auxiliary function, at step S608, a plurality of first-order partial derivatives (e.g., ) of the extended auxiliary function are calculated respectively in the Cartesian grid using the finite difference method and the scan conversion algorithm, a plurality of second-order partial derivatives ( and ) and the current gradient. Among them, the aforementioned plurality of second-order partial derivatives can be calculated based on the above formulas (3)-(5). At step S610, a loading function ρ (n) of the source density and the target density can be determined according to the aforementioned plurality of second-order partial derivatives, the current gradient, the source density, and the target density, such as shown in the above formula (8). Then, at step S612, a discrete Poisson equation is constructed based on the loading function

[0078] After obtaining the above discrete Poisson equation, at step S614, the discrete Poisson equation can be solved using the fast Fourier transform, and the auxiliary function of the next iteration is updated using its corresponding equation solution to obtain the auxiliary function of the next iteration at step 616 In one embodiment, solving the discrete Poisson equation using the fast Fourier transform can refer to step S502 of method 500 described above Figure 5 in this disclosure, which will not be elaborated here. Then, at step S618, the L distance between the auxiliary function of the next iteration and the current auxiliary function 2 is compared with a preset error value ε to determine whether holds. When holds, at step S620, is returned as the optimal transport map. When does not hold, return to the aforementioned step S606 to continue the iteration until holds and the iteration stops.

[0079] Figures 7 - 12 is an exemplary result graph showing the processing of different images using the embodiments of the present disclosure and other methods. As shown in Figures 7 - 12 respectively show the Buddha image, the human brain image, the David image, the old person image, the female face image, and the male face image, and Figures 7 - 12The upper left diagrams all show the original images of the corresponding images (i.e., the Buddha image, the human brain image, the David image, the old person image, the female face image, and the male face image), and the upper right diagrams all show the result diagrams of the corresponding images based on the Kantorovich potential energy (i.e., the auxiliary function). Further, Figures 7 - 12 The lower left diagrams all show the result diagrams of the corresponding images based on conformal mapping, and the lower right diagrams all show the result diagrams of the corresponding images obtained based on the embodiments of the present disclosure (i.e., based on the optimal transport mapping). It can be seen that the solution of the embodiments of the present disclosure can more accurately achieve the transmission from the source region to the target region, making the transmitted image information more complete.

[0080] In addition, the embodiments of the present disclosure also compare the calculation rates with other algorithms (such as the convex geometric optimization algorithm), thereby indicating that the embodiments of the present disclosure can improve the calculation efficiency. Specifically, by using the convex geometric optimization algorithm and the solution of the embodiments of the present disclosure to process the Figures 7 - 12 original images with a resolution of 512×512 and the original images with a resolution of 1k×1k in the middle respectively, and obtaining the corresponding calculation times, as shown in Table 1 and Table 2 below.

[0081] Table 1 Efficiency detection of the original image with a resolution of 512×512

[0082]

[0083] Table 2 Efficiency detection of the original image with a resolution of 1k×1k

[0084]

[0085] As can be seen from Table 1 and Table 2 above, compared with the convex geometric algorithm, using the solution of the embodiments of the present disclosure to obtain the optimal transport mapping, the calculation rate is increased by dozens to hundreds of times. In addition, the convex geometric algorithm cannot process the human brain image, while the solution of the embodiments of the present disclosure can, so the solution of the embodiments of the present disclosure not only has high calculation efficiency, but also has strong versatility, strong scalability and high stability.

[0086] Figure 13 is a block diagram showing a device 1300 for obtaining the optimal transport mapping within a two-dimensional region according to the embodiments of the present disclosure. It can be understood that the device implementing the solution of the present disclosure can be a single device (such as a computing device) or a multifunctional device including various peripheral devices.

[0087] Such as Figure 13As shown, the device of the present disclosure may include a central processing unit or central processing unit (“CPU”) 1311, which may be a general-purpose CPU, a dedicated CPU, or other information processing and program execution units. Further, the device 1300 may also include a mass storage 1312 and a read-only memory (“ROM”) 1313, where the mass storage 1312 may be configured to store various types of data, including various image regions and their densities, algorithm data, intermediate results, and various programs required to run the device 1300. The ROM 1313 may be configured to store data and instructions for power-on self-test of the device 1300, initialization of each functional module in the system, basic input / output drivers of the system, and data and instructions required to boot the operating system.

[0088] Optionally, the device 1300 may also include other hardware platforms or components, such as the shown tensor processing unit (“TPU”) 1314, graphics processing unit (“GPU”) 1315, field programmable gate array (“FPGA”) 1316, and machine learning unit (“MLU”) 1317. It can be understood that although various hardware platforms or components are shown in the device 1300, they are merely exemplary rather than restrictive, and those skilled in the art can add or remove corresponding hardware according to actual needs. For example, the device 1300 may only include a CPU, related storage devices, and interface devices to implement the method of the present disclosure for obtaining the optimal transport map within a two-dimensional region.

[0089] In some embodiments, to facilitate the transfer and interaction of data with an external network, the device 1300 of the present disclosure further includes a communication interface 1318, so that it can be connected to a local area network / wireless local area network (“LAN / WLAN”) 1305 through the communication interface 1318, and then can be connected to a local server 1306 or connected to the Internet (“Internet”) 1307 through the LAN / WLAN. Alternatively or additionally, the device 1300 of the present disclosure may also be directly connected to the Internet or a cellular network based on wireless communication technology through the communication interface 1318, such as based on the 3rd generation (“3G”), 4th generation (“4G”), or 5th generation (“5G”) wireless communication technology. In some application scenarios, the device 1300 of the present disclosure may also access a server 1308 and a database 1309 of an external network as needed to obtain various known image models, data, and modules, and may remotely store various data, such as various types of data or instructions for presenting, for example, iteration results, optimal transport map results, etc.

[0090] The peripheral devices of device 1300 may include a display device 1302, an input device 1303, and a data transmission interface 1304. In one embodiment, the display device 1302 may include, for example, one or more speakers and / or one or more visual displays, which are configured to provide voice prompts and / or display image videos for the transmission from the two-dimensional region to the target region of the present disclosure. The input device 1303 may include, for example, a keyboard, a mouse, a microphone, a gesture capture camera, and other input buttons or controls, which are configured to receive the input of the corresponding densities of the two-dimensional region and the target region and / or user instructions. The data transmission interface 1304 may include, for example, a serial interface, a parallel interface, or a Universal Serial Bus interface (“USB”), a Small Computer System Interface (“SCSI”), Serial ATA, FireWire, PCI Express, and a High-Definition Multimedia Interface (“HDMI”), etc., which are configured for data transmission and interaction with other devices or systems. According to the solution of the present disclosure, the data transmission interface 1304 may receive medical images collected by, for example, CT and two-dimensional region images collected by other image acquisitions, and transmit to device 1300 two-dimensional region images or various other types of data or results.

[0091] The above-mentioned CPU 1311, mass storage 1312, ROM 1313, TPU 1314, GPU 1315, FPGA 1316, MLU 1317, and communication interface 1318 of device 1300 of the present disclosure may be interconnected with each other through a bus 1319, and realize data interaction with the peripheral devices through this bus. In one embodiment, through this bus 1319, the CPU 1311 may control other hardware components in device 1300 and its peripheral devices.

[0092] The above combination Figure 13 has described the device that can be used to execute the method for obtaining the optimal transport mapping within a two-dimensional region of the present disclosure. It should be understood that the device structure or architecture here is only exemplary, and the implementation manner and implementation entity of the present disclosure are not limited by it, but may be changed without departing from the spirit of the present disclosure.

[0093] According to the above description with reference to the drawings, those skilled in the art can also understand that the embodiments of the present disclosure can also be implemented by software programs. Therefore, the present disclosure also provides a computer program product. This computer program product can be used to implement the method for obtaining the optimal transport mapping within a two-dimensional region described in the present disclosure in combination with the attached Figures 1 - 6 drawings.

[0094] It should be noted that although the operations of the method of the present disclosure are described in a specific order in the accompanying drawings, this does not require or imply that these operations must be performed in that specific order, or that all the operations shown must be performed to achieve the desired result. On the contrary, the steps depicted in the flowchart can be changed in the order of execution. Additionally or alternatively, certain steps may be omitted, multiple steps may be combined into one step for execution, and / or one step may be decomposed into multiple steps for execution.

[0095] It should be understood that when terms such as "first", "second", "third", and "fourth" are used in the claims, the specification, and the accompanying drawings of the present disclosure, they are only used to distinguish different objects and not to describe a specific order. The terms "comprising" and "including" used in the specification and claims of the present disclosure indicate the presence of the described features, wholes, steps, operations, elements, and / or components, but do not exclude the presence or addition of one or more other features, wholes, steps, operations, elements, components, and / or their combinations.

[0096] It should also be understood that the terms used in the specification of the present disclosure are only for the purpose of describing specific embodiments and are not intended to limit the present disclosure. As used in the specification and claims of the present disclosure, unless the context clearly indicates otherwise, the singular forms "a", "an", and "the" are intended to include the plural forms. It should be further understood that the term "and / or" used in the specification and claims of the present disclosure refers to any combination and all possible combinations of one or more of the associated listed items, and includes these combinations.

[0097] Although the embodiments of the present disclosure are as above, the above content is only examples adopted for the convenience of understanding the present disclosure and is not used to limit the scope and application scenarios of the present disclosure. Any person skilled in the art within the technical field of the present disclosure can make any modifications and changes in the form of implementation and details without departing from the spirit and scope disclosed by the present disclosure. However, the scope of patent protection of the present disclosure shall still be subject to the scope defined by the appended claims.

Claims

1. A method for obtaining an optimal transport map within a two-dimensional region, where the two-dimensional region is a region in a two-dimensional image, the method is executed by a computing device, and includes: Obtaining a source density of the two-dimensional region and a target density of a target region in the two-dimensional region; Determining a discrete Monge-Ampère equation associated with the optimal transport map based on the source density and the target density, where the discrete Monge-Ampère equation includes a Brenier potential function; Constructing a discrete Poisson equation according to the discrete Monge-Ampère equation; and Solving the discrete Poisson equation using the fast Fourier transform to determine the gradient of the Brenier potential function to obtain the optimal transport map from the source density to the target density, where constructing the discrete Poisson equation according to the discrete Monge-Ampère equation includes: Constructing an auxiliary function according to the discrete Monge-Ampère equation; Discretizing the target region into a Cartesian grid and initializing the auxiliary function; Calculating the auxiliary function of the current iteration of the current grid point in the Cartesian grid using the initialized auxiliary function; Copying the auxiliary function of the current iteration to a virtual cell region transformed via the current grid point to obtain an extended auxiliary function; Calculating a plurality of first-order partial derivatives and a plurality of second-order partial derivatives of the extended auxiliary function respectively using a finite difference algorithm on the Cartesian grid; Determining the current gradient of the extended auxiliary function based on the plurality of first-order partial derivatives using a scan conversion algorithm; Determining a loading function of the source density and the target density based on the plurality of second-order partial derivatives, the current gradient, the source density, and the target density; and Constructing the discrete Poisson equation based on the loading function.

2. The method according to claim 1, wherein determining the loading function of the source density and the target density based on the plurality of second-order partial derivatives, the current gradient, the source density, and the target density includes: Determining the loading function of the source density and the target density based on the extended auxiliary function according to the following formula: where ρ (n) represents a loading function of the source density and the target density based on the extended auxiliary function, represents the second-order partial derivative with respect to x, represents the second-order partial derivative with respect to y, represents the second-order mixed partial derivative with respect to xy, f represents the source density, and h represents the target density, represents the current gradient, and Id represents the identity transformation.

3. The method according to claim 2, wherein constructing the discrete Poisson equation based on the loading function includes: Constructing the discrete Poisson equation according to the following formula: where Δ represents the Laplace operator, denotes the auxiliary function for the next iteration, and ρ (n) represents the said loading function.

4. The method according to claim 3, where solving the discrete Poisson equation using the fast Fourier transform to determine the gradient of the Brenier potential function to obtain the optimal transport map from the source density to the target density includes: Solving the discrete Poisson equation using the fast Fourier transform to obtain a corresponding equation solution; Updating the auxiliary function of the next iteration according to the corresponding equation solution to obtain the auxiliary function of the next iteration; Compare the L 2 distance between the auxiliary function of the next iteration and the current auxiliary function with a preset error value; And The L between the auxiliary function for the next iteration and the current auxiliary function 2 is less than the preset error value, and the sum of the gradient of the auxiliary function for the next iteration returned and the identity transformation is determined as the gradient of the Brenier potential function, so as to obtain the optimal transport mapping from the source density to the target density.

5. The method according to claim 1, where solving the discrete Poisson equation using the fast Fourier transform to obtain a corresponding equation solution includes: Adding an offset to the loading function in the discrete Poisson equation so that the integral of the loading function is zero; Performing a discrete cosine transform on the loading function with the offset added using the fast Fourier transform to obtain a corresponding frequency domain representation; Performing a transformation operation on the corresponding frequency domain representation; And Performing an inverse discrete cosine transform on the transformed corresponding frequency domain representation using the fast Fourier transform to obtain the equation solution corresponding to the discrete Poisson equation satisfying the Neumann boundary condition.

6. A device for obtaining an optimal transport map within a two-dimensional region, including: A processor; And A memory stores program instructions for obtaining an optimal transport map within a two-dimensional region. When the program instructions are executed by the processor, the device implements the method according to any one of claims 1-5.

7. A computer-readable storage medium stores computer-readable instructions for obtaining an optimal transport map within a two-dimensional region. When the computer-readable instructions are executed by one or more processors, the method according to any one of claims 1-5 is implemented.

Images