4

I have two 2D arrays: $A$ is the original matrix that contains only 0s and 1s, and $B$ is the convolved matrix. I know the size of the convolution kernel $K$. Generally, it follows $B = A*K + S$, where $S$ is the additional noise with Gaussian distribution.

I want to estimate the kernel based on $A$ and $B$, how can I do it in Python?

Here is my example code:

import numpy as np
from scipy.signal import convolve2d,oaconvolve,fftconvolve

the dimension of the original matrix


N = 100


create an example original matrix with around 20% 1s


src_mat = np.random.binomial(n=1, p=0.2, size=(N, N))


a kernel function with a dimension of 2*N-1, the values decay outwards from center of the matrix, as a function of distance (here 1/r)


kernel = np.fromfunction(lambda x,y:np.sqrt((x-(N-1))2+(y-(N-1))2)*(-1) ,(2N -1,2*N-1))
kernel[N-1,N-1] = 0


the convolved matrix in the form of B=A*K+S


tgt_mat = oaconvolve(src_mat,kernel,mode='same') + np.random.normal(0,1,src_mat.shape)


```
lennon310
  • 3,590
  • 19
  • 24
  • 27
learner_lyf
  • 43
  • 5
  • 1
    Does this answer your question? Estimating Convolution Kernel from Input and Output Images – Gillespie Aug 25 '22 at 12:38
  • Not really, I don't understand some of the derivations in the accepted answer, I also don't have any knowledge in Matlab. I reached Royi and we had a bit discussion, then I decided to open a new question. – learner_lyf Aug 25 '22 at 12:53
  • I will have a look at it. But for sure don't go as is to the frequency domain :-). – Royi Aug 25 '22 at 13:41
  • @Royi, many thanks, hope this won't take too much effort. I also prefer a similar solution as you put there. – learner_lyf Aug 25 '22 at 14:16
  • Could you write the formula for the kernel and matrix as pure math in addition to the code? – Royi Aug 27 '22 at 18:17
  • Are you sure the kernel dimensions are as big as the input matrix? – Royi Aug 27 '22 at 18:20
  • It is better if you describe the values as an equation in the question. Also, usually Kernel are much smaller than the image, hence they can be sparse in their matrix representation. In this case they are not, it means that it will be a challenge for a large input. – Royi Aug 27 '22 at 21:18
  • sorry I made a mistake on the indices of the kernel matrix in my last reply. As the index of array in Python starts from 0, the central pixel in kernel $K$ is actually at $(N-1, N-1)$, now I added the correct numbers in the question. – learner_lyf Aug 27 '22 at 21:52
  • In the sense of image transforming, the kernel is usually much smaller than the image. But my question comes from some other application of the convolution method, so it has to be exactly that big. I've tried the original matrix of $2000 \times 2000$, the convolution is still very fast. But I don't know whether the size of the matrix will be an issue when solving my question. I hope it does not matter so much. – learner_lyf Aug 27 '22 at 22:01
  • The problem is the kernel, in the genera case of size $ 2 N $, has $ 4 {N}^{2} $ parameters while the output of the convolution has less. So the system is not only ill poised, it is under determined. – Royi Aug 28 '22 at 05:29
  • I am not sure what you mean. But it is better to an approximated smaller kernel. – Royi Aug 29 '22 at 18:41
  • 1
    @Royi, thanks for your help, a friend explained a bit what your matlab code is doing, I now understand your answer in another question. I just implemented it in Python and it works quite very when the kernel is smaller than the convolution matrix. For this question, I still see the potential of solving it in a similar way as you provided. I will continue exploring it. – learner_lyf Sep 01 '22 at 22:47

1 Answers1

1

The issue is with the kernel being too big and creating a case that not only we're trying to solve ill poised problem, we also have more parameters to estimates than measurements.

One way to handle this would be a lower rank approximation for the kernel. By using the approach in Estimating Convolution Kernel from Input and Output Images one can chose a big support for the kernel which still have lower number of parameters than measurements. Yet the approximation is $ {L}_{2} $ based.

One might try other regularizations on the kernel which can be solved using iterative methods such as ADMM or Accelerated Proximal Gradient Descent.

One could even dare to try solving the undetermined system but probably the results won't be satisfactory.

Royi
  • 19,608
  • 4
  • 197
  • 238