Image noise reduction¶

Problem Overview¶

We use an annealing machine to restore the original image from an image with added noise.

We attempt to reduce noise based on the following assumptions:

The original image and the image with noise have relatively large overlaps.
In the original image, neighboring pixels often have the same color.

We use black and white images here as a simple example. Since the data of a pixel can be represented by a binary value of black and white, the value of each pixel can be expressed a binary variable. By formulating an objective function which represents the interaction between pixels expressing the above assumptions, the original image can be estimated by finding the optimal value of this function.

Constructing the Objective Function¶

We let $V$ be the set of pixels, and let $i\in V$ be the index representing each pixel. We then let the binary valued data representing the input pixels with noise be $y$, where the value of each pixel is expressed as follows:

$$ y_{i} = \left\{ \begin{align} &+1 \quad\text{(white)}\\ &-1 \quad \text{(black)} \end{align} \right. \quad i\in V\\ $$

Also, the binary Ising variables corresponding to the output pixels are represented as follows:

$$ s_{i} = \left\{ \begin{align} &+1 \quad\text{(White)}\\ &-1 \quad \text{(Black)} \end{align} \right. \quad i\in V\\ $$

Based on the assumption that the input and output images are are relatively close (i.e., there is not much noise), we need the effect that the input and output pixels have the same values. In other words, we introduce an objective function such that $s_i$ and $y_i$ become smaller when they have the same value. For example, it is given as follows:

$$ f_1 = - \sum_{i\in V} y_{i} s_{i} $$

Since the value of the above objective function decreases when $y_{i}$ and $s_{i}$ have the same value and increases when they have different values, $f_1$ takes the minimum value when $y_{i} = s_{i}$ for all $i\in V$. However, since the input image has noise on it, the noise cannot be reduced if the output image is the same as the input image.

We thus consider the assumption that neighboring pixels tend to be the same color. In other words, we introduce an objective function that reduces the value when neighboring output pixels have the same value. For example, it is given as follows:

$$ f_2 = -\sum_{(i,j)\in E} s_i s_j $$

Here, the set of adjacent pixel pairs is defined as $E$. If all the output pixels have the same value, $f_2$ takes the smallest value. However, if all pixels have the same value, every pixel will turn out to be white or black, and the information in the original image will be lost.

Therefore, by appropriately adding $f_1$ and $f_2$ together, we try to remove the pixels that are considered to be noise while making the output image close to the input image.

$$ \begin{align} f & = f_1 + \eta f_2\\ &=- \sum_{i\in V}y_is_i - \eta \sum_{(i,j)\in E}s_i s_j \end{align} $$

Here, we have introduced the parameter $\eta>0$. This allows us to adjust the relative strength of $f_1$ and $f_2$. The larger the value of $\eta$ is, the stronger the effect of the noise reduction term is.

By minimizing this objective function and interpreting the values of the Ising variables $s$ as the values of the pixels, we can obtain an image with reduced noise.

Reference¶

Loading an Image¶

First, we define a function to download an image data and a function to convert the downloaded image to an array of Ising variables.

In [ ]:

                  import os
from PIL import Image
import numpy as np
    
%matplotlib inline
import matplotlib.pyplot as plt
    
# Create an Ising array of the original image
img = Image.open("sample.png")
x = np.where(np.array(img) >= 128, 1, -1)
plt.imshow(x, cmap="gray")
plt.show()
print(x.shape)
print(x)

                

Creating an Image with Noise¶

Next, we define a function that represents noise by randomly selecting pixels and inverting their values.

In [ ]:

                  # Function to add noise to an image
def get_noisy_img_array(img_array):
    # Convert a two-dimensional array to a one-dimensional array for easier handling
    img_shape = img_array.shape
    flattened_img = img_array.flatten()
    
    # Define a function to swap the maximum and minimum values
    min_v = min(flattened_img)
    max_v = max(flattened_img)
    
    def invert_value(v):
        return min_v + max_v - v
    
    # Percentage of noise
    ratio = 0.02
    
    # Randomly select an image to have noise and invert pixels
    for idx in np.random.choice(len(flattened_img), int(ratio * len(flattened_img))):
        flattened_img[idx] = invert_value(flattened_img[idx])
    
    # Convert to the original array form
    return flattened_img.reshape(*img_shape)

                

In [ ]:

                  # Create an Ising array of image with noise
y = get_noisy_img_array(x)
plt.imshow(y, cmap="gray")
plt.show()
print(y.shape)
print(y)

                

Creating the Ising Variable Array¶

Next, we create an array s of Ising variables. If the input image data y is a $h\times w$ two-dimensional array, the Ising variable s corresponding to the output image is also a $h\times w$ two-dimensional array.

We use gen_symbols to generate the variables. We specify the type of the variable as IsingPoly here because the polynomial of the Ising variables will be the objective function in the end. Since gen_symbols can create an array of variables in the same form as the array of the input image data y, we use this feature as follows.

In [ ]:

                  from amplify import gen_symbols, IsingPoly, sum_poly
    
# Obtain the height (h) and width (w) of an image
h, w = y.shape
    
# Create an Ising variable in the form of an h x w array
s = gen_symbols(IsingPoly, h, w)

Objective Function¶

We construct the objective function using the array of the input image data $y$ and the Ising variable array $s$ corresponding to the output image.

In [ ]:

                  # Strength parameter
eta = 0.333
    
# Calculate the objective function f
    
# - \sum_{i\in V} y_{i} s_{i}
f1 = sum_poly(h, lambda i: sum_poly(w, lambda j: -s[i][j] * y[i][j]))
    
# -\sum_{(i,j)\in E} s_i s_j
f2 = sum_poly(
    h - 1, lambda i: sum_poly(w, lambda j: -s[i][j] * s[i + 1][j])
) + sum_poly(h, lambda i: sum_poly(w - 1, lambda j: -s[i][j] * s[i][j + 1]))
    
f = f1 + eta * f2

                

Setting Up the Client and Running the Machine¶

Next, we set up the client and search for the solution corresponding to the minimum value of the objective function, by an Ising machine.

In [ ]:

                  from amplify.client import FixstarsClient
from amplify import Solver
    
# Set up the client
client = FixstarsClient()
client.parameters.timeout = 1000  # Timeout is 1 second
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"  # If you use it in a local environment, please enter the access token for Amplify AE
    
# Setting up the solver and obtaining the results
solver = Solver(client)
result = solver.solve(f)

                

Obtaining the Solutions and Displaying the Results¶

Finally, we substitute the found solution to the Ising variable $s$ to obtain the data of the output image.

Comparing with the input image, we can see that the noise has been reduced.

In [ ]:

                  from amplify import decode_solution
    
# Obtain the solution
values = result[0].values
    
# Assign the solution to the Ising variable
output = decode_solution(s, values, 1)
    
plt.imshow(output, cmap="gray")  # Restored image
plt.show()
    
plt.imshow(x, cmap="gray")  # Original image
plt.show()
    
plt.imshow(y, cmap="gray")  # Image with noise
plt.show()

                