An Analytical Approach for Fast Recovery of the LSI Properties in Magnetic Particle Imaging

Linearity and shift invariance (LSI) characteristics of magnetic particle imaging (MPI) are important properties for quantitative medical diagnosis applications. The MPI image equations have been theoretically shown to exhibit LSI; however, in practice, the necessary filtering action removes the first harmonic information, which destroys the LSI characteristics. This lost information can be constant in the x-space reconstruction method. Available recovery algorithms, which are based on signal matching of multiple partial field of views (pFOVs), require much processing time and a priori information at the start of imaging. In this paper, a fast analytical recovery algorithm is proposed to restore the LSI properties of the x-space MPI images, representable as an image of discrete concentrations of magnetic material. The method utilizes the one-dimensional (1D) x-space imaging kernel and properties of the image and lost image equations. The approach does not require overlapping of pFOVs, and its complexity depends only on a small-sized system of linear equations; therefore, it can reduce the processing time. Moreover, the algorithm only needs a priori information which can be obtained at one imaging process. Considering different particle distributions, several simulations are conducted, and results of 1D and 2D imaging demonstrate the effectiveness of the proposed approach.


Introduction
Magnetic particle imaging (MPI) is a new method for imaging the spatial distribution of magnetic nanoparticles, as tracers, with high resolution. The method was proposed by Gleich and Weizenecker [1] and exploits the nonlinear magnetization response of the nanoparticles to a timevariable magnetic field and allows for fast image acquisition. MPI has many applications, especially in medical diagnosis such as blood flow visualization for coronary artery diseases, cancer detection [2,3], stem cell tracking [4], and molecular imaging [5].
MPI uses an oscillating drive field (excitation field) of sufficient amplitude to change the magnetization of the nanoparticles, which induces a voltage signal in the receive coils. To enable spatial encoding of the information, a static magnetic gradient field, also known as the selection field, is utilized in MPI. This field contains a spatial location named a field-free point (FFP) that has zero field magnitude, and only the particles located at the FFP induce the MPI signal in the receive coils [6][7][8].
Image reconstruction in MPI includes two main approaches [9]. The first approach makes use of the system matrix. This method, also known as frequency space reconstruction, employs a large system matrix and requires its inversion and some postprocessing to deal with poor conditioning. The approach has high computational load and requires the system matrix to be estimated before imaging [10][11][12]. The second approach is the -space method, which was introduced by Goodwill et al. [13][14][15]. This approach does not require matrix inversion or precharacterization and hence provides a robust reconstruction algorithm with a potential for real-time imaging.
Linearity and shift invariance (LSI) are important characteristics of most medical diagnostic systems. For example, since X-ray computed tomography (CT) is LSI, the CT images 2 International Journal of Biomedical Imaging of tissue attenuation coefficient maps provide quantitative lumen diameters for cardiovascular diagnosis [16]. It has been theoretically proven that the MPI is an LSI system [17]. However, in practice, to detect MPI signal, it is necessary to utilize a filter to remove the drive field signal at its fundamental excitation frequency, which is an unavoidable phenomenon in current MPI techniques. This filtering action also removes the induced particle signal at the excitation frequency and hence the complete MPI signal is not available in practice which destroys the LSI properties of the MPI image.
To recover the complete MPI signal to facilitate the use of the fast -space imaging method in qualitative medical diagnosis, an algorithm is presented in [16] and is optimized in [18]. The approaches use the partial field of views (pFOVs), which are primarily generated for enlarging the FOV. It is shown that the lost information can be recovered by matching the two successive overlapped pFOVs, provided the particle concentration at one reference location is known. Although the methods provide satisfactory results, the computational load is reported in [19] to be relatively high, which increases the imaging time and so this cannot be used for ultrafast imaging.
In this paper, the effect of the filtered MPI image is analyzed and it is shown that its value is constant provided that the trajectory of the FFP is generated by a harmonic function, a property that has to be satisfied in both one-dimensional (1D) and multidimensional imaging. Assuming a constant image loss, a fast analytical algorithm is proposed to restore the lost information. The approach can be real time since its complexity depends on solution of a small-sized system of linear equations and does not require overlapping of pFOVs. The method utilizes the derivative of the magnetization curve of the particles and a mathematical model of the -space image. It also requires a priori information regarding the system, which can be obtained experimentally or estimated theoretically prior to imaging. This is another advantage of the proposed approach in comparison to previous works, which require a boundary condition to be obtained at the start of imaging. Several 1D and 2D imaging simulations are presented to demonstrate the effectiveness of the proposed algorithm.

Methods
The developed recovery algorithm is based on the mathematical model of the MPI signal. In this section, the signal model of the -space method is presented and a model of the term removed from the signal as a result of filtering the fundamental frequency of the excitation signal is derived. Then, a new algorithm is proposed to recover the lost image.

Mathematical Model of a 1D MPI Signal.
In this section, the mathematical model of the 1D MPI signal is derived which is used in Section 2.2 to model the -space reconstruction method and to analyze the LSI properties of the -space.
When the excitation field in MPI is periodic, the voltage signal V( ), induced by the particles' magnetization, is periodic as well, and therefore it can be expanded into a Fourier series as follows: where are the Fourier coefficients and 0 = 2 0 with 0 as the frequency of the excitation field. The coefficients can be written as where = 1/ 0 . Assuming a 1D distribution of the particles in thedirection and a constant sensitivity for the receive coil in this direction, the induced voltage by the time-varying magnetization of the particles can be written as follows [17]: where 0 denotes the permeability of free space and is the particle magnetization, which depends on the magnetic field ( , ). The particle magnetization can be modeled using the Langevin function L(⋅) as follows: where [Am 2 ] is the particle magnetic moment, [mA −1 ] is a property of the magnetic particle, and ( ) [particles/m 3 ] is the particle concentration that has to be measured in an MPI image. According to (3), the time derivative of magnetization is of interest and can be written as On the other hand, the magnetic field that particles experience in MPI is as follows: where is the selection field and is the excitation or drive field. Therefore, using (6), (5) can be written as whereL( ) = L( )/ . Now, using (3) and (7), (2) can be written as follows: International Journal of Biomedical Imaging 3 By defining the function (8) can be written as ( ) are the so-called system function, which in fact are the Fourier coefficients of an induced voltage by a point-like distribution of particles at position . (3) is "object" if the particles are ideal (i.e., have a step-like magnetization cure). When the particle magnetization is modeled by the Langevin function, the field of view is slightly different from the object to be imaged. Therefore, in the remaining equations, the region of integration is replaced by "FOV."

Remark 1. The region of integration in
For a homogeneous harmonic drive field with cosine waveform as and a linear selection field ( ) = , where is the gradient of the selection field in the -direction, it has been shown in [20] that ( ) can be written by a convolution as where ( ) = sin (( + 1) arccos ( )) sin (arccos ( )) is the Chebyshev polynomial of order . Now, substituting (12) in (10), one has where ideal represent the system function for the ideal particles, which are defined as Now, one can rewrite (14) as follows: wherê( From (15) and (16), one can conclude that correspond to coefficients of a Chebyshev series (see, e.g., [21]), and thereforê( ) can be written based on the Chebyshev series as follows:̂( This equation is used in the following subsection to describe the 1D -space reconstruction method.

-Space Imaging and LSI Characteristics.
In the -space reconstruction method, the image is defined to bê( ( )) at position ( ) = FFP ( ), where FFP ( ) denotes the position of the FFP.̂( FFP ( )) can be measured from the induced voltage V( ) through the following relation [13]: wherėF FP ( ) is the time derivative of FFP ( ), that is, the velocity of the FFP. The formulation of̂( ) by a convolution in (17) demonstrates that the -space MPI image is a linear and shift-invariant image, properties which are characteristic of most clinical imaging systems, for example, ultrasound, CT, and magnetic resonance imaging (MRI). However, to acquire an MPI image while retaining its LSI properties, one needs the measurement of the true value of V( ), which is not available in practice. In fact, a band-stop filter is used in practice to segregate V( ) from the signal induced by the drive field (11) in the receive coil. The filter suppresses the first harmonic of the signal at the frequency of the drive field and leaves the signal at higher harmonics. According to (13) and (18), the filtered first harmonic can be written as [16]̂( which is a constant value in the image. It should be noted that since 1 is the Fourier coefficient of an odd real function (as a result of the sinusoidal drive field, V( ) is an odd function), its value is imaginary and therefore, according to (20), the lost image is a real constant value. Using (15) and (16), one can write (20) as follows: which shows that the lost value varies when̂( ) is shifted due to the velocity term √1 − ( / ) 2 , and hence its envelope at different positions resembles the sinusoidal excitation pattern (note that sin(arccos ) = √ 1 − 2 ). This means that the lost value is maximal at the center of the FOV and minimal at its edges. In the next subsections, the effect of this constant loss is analyzed and a compensation algorithm is proposed.

LSI Analysis.
According to the definition, the MPI image is linear when the image intensity is linearly proportional to the concentration of the particles and is shift-invariant when the image intensity is independent of the location of the particles. Considering the expression of the lost image in (21) and the fact that the image loss depends on the location of the particles, it can be verified that the lost information destroys the LSI properties of the MPI image [16].
To visualize the abovementioned phenomenon, consider Figure 1 which shows the -space MPI images of a set of magnetic nanoparticles located at different positions of the FOV with concentration 1 ( ). As demonstrated, there is a difference in the maximum value of the acquired image intensity for the same particles at different locations. Also, it is clear from this figure that the lost constant value is larger when the particles are located at the center of the FOV, since the velocity of the FFP is maximum at the center.

Recovery Algorithm.
Based on the results in the previous section, it can be verified that an -space MPI without first harmonic information does not exhibit LSI and, to provide an image for quantitative medical diagnosis, the lost constant image should be recovered. An algorithm is presented in [16,18], which can recover the lost image by dividing the FOV mechanically or electronically, into several pFOVs having sufficient overlap to reliably match the successive pFOVs. This overlapping requirement increases the imaging time to cover the entire object. The algorithm also requires an image of a location having no particle, for example, a location outside the patient. Although the method fully recovers the lost image, its computational load may be prohibitive for fast and real-time applications.
If the FOV includes a zero-particle area (the reference location for the previous algorithms), finding the lost image is a straightforward task. For example, in MPI images of Figure 1, FOVs have at least one edge at which no particle is located. Therefore, in these cases, finding the global minimum of the image is enough to completely recover the lost image (see Figure 2 for an example). This method, however, fails to completely recover the lost image in general. For example, for the sets of particles shown in Figure 8, only a partial image loss can be compensated by using this method.
International Journal of Biomedical Imaging 5 The proposed method in this paper is a model-based compensation algorithm, which can be applied in general cases. The approach does not require overlapping of successive FOVs to cover the object, and accordingly it can reduce the imaging time. The method utilizes the model of the point spread function (PSF) of 1D -space image (denoted bẏ L(⋅)), which, as noted in [14], can be obtained theoretically with precise matching with experiment. The algorithm also assumes that the particle distribution can be described as the sum of discrete point-like particles with known locations in the FOV. This is achievable through deconvolution of the original image and looking at the local maxima of the signal to determine the location of particles (a common approach to find the local maxima is to find the first-order difference information of the signal; whenever a smooth signal is available, it is easy to use the available basic algorithms to find the particles, which are very robust; in the simulation results, the "findpeaks" command of MATLAB is used to find the particles and their positions). A deconvolution method is studied for MPI signal in [15]. The reported experimental results of MPI signal show an acceptable noise level to roughly deconvolve the signal for the proposed algorithm [13,22]. In addition, the approach requires a priori information, which can be obtained theoretically or from a reference image, having the same property as used in [16] for initialization. However, unlike the previous studies, this information can be obtained before imaging.
Based on the above explanation, the following assumptions, besides the general assumptions of the -space reconstruction [13], turn out to be crucial within the development of the recovery algorithm.   In the following, first, an approach is developed for the case in which the particles in the FOV are not in close proximity; that is, the tails of the image for one set do not affect the image of the other sets. Then, the approach is modified for the more general case.
According to (17), the image for the set of point-like particles can be written aŝ where the subscript denotes the th particle. Also, it is obvious from (21) that the image loss for this set of particles (voxel) is proportional to the particle concentration and an integral term depending onL and the velocity of the FFP. For point-like particles, the derivative of the Langevin function is nonzero only in a small region, and therefore the velocity term can be assumed to be constant in this region (however, for a case in which the PSF is not very narrow, this assumption is valid up to a certain accuracy). As a consequence, the lost image for the set of particles with concentration located at can be approximated aŝ where includes the constant terms in (21) and also the integral of the derivative of the Langevin function (which is the PSF of the image). Equation (23) indicates that the image loss for a set of particles is linearly proportional to the velocity of the FFP at the position of the particles, and hence its envelope, throughout the FOV, resembles a sinusoidal pattern.
Using (22), one can write (23) as in the following: where ≜ /L( ) is a constant value for all positions in the FOV (the value of PSF,L, is maximum and constant at ). This value, which is required in the subsequent development, can be identified theoretically or experimentally through the imaging of a simple set of particles, whose image loss can be easily obtained by finding the global minimum of the signal. In other words, a boundary condition as satisfied in Figure 2, similar to the previous pFOV methods [16], is enough to experimentally measure this value. It should also be noted that, unlike the pFOV methods, this measurement is only required to be done once.
Based on the above result, assuming that the -space imaging conditions are satisfied, the image loss for the particles with concentration located at position can be related to the image loss of the particles with concentration at position , as follows: Since imagê( ) is linear, it can be considered as the sum of the images of the point-like particles distributed throughout the imaging axis [20]. The lost imagê( ) lost associated witĥ( ), therefore, can also be considered as the sum of the image losses associated with each set of particles. Then, assuming that̂( ) is the image of the sets of particles located in the FOV, it can be written aŝ Also, for each individual set, one haŝ wherê( ) read is the measured value for the particles noted by index at position . Therefore, it can be verified that, for two particle sets of and , the subtraction of measured values, defined by , is equal to the subtraction of the real image values at and , respectively; that is, International Journal of Biomedical Imaging which can be measured from the filtered image, based on the assumption that the particles can be identified in the FOV. Now, using (24) and (25), one haŝ Therefore, considering an initial particle set, noted here by index 1, the total image loss can be written based on̂1( 1 ) lost as follows: Now, substituting (27) in (30), the image loss can be measured based on available information as follows: Equation (31) can be applied for cases in which the discrete particles are sufficiently distant from each other such that the tails of the image for one set do not affect the images of other particles. However, when the particles are close to each other, the algorithm should be modified. This issue is discussed below for the case of two particles; however, it can be easily extended to cases involving multiple particles.
Consider the particles and to be very close to each other; therefore, for these particles, (27) will be as follows: which is a system of linear equations for the variableŝ( ) and̂( ). The solution of this system gives the values of interaction terms which are a function of̂( ) lost . Finally, these terms should be subtracted from the "read" values in (28) and (31), and the total lost image can be obtained from (31). In the case of multiple particles with multiple interactions, a similar system of equations as in (35) can be developed, where the number of equations is equal to the number of interaction terms. The image loss recovery algorithm is summarized as follows.

Image Loss Recovery for 1D -Space MPI
Require. The model of the 1D PSF: (1) In an initialization scanning, find the image loss for one set of particles by measuring the global minimum of the image signal, and then compute in (24). This step can be omitted when the value of is measured theoretically.
(2) Deconvolve the original image to find the particles and their positions in the FOV through a searching algorithm, and then compute using (23) and ( ) read for each set of particles.
(3) Compute the interaction values of each particle set on the neighbors by solving the system of linear equations, such as (35).
(4) Consider a particle set (indexed 1) as the reference set and measure 1 for all sets using (28), taking into account the effect of the neighbor particles.

Remark 5.
An important note which should be mentioned is that the algorithm requires knowing as narrow convolution kernel (PSF) model (and its image loss to measure ) as possible, which can be used to deconvolve the MPI signal to roughly represent it as the sum of discrete point-like particles. Apparently, when the noise level of the signal allows using a narrower PSF, then the error of approximation in (23) would be less, and the overall performance of the algorithm would be improved.

Remark 6.
From the theoretical point of view, any particle set can be considered as the initial set, indexed by 1 in (31). However, in case that a nonideal filter is implemented in practice, the initial set may affect the accuracy of the algorithm. Selecting this set from those in the center of FOV with a medium concentration can be a proper choice. But, in general, the best choice, which depends on the performance of the filter, can be determined from the results of the first experiments (which can be considered as a calibration procedure).

Remark 7.
Considering the case of a nonideal filter, the same procedure as in Remark 6 can be applied to measure the value of ; that is, an appropriate choice is to measure it for the particle at the center of trajectory.
Remark 8. The approach of [16] uses the average of several scans to find the image loss, and it is guaranteed to be robust with respect to the noise, provided that the image noise has zero mean. Similarly, if this condition is satisfied, the assumption of having smooth signal (introduced in Assumption 3), and accordingly the robustness property, for the proposed approach can be readily achieved by using multiple scans.

Evaluation for Multidimensional MPI.
Before studying multidimensional MPI, it is worth mentioning that the reason that the lost image is constant in 1D MPI is that the Chebyshev polynomial of order zero 0 ( ) is constant and does not depend on the spatial position. On the other hand, thespace image can be written as the Chebyshev series (18) since the system function is expressed by the Chebyshev polynomials. This notation for the system function is the direct effect of a cosine (harmonic) drive field. This means that, for other waveforms of the drive field, the system function cannot be written as the Chebyshev polynomials. Based on the above discussion, it can be verified that, in general, the image loss in multidimensional MPI is not a dc value. However, one method of maintaining a constant image loss in multidimensional -space is to combine the 1D images of the object to generate a 2D/3D image. This method is experimentally studied in [23], which requires only one pair of drive coils and one receive coil.
Using the above method, the image loss would be constant and recoverable for 1D MPI images, but the signal equation of multidimensional MPI differs from that of 1D MPI. Then, a question arises of whether it is possible to use the proposed recovery algorithm for 2D/3D MPI. This question is investigated below by studying the multidimensional MPI signal.
The generalized 3D MPI signal is given by [23] k ( ) = (r) (r) * * * h (r)ṙ FFP r=r FFP ( ) , where r = [ ] ⊤ and h(r) is a matrix PSF which is convolved in 3D (denoted by * * * ) with the particle distribution. When the multidimensional image is reconstructed from 1D images, which is by a collinear arrangement of detection coil and drive coil trajectory, then the collinear PSF ℎ ‖ (r) should be considered. Assuming the drive field trajectory is in the -direction, one has [23] ℎ ‖ (r) =L ( (r)) with (r) = √( ) 2 + ( ) 2 + ( ) 2 . Supposing the particles are located at line = = 0, the second term in (37) will disappear when a drive field trajectory is on this line. The remaining PSF is identical to the 1D PSF, which depends on the time derivative of the Langevin function, and is assumed to be known by the proposed algorithm. When the trajectory moves away from this line, the resultant PSF depends not only on the derivative of the Langevin function, but also on the function itself. Therefore, in different spatial positions for the drive field, the PSF would no longer be the same as the one used by the algorithm. However, the collinear PSF is similar to Gaussian curves, which has a narrow shape at line = = 0 and becomes broader with distance from this line. Figure 3 illustrates a collinear PSF for a point-like sample of particles at different positions of the drive field.
It is well investigated in [22] that it is always possible to reconstruct a broad Gaussian-like curve from linear superposition of a narrow one, which is why multicolor imaging is not possible in 1D MPI of Langevin particles with FFP motion being a single line. This property provides the possibility of describing an MPI signal, obtained from 8 International Journal of Biomedical Imaging a linear FFP trajectory, based on the linear superposition of 1D MPI images resulting from the spatial distribution of the particles along the trajectory line. This distribution can be obtained by deconvolution of the 1D -space PSF from the obtained signal. To clarify this property, an example is illustrated in Figure 4. In this figure, the dashed MPI signal is the output of the band-stop filter, which is induced by the particles at the center of imaging area and the available large set of particles at the border, due to a linear trajectory of the drive field at the center. One should note that although the drive field in Figure 4 does not pass through the large particle set, this set also induces a signal in the receive coil because the FFP is not only a small area for medium-level selection field power, which is generally used in practice.
Using a deconvolution algorithm, the solid line in the figure shows a new representation of the signal based on the 1D PSF, which can be interpreted as the superposition of 1D MPI images of several particles with the same PSF. Therefore, this representation property shows that the proposed recovery algorithm is applicable in multidimensional MPI.
Remark 9. Image of particles with broad PSF, similar to the case where they are close to each other or, in the extreme case, cover the entire FOV, can be expressed as the sum of images of particles with narrow PSF. In fact, the algorithm, regardless of the dimension of particle distribution or imaging, simply tries to find the level of the received signal by expressing the signal as the sum of images of point-like particles (like Figure 4), and then, knowing the effect of a sample set of particles ( in (24)), the algorithm estimates the total image loss.
Remark 10. Nonlinearity of the gradient field or the relaxation effects may cause an asymmetry in the PSF. Considering the point in Remark 9, as long as a model of the PSF is available and the image loss is quasi-constant, the algorithm can still be applied regardless of the asymmetry of the signal. Keep in mind that when the magnetic fields are far from being ideal (ideal magnetic fields: a linear selection field and a homogeneous drive field), the MPI image cannot be formulated as a convolution like (17), and therefore the assumption of constant image loss is not also valid.

Results and Discussion
In this section, simulation results are presented to evaluate the effectiveness of the recovery algorithm proposed in Section 2.4. In the following, simulations are first presented for 1D imaging, and then the result for application of the algorithm to a 2D MPI image is illustrated. MATLAB codes are used to simulate both 1D and 2D images. In 1D simulations, point-like particles are considered in the FOV, and since the algorithm is derived based on the shape/model of the PSF, all images are normalized according to the maximum level of concentration; this means the results are reported regardless of the concentration level of the particles. 1D simulations are conducted in three cases for the distribution of particles, and two different particle sizes are analyzed: particles with a core diameter of 30 nm which have a narrow PSF (considered in the first two cases) and particles with a core diameter of 22 nm (considered in the third case). This is to evaluate the error of the approximation used to derive (23). The PSFs for these particles are illustrated in Figure 5.

Case 1.
As the first study, three sets of particles with different concentrations are considered in the FOV, having no interaction. Since the image of particles does not affect the image of other particles, then the location of particles can be easily determined, and step (3) of image loss recovery for 1Dspace MPI is not required. However, since the particles do not have similar concentrations, must be computed. The result for this case is shown in Figure 6. As the figure demonstrates, in this case, the algorithm can detect the image loss with high accuracy.   Case 2. In the second case, four sets of particles are considered in the FOV, where particles indexed by "2" and "3" are very close to each other and interact. By deconvolving the image, it is possible to easily determine the position of particles. For this case, it is required to use (35) to cancel the interaction effect. The result for this case is shown in Figure 7, which demonstrates the proper compensation of the lost information. It should be noted that if the particles were closer, then deconvolving could not segregate the particles and they would be considered as one set of particles with a higher concentration.
Case 3. In the final case, four sets of particles with a wide PSF are considered in the FOV. In this simulation, the images of particles indexed by "1," "2," and "3" undergo interactions with each other. Also, particles "3" and "4" experience interactions. There are a total of eight interactions in the image; therefore, the equivalent model for (35) includes a system of eight linear equations. To simulate the measurement errors, the level of interactions, obtained by (34), is considered to be 90% of its true value. The result of using the algorithm for this case is illustrated in Figure 8, which shows satisfactory recovery of the image, with a small deviation, despite the wide PSF and error in the measurements. To simulate the algorithm for 2D imaging, the MPI simulation toolbox developed in [24] is utilized. This MATLAB toolbox simulates the magnetic fields and magnetization of the particles and uses the system function method to reconstruct the image. This program is appropriately modified for the -space imaging and the FFP trajectory is changed to meet the requirements for constant image loss in multidimensional imaging. Measurement error, similar to "Case 3," is also considered in this simulation. Figure 9 illustrates the result of applying the recovery algorithm to a 2D image. It is worth mentioning that, in the result, one may expect that Figure 9(c) is better that Figure 9(d), but in fact the concentration of particles shown by Figure 9(d) is more close to the true value. In terms of the image quality, it can be easily improved by combining the images of two orthogonal line-scan drive fields [25].

Conclusion
This paper analyzes the LSI characteristics of the -space MPI image and proposes an algorithm to recover these characteristics, which are important for quantitative medical diagnosis. Although the theoretical model of the -space image is known to exhibit LSI, it is shown that filtering the fundamental excitation frequency, which is necessary in practice to manifest the particle signal, destroys the LSI properties. The lost image could be constant and hence recoverable when a 1D harmonic drive field is applied in both 1D and multidimensional imaging. To recover the lost information, an algorithm is proposed using -space PSF model of the particles and the properties of the mathematical model of the image. The algorithm tries to find the level of the received signal by expressing it as the sum of images of point-like particles, and then, knowing the image loss of a sample set of particles, the algorithm estimates the total image loss. The initial image loss information can be obtained offline and is enough to be measured for one time. In this regard, the proposed approach has an advantage over the previous methods, where a priori information is required to be measured online at the beginning of each test. The complexity of the proposed approach depends on a smallsized system of linear equations, and therefore it has low computational load, which preserves the real-time advantage of -space imaging. This is another advantage of the proposed approach with respect to current available approaches, where overlapping of successive FOVs is required which increases the imaging time to cover the entire object. The results of simulations for 1D and 2D imaging, for both narrow and broad PSFs, demonstrate the effectiveness of the algorithm. The algorithm is developed based on the assumptions expressed in Assumptions 2-4. Although in simulation study these assumptions are released up to some level and the results are still promising, they seem to be the most challenging parts of implementing the algorithm in practice. In this regard, the future work is devoted to the study, through experiments, of the sensitivity of the algorithm to derive a robustness level with respect to the assumptions. The future work is also devoted to the comparison of the processing times of the proposed algorithm with the previous ones.