3D multi-plane phase retrieval method based on compressed support estimation, medium and device

Abstract

The application discloses a 3D multi-plane phase retrieval algorithm based on compressed support estimation, a medium and equipment, and comprises the following steps: firstly, a compressed support estimation (CSD) algorithm is used to adaptively estimate the region support of each plane from a single diffraction intensity, and prior knowledge is not required to be predefined; then, a CSD-MIPR-HIO algorithm introduces a HIO cycle as a restriction condition, and the performance of phase retrieval is improved by adjusting the constraint condition of the Pth plane; in the iteration process, the transmission function of each plane is respectively subjected to a modulus constraint in a Fourier domain and a support constraint of a target region, and finally converges to the real distribution of an object quickly. The method disclosed by the application breaks the limitation of prior conditions of a traditional 3D multi-plane phase retrieval, and shows high robustness in a noise environment. Numerical and optical experimental results show that the CSD-MIPR-HIO method is feasible, superior and robust to noise.

Description

Technical Field

[0001] This invention relates to the field of 3D multiplane phase retrieval technology, and more specifically to a 3D multiplane phase retrieval method based on compression support estimation. Background Technology

[0002] With the rapid development of optical imaging applications, more and more 3D information devices are entering the public eye. For 3D objects, their phase carries important information such as the object's shape, internal structure, and propagation depth. Therefore, 3D phase reconstruction has become a key research step in many fields. Currently, 3D phase retrieval technology has been widely applied in various optical imaging fields, significantly impacting the improvement of application technology levels across industries. Among these applications, X-ray crystal structure analysis, used to generate 3D quantitative density maps of molecules, has made significant contributions to the field of routine medical imaging. Improving imaging resolution has always been a driving force for 3D phase retrieval. However, multiple measurements from multiple angles result in high experimental costs and radiation damage to samples. Furthermore, the support requirements for standardizing 3D phase retrieval remain a challenging issue.

[0003] The most widely used iterative methods for recovering phase from measured intensity are the Gerchberg-Saxton (GS) and hybrid input-output (HIO) algorithms. Modulus constraints in the Fourier domain and support constraints on the target field are necessary conditions for convergence of the iterative algorithm. The support mask refers to the smallest closed region outside the target that is zero. When the phase is unknown, satisfying the support constraints on the target field ensures the consistency and uniqueness of the recovered phase. Prior knowledge not only helps to standardize the reconstruction process but also compensates for the loss of other information. However, the accuracy of the prior knowledge assumptions has a significant impact on the reconstruction results; generally, stricter (more detailed) support constraints lead to faster convergence.

[0004] In the field of optical imaging, 3D phase retrieval technology has a wide range of scientific and industrial applications. Traditionally, the most widely used approach for determining 3D structures involves acquiring multiple measurements from different sample orientations, such as tomography and microscopy. Traditional tomography algorithms for 3D phase retrieval typically require calculating the 2D phase at each angle and then inputting the data into the tomography algorithm, thus necessitating multiple measurements at different angles. In 2015, Tian et al. proposed an improvement to the 3D tomography algorithm, enabling 3D sample reconstruction by capturing only a single intensity image at each angle. However, multi-angle scanning places stringent demands on experimental equipment. Furthermore, for samples with sensitive characteristics, each measurement causes radiation damage, altering the sample structure at the atomic scale. High-resolution reconstruction of thick objects has always been challenging due to the potential for multiple scattering within the sample. In recent years, lensless imaging has made some progress in phase retrieval. To recover 3D objects using coherent diffraction microscopy, multiple diffraction patterns are required in different sample orientations. Precise mechanical tilting and high experimental costs are major limiting factors. Typically, traditional techniques require multiple, multi-angle measurements. However, due to limitations in experimental setups and the radiation-induced sensitivity of samples, multiple scans or rotations of the samples are difficult to achieve. Therefore, recovering object information from a single intensity measurement is an important research direction, enabling high-speed, high-resolution acquisition and providing significant reference value for real-time imaging.

[0005] In 2021, Latychevskaia proposed a Multislice Iterative Phase Retrieval (MIPR) method. This algorithm proposes that when an object is probed with coherent waves, a single 2D intensity measurement contains all the information about the 3D sample distribution. The MIPR method uses alternating reconstructions of the transmission function and the incident wavefront during both forward and backward propagation to recover the 3D object distribution from a single 2D intensity measurement in the far field. This algorithm overcomes the limitation of requiring multiple measurements for 3D objects, achieving single reconstruction from far-field measurements. However, far-field intensity is often the only measurement in phase retrieval problems. While this technique is effective, it is limited by prior knowledge of the region's support, reducing the flexibility of reconstruction and exhibiting significant limitations. Furthermore, due to the sensitivity of the GS algorithm to noise, MIPR demonstrates poor robustness. Summary of the Invention

[0006] The present invention proposes a 3D multiplane phase recovery method based on compression support estimation, which can at least solve one of the above-mentioned technical problems.

[0007] To achieve the above objectives, the present invention adopts the following technical solution:

[0008] A 3D multiplane phase retrieval method based on compressive support estimation includes:

[0009] Includes the following steps,

[0010] First, the region support of each plane is adaptively estimated from the single diffraction intensity using the Compressed Support Estimation (CSD) algorithm, without the need for predefined prior knowledge;

[0011] Then, the CSD-MIPR-HIO algorithm introduces the HIO loop as a constraint condition to improve the phase recovery performance by adjusting the constraint condition of the Pth plane. During the iteration process, the transmission function of each plane is subject to the modulus constraint of the Fourier domain and the support constraint of the target region, and finally converges quickly to the true distribution of the object.

[0012] Furthermore, this includes the following steps:

[0013] S1: Mask estimation: The single 2D diffraction intensity is used as the input of the CSD module to adaptively estimate the support mask of each plane based on the diffraction information;

[0014] S2: Initial Reconstruction: Perform a central inverse Fourier transform on the input single-frame 2D diffraction intensity to obtain the transmission function of the P-th plane. The angular spectrum method (ASM) is used to simulate the propagation process between planes; the obtained... Backpropagation is performed, and support constraints are applied to filter the reconstructed 2D distributions of each plane. At the end of propagation, a set of initial transmission functions is obtained. , , … ;

[0015] S3: Forward Propagation: Initial Transmission Function , , … Forward propagation, where the transmitted wavefront of the first plane Initialize to And use ASM simulation Propagation to obtain the second plane ;

[0016] Similarly, forward propagation yields... ,right Diffraction intensity is obtained using the central Fourier transform. And updated to using modulus constraints ;

[0017] S4: Backpropagation: Updated diffraction information Backpropagation, for the obtained HIO algorithm constraints are as follows And propagate backwards to update the transmission function of each plane, where , , … This is the reconstruction result of this iteration.

[0018] S5: Iterative Loop: Repeat S3 and S4 continuously until the algorithm converges. The final reconstruction result of each plane is then obtained. , , … .

[0019] In another aspect, the present invention also discloses a computer-readable storage medium storing a computer program, which, when executed by a processor, causes the processor to perform the steps of the method described above.

[0020] In another aspect, the present invention also discloses a computer device, including a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor performs the steps of the method described above.

[0021] As can be seen from the above technical solutions, traditional 3D multi-plane phase retrieval algorithms face the dual limitations of predefined supports and iterative stagnation. To reduce the need for known priors, this invention proposes a 3D multi-plane phase retrieval algorithm based on compressed support estimation (CSD-MIPR-HIO). This algorithm uses compressed sensing to adaptively estimate the supports of each plane, requiring no prior knowledge of the region, and recovers the 3D object only from a single 2D diffraction intensity. Simultaneously, CSD-MIPR-HIO employs a hybrid input-output (HIO) iterative algorithm, improving upon the GS loop's tendency to stagnate in traditional multi-plane algorithms, thereby enhancing the algorithm's convergence. The method proposed in this invention overcomes the limitations of traditional 3D multi-plane phase retrieval based on prior conditions and exhibits high robustness in noisy environments. Numerical and optical experimental results demonstrate the feasibility, superiority, and noise-resistant robustness of the CSD-MIPR-HIO method. Attached Figure Description

[0022] Figure 1 It is a 3D multi-planar object coherent diffraction imaging optical path diagram;

[0023] Figure 2 This is a schematic diagram of the MIPR algorithm;

[0024] Figure 3 This is a schematic diagram of the compression support estimation principle;

[0025] Figure 4This is a flowchart of the CSD-MIPR-HIO method according to an embodiment of the present invention;

[0026] Figure 5 This is a comparison diagram between the MIPR method and the CSD-MIPR-HIO method of this invention embodiment;

[0027] Figure 6 This is a schematic diagram showing two different 3D objects selected as test targets in the experiment.

[0028] Figure 7 This is a CSD-MIPR-HIO reconstruction diagram of the invention embodiment. Detailed Implementation

[0029] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are some embodiments of the present invention, but not all embodiments.

[0030] To overcome the dual limitations of predefined support and iterative stagnation in the Multi-Slice Iterative Phase Retrieval (MIPR) method, this invention proposes a 3D multi-plane phase retrieval method based on compressed support estimation (CSD-MIPR-HIO). This method first uses the Compressed Support Detection (CSD) algorithm to adaptively estimate the regional support of each plane from a single diffraction intensity, requiring no predefined prior knowledge. Then, CSD-MIPR-HIO introduces an HIO loop to impose dual constraints on the modulus (squared modulus of the Fourier transform) and the region during the algorithm's iteration. The HIO loop escapes local minima, solving the problem of slow convergence and low robustness of the MIPR method due to GS iteration. The proposed method starts with a single intensity measurement in the far field, adaptively estimating regional support, and accurately reconstructing 3D multi-plane objects in noisy environments. Compared to the MIPR algorithm, CSD-MIPR-HIO significantly improves the flexibility and robustness of 3D multiplane phase retrieval.

[0031] The following are explanations:

[0032] First, let's introduce the basic principles of traditional MIPR:

[0033] Imaging optical path introduction:

[0034] The imaging principle of 3D phase retrieval is as follows: Figure 1 As shown, a 3D multiplane object is illuminated with a light beam, and diffraction intensity is measured on the detector plane in the far field. The 3D multiplane phase retrieval algorithm (MIPR) is used to accurately reconstruct the real spatial image of the 3D multiplane object from a single 2D diffraction intensity.

[0035] A 3D multi-plane object is composed of P 2D planes, and the object characteristics of each plane are determined by the transmission function. This indicates that the incident wave propagates through P planes and is modified after passing through each plane. The wavefront, after propagating through each plane, can be divided into the incident wavefront and the wavefront itself. and transmitted wavefront The propagation calculation formula is as follows:

[0036] (1)

[0037] The propagation process between planes is simulated using the angular spectrum method (ASM), and the propagation formula is:

[0038] (2)

[0039] in, Represents the wavefront distribution on plane p, and is composed of It was obtained using the angle spectrum method. Represents Fourier plane coordinates. and These are the Fourier transform and the inverse Fourier transform, respectively. For wavelength, It represents the distance between two planes.

[0040] Introduction to MIPR Iteration Principle

[0041] The principle of the MIPR algorithm is as follows: Figure 2 As shown. The known information for the MIPR algorithm includes the single-image 2D diffraction intensity distribution and the supporting masks for each plane. The algorithm's iterative process consists of four steps:

[0042] 1) Initial Reconstruction: Initial reconstruction is performed using a single intensity diffraction pattern and a predefined support mask to obtain the initial features of each plane. ;

[0043] 2) Forward propagation: propagating the first plane... Initialize to And by simulating angular spectrum propagation (Formula 2), we obtain... From initial features Forward propagation yielded the following results: and , will the P-th plane The diffraction intensity is obtained by central Fourier transform and amplitude constraint. ;

[0044] 3) Backpropagation: by Perform a central inverse Fourier transform and update the features of each plane through angular spectrum propagation. ;

[0045] 4) By iterating through steps 2) and 3), the final reconstruction results of each plane are obtained. .

[0046] The MIPR algorithm can achieve high-precision reconstruction of multi-plane 3D objects from a single diffraction pattern. However, this method relies heavily on prior knowledge—pre-defined support masks for each plane of the 3D object; furthermore, because the MIPR algorithm uses a multi-plane GS loop, the iterative process is prone to stagnation, resulting in poor convergence.

[0047] The following is the method CSD-MIPR-HIO according to an embodiment of the present invention:

[0048] Compression support estimation

[0049] Forward measurement model:

[0050] Under coherent light illumination conditions Figure 1 The imaging process can be modeled as follows:

[0051] (3)

[0052] in The holographic intensity recorded by the detector. It includes a DC term caused by the reference light, and a vectorized term. ,here It is the scattering field of the surface of a 3D object. This represents additional noise. ,in This represents a two-dimensional inverse DFT matrix. , Represents a block diagonal matrix. This represents the 2D Discrete Fourier Transform, and its size is... ,here and These represent the number of detector pixels in the x and y directions, respectively. ,and

[0053] (4)

[0054] Its representation matrix No. Line number The element values ​​of the column, It is the wavenumber of light. This represents the sampling interval on the z-axis. This represents the sampling interval in the k-space.

[0055] Supporting estimation algorithm:

[0056] Equation (1) can be transformed into the following problem by decompressing the interference and applying a TV (Total Variation) constraint.

[0057] (5)

[0058] (6)

[0059] in Defined as

[0060] (7)

[0061] in, This represents a 3D multi-plane sample. This represents a corrosion operation with a corrosion coefficient of W. A mask corresponding to an approximate 3D sample. This represents the final estimated mask. Finally, a two-step iterative shrinkage / thresholding (TwIST) algorithm is used to solve this optimization problem.

[0062] To overcome the limitation of MIPR requiring a known object support mask, this invention introduces a compressed support estimation algorithm to adaptively estimate the support mask for different planes of a 3D object. For example... Figure 3 As shown, for a given 3D multi-plane object, the method of this invention adaptively estimates the support mask of different planes of the 3D object from a single 2D diffraction intensity. Specifically, it consists of a forward measurement model and a support estimation algorithm: the forward measurement model obtains a single 2D diffraction intensity map by propagating from each sample plane to the detector plane. The support estimation algorithm involves backpropagating from the detector plane to each sample plane to ultimately obtain multiple sets of data. The specific steps supporting the estimation algorithm are as follows:

[0063] Algorithm Support Estimation Algorithm

[0064] Input: A captured single 2D diffraction intensity map

[0065] Output: Estimated 3D sample mask

[0066] The TwIST algorithm is used to reconstruct formula (5) to obtain an approximate 3D sample. .

[0067] By setting the corrosion coefficient in formula (6), multiple sets of [various methods] can be obtained. This serves as the basis for the final estimation of the 3D samples. End

[0068] For reconstructing 3D multiplanar objects from a single 2D diffraction intensity map, accurately constraining the supports of each plane during the iterative process is crucial. The CSD method adaptively estimates the masks of each 3D plane based on diffraction intensity and obtains different support masks by setting different erosion coefficients W. Its accurate estimation of the support masks can improve the convergence of the iterative algorithm. The CSD method provides a simple way to adaptively estimate the supports of all planes, overcoming the limitation of the original MIPR algorithm that relies on prior knowledge of the support masks.

[0069] Multiplane iteration

[0070] To address the issues of local optima and noise sensitivity inherent in the GS loop of the MIPR algorithm, this invention proposes a 3D multi-plane phase retrieval method based on compressed support estimation (CSD-MIPR-HIO). This method reconstructs the 3D object distribution from a single 2D intensity measurement without predefined support masks. Simultaneously, the CSD-MIPR-HIO algorithm introduces an HIO loop as a constraint, improving phase retrieval performance by adjusting the constraints of the P-th plane. During iteration, the transmission function of each plane is subject to modulus constraints in the Fourier domain and support constraints in the target region, ultimately converging rapidly to the true object distribution. CSD-MIPR-HIO improves the reconstruction quality of the traditional MIPR algorithm, overcoming the drawbacks of stagnation and slow convergence of the traditional GS iterative algorithm. It also overcomes the limitation of requiring predefined target supports in 3D multi-plane phase retrieval, while enhancing the algorithm's robustness. The principle of the CSD-MIPR-HIO algorithm is as follows: Figure 4 As shown.

[0071] The CSD-MIPR-HIO algorithm uses a single 2D diffraction intensity map as input, does not require a known support mask, and consists of four steps:

[0072] S1: Mask estimation: The single 2D diffraction intensity is used as the input of the CSD module to adaptively estimate the support mask of each plane based on the diffraction information.

[0073] S2: Initial Reconstruction: Perform a central inverse Fourier transform on the input single-frame 2D diffraction intensity to obtain the transmission function of the P-th plane. The angular spectrum method (ASM) is used to simulate the propagation process between planes. The obtained... Backpropagation is performed, and support constraints are applied to filter the reconstructed 2D distributions of each plane. At the end of propagation, a set of initial transmission functions is obtained. , , … .

[0074] S3: Forward Propagation: Initial Transmission Function , , … Forward propagation, where the transmitted wavefront of the first plane Initialize to And use ASM simulation Propagation to obtain the second plane Similarly, forward propagation yields... ,right Diffraction intensity is obtained using the central Fourier transform. And updated to using modulus constraints .

[0075] S4: Backpropagation: Updated diffraction information Backpropagation, for the obtained HIO algorithm constraints are as follows And propagate backwards to update the transmission function of each plane, where , , … This is the reconstruction result of this iteration.

[0076] S5: Iterative Loop: Repeat S3 and S4 continuously until the algorithm converges. The final reconstruction result of each plane is then obtained. , , … .

[0077] Compared to MIPR, the CSD-MIPR-HIO algorithm differs in the following ways: 1. The CSD-MIPR-HIO algorithm requires only a single 2D diffraction intensity map and uses compressed sensing to adaptively estimate the multi-plane mask. 2. The CSD-MIPR-HIO algorithm employs the HIO iterative algorithm, improving the robustness of the 3D multi-plane phase retrieval algorithm to noise. A comparison of the MIPR and CSD-MIPR-HIO methods is shown in the figure below. Figure 5 As shown.

[0078] To verify the effectiveness and superiority of the CSD-MIPR-HIO method proposed in this invention, a feasibility experiment of the CSD-MIPR-HIO method was designed in this part. Phase retrievals were performed on two different groups of 3D multi-planar objects respectively, verifying that the method of this invention can effectively reconstruct the phases of 3D multi-planar objects.

[0079] Experiment A: Feasibility Experiment

[0080] To verify the effectiveness of the CSD-MIPR-HIO method for reconstructing the phases of 3D multi-planar objects, two different 3D objects were selected as test targets in this group of experiments. The 3D objects are as Figure 6 shown. This group is a numerical simulation experiment based on MatlabR2021b, in which the beam wavelength is set to 532 nm.

[0081] Figure 6 . Schematic diagram of 3D objects (a) The 3D object is composed of four stacked 400×400 two-dimensional planes. The main information of the four planes is ’α’, ’β’, ’γ’, ’δ’. The interval Δz between adjacent two planes is 50 μm. (b) The 3D object is composed of four stacked 400×400 two-dimensional planes. The main information of the four planes is “An”, “Hui”, “Da”, “Xue”. The interval Δz between adjacent two planes is 50 μm.

[0082] As Figure 7 shown, two groups of 3D objects in Figure 6 were selected as test objects in this group of experiments. First, the support masks of each layer were estimated by the CSD method. Three different masks were obtained by setting the erosion coefficient W to 1, 3, and 5 respectively. Different masks were selected for support during the HIO cycle for reconstruction. The experimental iteration number was set to 1000, and the parameter in the HIO algorithm was set to 0.9. The reconstruction results of the 3D objects under different support masks are as Figure 7 shown, where the performance evaluation indexes are CC value, PSNR and SSIM.

[0083] As Figure 7 shown, the reconstruction results of the first group of objects and the second group of objects both show that the reconstruction effects under different mask constraints are different, and the performance of the support mask constraint reconstruction with the erosion coefficient W set to 1 is the best. Under the masks with the same erosion coefficient, since the information of each sample plane of the second group of objects is more complex, the reconstruction effect of the first group of objects is better than that of the second group of objects. The numerical values of each reconstruction index all show that the CSD-MIPR-HIO method can accurately reconstruct 3D objects.

[0084] In another aspect, the present invention also discloses a computer-readable storage medium storing a computer program, which, when executed by a processor, causes the processor to perform the steps of any of the methods described above.

[0085] In another aspect, the present invention also discloses a computer device, including a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor performs the steps of any of the methods described above.

[0086] In another embodiment provided in this application, a computer program product containing instructions is also provided, which, when run on a computer, causes the computer to perform the steps of any of the methods described in the above embodiments.

[0087] It is understood that the system provided in the embodiments of the present invention corresponds to the method provided in the embodiments of the present invention, and the explanation, examples and beneficial effects of the relevant content can be referred to the corresponding parts of the above methods.

[0088] Those skilled in the art will understand that all or part of the processes in the methods of the above embodiments can be implemented by a computer program instructing related hardware. The program can be stored in a non-volatile computer-readable storage medium, and when executed, it can include the processes of the embodiments of the above methods. Any references to memory, storage, databases, or other media used in the embodiments provided in this application can include non-volatile and / or volatile memory. Non-volatile memory can include read-only memory (ROM), programmable ROM (PROM), electrically programmable ROM (EPROM), electrically erasable programmable ROM (EEPROM), or flash memory. Volatile memory can include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM is available in various forms, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous DRAM (SDRAM), dual data rate SDRAM (DDRSDRAM), enhanced SDRAM (ESDRAM), synchronous link DRAM (SLDRAM), Rambus direct RAM (RDRAM), direct memory bus dynamic RAM (DRDRAM), and memory bus dynamic RAM (RDRAM), etc.

[0089] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0090] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A 3D multiplane phase retrieval method based on compression support estimation, characterized in that, Includes the following steps, First, the region support of each plane is adaptively estimated from the single diffraction intensity using the Compressible Support Estimation (CSD) algorithm, without the need for predefined prior knowledge. Then, the CSD-MIPR-HIO algorithm introduces the HIO loop as a constraint condition to improve the phase recovery performance by adjusting the constraint condition of the Pth plane. During the iteration process, the transmission function of each plane is subject to the modulus constraint of the Fourier domain and the support constraint of the target region, and finally converges quickly to the true distribution of the object. The CSD-MIPR-HIO algorithm includes the following steps: S1, Mask estimation: The single 2D diffraction intensity is used as the input of the CSD module, and the support mask of each plane is adaptively estimated based on the diffraction information; S2. Initial Reconstruction: Perform a central inverse Fourier transform on the input single-frame 2D diffraction intensity to obtain the transmission function of the P-th plane. The angular spectrum method (ASM) is used to simulate the propagation process between planes; the obtained... Backpropagation is performed, and support constraints are applied to filter the reconstructed 2D distributions of each plane. At the end of propagation, a set of initial transmission functions is obtained. , , … ; S3, Forward Propagation: Initial Transmission Function , , … Forward propagation, where the transmitted wavefront of the first plane Initialize to And use ASM simulation Propagation to obtain the second plane ; Similarly, forward propagation yields... ,right Diffraction intensity is obtained using the central Fourier transform. And updated to using modulus constraints ; S4, Backpropagation: Updating diffraction information Backpropagation, for the obtained Apply HIO algorithm constraints as follows And propagate backwards to update the transmission function of each plane, where , , … This is the reconstruction result of this iteration; S5. Iterative Loop: Repeat S3 and S4 continuously until the algorithm converges. The final reconstruction result of each plane is then obtained. , , … .

2. A computer-readable storage medium storing a computer program, which, when executed by a processor, causes the processor to perform the steps of the 3D multiplane phase retrieval method as described in claim 1.

3. A computer device comprising a memory and a processor, the memory storing a computer program that, when executed by the processor, causes the processor to perform the steps of the 3D multiplane phase retrieval method as claimed in claim 1.

Images