# Photoacoustic Imaging with a Finite-Size Circular Integrating Detector
## 3 Algorithms
### Linear Equation Formulation of the Blurring Kernel Extraction Problem  
This section elaborates on the blurring kernel estimation method based on polynomial algebraic common factor extraction. Initially, we transform Equation (8) into the frequency domain to simplify the convolution operation. Consequently, the relationship between the blurred and initial pressure distributions in the frequency domain can be expressed as:  
\[P_{j}^{\theta}(\omega)=P_{j}^{0}(\omega)\cdot H^{\theta}(\omega) \tag{9}\]where \(P_{j}^{\theta}(\omega)\) represents the blurred pressure distribution \(p_{j}^{\theta}(\phi)\) at a certain radial position \(r_{j}\) in the frequency domain, \(P_{j}^{\theta}(\omega)\) represents the initial pressure distribution \(p_{j}^{0}(\phi)\) in the frequency domain, and \(H^{\theta}(\omega)\) represents the blurring kernel \(h_{\theta}(\phi)\) in the frequency domain. Since the blurring kernel \(H^{\theta}(\omega)\) of the circular detector is independent of radial positions, the blurring kernel is consistent across all radial positions. Thus, we can derive the following relationship:  
\[H^{\theta}(\omega)=\frac{P_{1}^{\theta}(\omega)}{P_{1}^{0}(\omega)}=\frac{P_{ 2}^{\theta}(\omega)}{P_{2}^{0}(\omega)}=\cdot\cdot\cdot=\frac{P_{j}^{\theta}( \omega)}{P_{j}^{0}(\omega)} \tag{10}\]  
where \(j=1,2,\cdot\cdot\cdot,J\) indexes the radial positions in the polar coordinate image. Each position corresponds to an independent blurring problem. Taking two radial positions \(r_{1}\) and \(r_{2}\) as an example, we have the following proportionality relationship:  
\[H^{\theta}(\omega)=\frac{P_{1}^{\theta}(\omega)}{P_{1}^{0}(\omega)}=\frac{P_{ 2}^{\theta}(\omega)}{P_{2}^{0}(\omega)} \tag{11}\]  
We can express this as:  
\[P_{1}^{\theta}(\omega)\cdot P_{2}^{0}(\omega)-P_{2}^{\theta}(\omega)\cdot P_{ 1}^{0}(\omega)=0 \tag{12}\]  
The frequency domain multiplication can be transformed into the time domain convolution form through the inverse Fourier transform:  
\[P_{1}^{\theta}(\phi)*P_{2}^{0}(\phi)-P_{2}^{\theta}(\phi)*P_{1}^{0}(\phi)=0 \tag{13}\]  
where \(*\) denotes the convolution operation. We can discretize Equation (13):  
\[P_{1}^{\theta}(k)*P_{2}^{0}(k)-P_{2}^{\theta}(k)*P_{1}^{0}(k)=0 \tag{14}\]  
Convolutions in the time domain can be represented as matrix operations. We derive the following matrix form:  
\[\mathbf{P_{1}^{0}P_{2}^{0}-P_{2}^{0}P_{1}^{0}=0} \tag{15}\]  
where \(\mathbf{P_{i}^{0}}\) is the discretized vector of the initial pressure distribution \(p_{i}^{0}(\phi)\) at radial position \(r_{i}\). The matrix \(\mathbf{P_{j}^{0}}\) is a special Toeplitz matrix, constructed as follows: the first column consists of the discrete values of the blurred pressure distribution \([p_{j}^{\theta}(\phi_{1}),p_{j}^{\theta}(\phi_{2}),\cdot\cdot\cdot,p_{j}^{ \theta}(\phi_{M})]\), followed by \(K-1\) zeros. Each subsequent column is a cyclic shift of the first column:  
\[\mathbf{P_{j}^{0}}(0)= \begin{pmatrix}p_{j}^{\theta}(\phi_{1})&0&\cdot\cdot \cdot&0&0\\ p_{j}^{\theta}(\phi_{2})&p_{j}^{\theta}(\phi_{1})&\cdot\cdot\cdot&0&0\\ \cdot\cdot&\cdot\cdot&\cdot\cdot&\cdot\cdot&\cdot\cdot\\ 0&0&\cdot\cdot&p_{j}^{\theta}(\phi_{M})&p_{j}^{\theta}(\phi_{M})\\ 0&0&\cdot\cdot\cdot&0&p_{j}^{\theta}(\phi_{M})\\ \end{pmatrix}  \tag{16}\]  
The matrix dimensions satisfy \(K=M-L+1\), where \(M\) is the discrete length of the blurred pressure distribution, \(K\) is the discrete length of the initial pressure distribution, and \(L\) is the discrete length of the blurring kernel.