Login| Sign Up| Help| Contact|

Patent Searching and Data


Title:
METHOD, DEVICE AND STORAGE MEDIUM FOR RECONSTRUCTING A CORRECTED MOSAIC IMAGE FROM BURST IMAGES
Document Type and Number:
WIPO Patent Application WO/2024/013538
Kind Code:
A1
Abstract:
The invention relates to a computer-implemented method for reconstructing a corrected primary mosaic image (MI2) from a set of primary burst images (Ii) of a ground surface area, wherein the method comprises inputting the primary burst images (Ii) and primary metadata, defining a primary mosaic image (Ml), choosing points (Q) in a digital elevation model (DEM) covering the primary mosaic image (Ml), projecting the points (Q) into primary projected points (Q(p)), calculating to which primary burst image (Ii) belongs each primary projected point (Q(p)), applying primary corrections to the primary coordinates (R, COL) of the primary projected points (Q(p)) in the primary mosaic image (Ml), calculating (S6) a primary warp function which transforms the primary coordinates (R, COL) into the primary corrected coordinates (Formula I), transforming each primary burst image (Ii) by the primary warp function (Αi(p)) and assembling them.

Inventors:
AKIKI ROLAND (FR)
ANGER JÉRÉMY (FR)
DE FRANCHIS CARLO (FR)
FACCIOLO GABRIELE (FR)
MOREL JEAN-MICHEL (FR)
Application Number:
PCT/IB2022/000511
Publication Date:
January 18, 2024
Filing Date:
July 15, 2022
Export Citation:
Click for automatic bibliography generation   Help
Assignee:
KAYRROS (FR)
International Classes:
G01S13/90; G06T3/40
Foreign References:
US20210134055A12021-05-06
Other References:
WEGNÜLLER URS ET AL: "Sentinel-1 Support in the GAMMA Software", PROCEDIA COMPUTER SCIENCE, vol. 100, 7 October 2016 (2016-10-07), AMSTERDAM, NL, pages 1305 - 1312, XP093020498, ISSN: 1877-0509, DOI: 10.1016/j.procs.2016.09.246
YAGUE-MARTINEZ NESTOR ET AL: "Interferometric Processing of Sentinel-1 TOPS Data", IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, IEEE, USA, vol. 54, no. 4, 1 April 2016 (2016-04-01), pages 2220 - 2234, XP011602175, ISSN: 0196-2892, [retrieved on 20160309], DOI: 10.1109/TGRS.2015.2497902
SANTORO MAURIZIO: "Sentinel-1 processing with GAMMA software", 31 May 2015 (2015-05-31), XP093020720, Retrieved from the Internet [retrieved on 20230203]
Attorney, Agent or Firm:
REGIMBEAU (FR)
Download PDF:
Claims:
CLAIMS 1. A computer-implemented method for reconstructing a corrected primary mosaic image (MI2) from a set of primary burst images (Ii) of a ground surface area, wherein the method comprises : inputting (S0) the set of primary burst images (Ii) of the ground surface area, having been acquired by at least one acquisition swath carried out by a synthetic aperture radar carried by a moving satellite, and primary metadata (Bi(p), hi, li), which have been provided by the at least one acquisition swath carried out by the synthetic aperture radar, which are associated to the primary burst images (Ii) and which comprise a first primary origin (Bi(p)) of each primary burst image (Ii), a height (hi) of each primary burst image (Ii) and a width (li) of each primary burst image (Ii), wherein the primary burst images (Ii) are partially overlapping from one to another and have primary pixels (P), defining (S1) from the primary metadata (Bi(p), hi, li) and from the primary burst images (Ii) a primary mosaic image (MI) containing the primary burst images (Ii) and having primary coordinates (R, COL) of the primary pixels (P) defined relative to a second primary origin (O(p)) calculated in the primary mosaic image (MI), choosing (S2) points (Q) having each a longitude (Long), a latitude (Lat) and an altitude (h) in a digital elevation model (DEM) covering the primary mosaic image (MI), projecting (S3) the points (Q) into primary projected points (Q(p)) having the primary coordinates (R, COL) of some of the primary pixels (P) in the primary mosaic image (MI), calculating (S4) to which primary burst image (Ii) of the mosaic image (MI) belongs each primary projected point (Q(p)), applying (S5) primary corrections (C(Q(p), i, θ)) to the primary coordinates (R, COL) of the primary projected points (Q(p)) e primary mosaic image (MI) for calculating primary corrected coordinates (^^(^)) of the primary projected points (Q(p)) in the primary mosaic image (MI), calculating (S6), for each primary burst image (Ii) to which at least one (Qi(p)) of the primary projected points (Q(p)) belongs, a primary warp function (Ai(p)), which is different from the primary corrections (C(Q(p), i, θ)) and which transforms the primary coordinates (R, COL) of said at least one (Qi(p)) of the primary projected points (Q(p)) into the primary corrected co (^) ates (^^ ) of said at least one (Qi(p)) of the primary projected points (Q(p)) in the primary mosaic image (MI), transforming (S7) each primary burst image (Ii) by the primary warp function (Ai(p)) having been calculated for the primary burst image (Ii) and assembling (S7) into a corrected primary mosaic image (MI2) the primary burst images (Ii) having been transformed by the primary warp functions (Ai(p)) having been calculated. 2. The method of claim1, further comprising inputting (S10) a secondary set of secondary burst images (I’i) of the ground surface area, having been acquired by at least one acquisition swath carried out by the synthetic aperture radar carried by the moving satellite, and secondary metadata (Bi(s), h’i, l’i), which have been provided by the at least one acquisition swath carried out by the synthetic aperture radar, which are associated to the secondary burst images (I’i) and which comprise a first secondary origin (Bi(s)) of each secondary burst image (I’i), a height (h’i) of each secondary burst image (I’i) and a width (l’i) of each secondary burst image (I’i), wherein the secondary burst images (I’i) are partially overlapping from one to another and have secondary pixels (P’), defining (S11) from the secondary metadata (Bi(s), h’i, l’i) and from the secondary burst images (I’i) matching the primary burst images (Ii) a secondary mosaic image (MI’) containing the secondary burst images (I’i) and having secondary coordinates (R’, COL’) for the secondary pixels (P’) defined relative to a second secondary origin (O(s)) calculated in the secondary mosaic image (MI’), projecting (S13) the points (Q) into secondary projected points (Q(s)) having the secondary coordinates (R’, COL’) of some of the secondary pixels (P’) in the secondary mosaic image (MI’), calculating (S14) to which secondary burst image (I’i) of the secondary mosaic image (MI’) belongs each secondary projected point (Q(s)) by finding the secondary burst images (I’i) that match the primary burst images (Ii), applying (S15) secondary corrections (C(Q(s), i, θ)) to the secondary coordinates (R’, COL’) of the secondary projected points (Q(s)) in the secondary mosaic image (MI’) for calculating secondary corrected coordinates (^^(^)) of the secondary projected points (Q(s)) in the secondary mosaic image (MI’), calculating (S16), for each secondary burst image (I’i) to which at least one (Qi(s)) of the secondary projected points (Q(s)) belongs, a secondary warp function (Ai(s)), which transforms the primary coordinates (R, COL) of said at least one (Qi(s)) of the primary projected points (Q(p)) into the secondary corrected coordinates (^^(^)) of said at least one (Qi(s)) of the secondary projected points (Q(s)) in the secondary mosaic image (MI’), transforming (S17) each secondary burst image (I’i) by the secondary warp function (Ai(s)) having been calculated for the secondary burst image (I’i) and assembling (S17) into a secondary corrected mosaic image (MI2’) the secondary burst images (I’i) having been transformed by the secondary warp functions (Ai(s)) having been calculated. 3. The method according to claim 1 or 2, wherein the primary warp function (Ai(p)) is a primary affine resampling matrix for each primary burst image (Ii) and/or the secondary warp function (Ai(s)) is a secondary affine resampling matrix for each secondary burst image (I’i). 4. The method according to claim 3, wherein the primary warp function (Ai(p)) is defined as the primary affine resampling matrix transforming a first primary difference (Q(p)-Bi(p)) between the primary coordinates (R, COL) of said at least one (Qi(p)) of the primary projected points (Q(p)) and the first primary origin (Bi(p)) of the primary burst image (Ii) to which said at least one the primary projected points (Q(p)) belongs, in second primary difference (^^(^)- Bi(p)) between the primary corrected coordinates (^^(^)) of said at least one (Qi(p)) of the primary projected points and the first primary origin (Bi(p)) of the primary burst image (Ii) to which said at least one (Qi(p)) of the primary projected points (Q(p)) belongs. and/or the secondary warp function (Ai(s)) is defined as the secondary affine resampling matrix transforming a first secondary difference (Q(p)-Bi(p)) between the primary coordinates (R, COL) of said at least one (Qi(p)) of the primary projected points (Q(p)) and the first primary origin (Bi(p)) of the primary burst image (Ii) to which said at least one (Qi(s)) of the secondary projected points (Q(s)) belongs, into a second s ( difference (^^ ^)-Bi(s)) between the secondary co coordinates (^^(^)) of said at least one (Qi(s)) of the secondary projected points and the first secondary origin (Bi(s)) of the secondary burst image (Ii) to which said at least one (Qi(s)) of the secondary projected points (Q(s)) belongs. 5. The method according to claim 1, wherein the primary warp function (Ai(p)) is a primary homography resampling matrix for each primary burst image (Ii) and/or the secondary warp function (Ai(s)) is a secondary homography resampling matrix for each secondary burst image (I’i). 6. The method according to claim 1, wherein the primary warp function (Ai(p)) is a primary polynomial function for each primary burst image (Ii) and/or the secondary warp function (Ai(s)) is a secondary polynomial function for each secondary burst image (I’i). 7. The method according to any one of the preceding claims, wherein the primary corrections (C(Q(p), i, θ)) of the primary coordinates (R, COL) of the primary projected points (Q(p)) in the primary mosaic image (MI) varies depending on at least one primary parameter (cosαinc; Rank, τmid, IW2, EW3, PRI ; Kr, Nlburst, Fa, t0,i, fηc,i(τi), ka,i(τi), vs,i, fc, kψ,i, Nsswath) calculated from the primary metadata or contained in the primary metadata, and/or the secondary corrections (C(Q(s), i, θ)) of the secondary coordinates (R’, COL’) of the secondary projected points (Q(s)) in the secondary mosaic image (MI’) varies depending on at least one secondary parameter (cosα’inc; Rank’, τ’mid, IW2, EW3, PRI’ ; K’r, N’lburst, F’a, t’0,i, f’ηc,ii), k’a,ii), v’s,i, f’c, k’ψ,i, Ns’swath) calculated from the secondary metadata or contained in the secondary metadata.

8. The method according to claim 7, wherein the primary corrections (C(Q(p), i, θ)) of the primary coordinates (R, COL) of the primary projected points (Q(p)) in the primary mosaic image (MI) comprise a primary function (δτAPD, ψtropo) of correction of an atmospheric path delay between the synthetic aperture radar and the ground surface area, wherein the primary function (δτAPD, ψtropo) of correction of the atmospheric path delay varies at least with the altitude (h) of the points (Q) corresponding to the primary projected points (Q(p)) and with a first primary parameter (cosαinc) calculated from the primary metadata, and/or the secondary corrections (C(Q(s), i, θ)) of the secondary coordinates (R’, COL’) of the secondary projected points (Q(s)) in the secondary mosaic image (MI’) comprise a secondary function (δτ’APD, ψ’tropo) of correction of an atmospheric path delay between the synthetic aperture radar and the ground surface area, wherein the secondary function (δτ’APD, ψ’tropo) of correction of the atmospheric path delay varies at least with the altitude (h) of the points (Q) corresponding to the secondary projected points (Q(s)) and with a first secondary parameter (cosα’inc) calculated from the secondary metadata. 9. The method according to claim 8, wherein the primary function (δτAPD, ψtropo) of correction of the atmospheric path delay comprises a primary polynomial function (ψtropo,zenith) of the altitude (h) of the points (Q) corresponding to the primary projected points (Q(p)), which has a degree equal to 2, the first primary parameter is cosαinc, wherein αinc is a primary incidence angle between a line of sight of the synthetic aperture radar and a local normal to the ground surface area, the primary function (δτAPD, ψtropo) of correction of the atmospheric path delay is proportional to the primary polynomial function (ψtropo,zenith) divided by cosαinc, and/or the secondary function (δτ’APD, ψ’tropo) of correction of the atmospheric path delay comprises a secondary polynomial function (ψ’tropo,zenith) of the altitude (h) of the points (Q) corresponding to the secondary projected points (Q(s)), which has a degree equal to 2, the first secondary parameter is cosα’inc, wherein α’inc is a secondary incidence angle between a line of sight of the synthetic aperture radar and a local normal to the ground surface area, the secondary function (δτ’APD, ψ’tropo) of correction of the atmospheric path delay is proportional to the secondary polynomial function (ψ’tropo,zenith) divided by cosα’inc. 10. The method according to any one of claims 7 to 9, wherein the primary corrections (C(Q(p), i, θ)) of the primary coordinates (R, COL) of the primary projected points (Q(p)) in the primary mosaic image (MI) comprise a primary function (δtbistatic) of bistatic correction of an azimuth time (t) of the synthetic aperture radar along an azimuth direction (T), along which the moving satellite moves, between transmision of a pulse by the synthetic aperture radar and reception of the pulse by the synthetic aperture radar, and/or the secondary corrections (C(Q(s), i, θ)) of the secondary coordinates (R’, COL’) of the secondary projected points (Q(s)) in the secondary mosaic image (MI’) comprise a secondary function (δt’bistatic) of bistatic correction of an azimuth time (t’) of the synthetic aperture radar along an azimuth direction (T), along which the moving satellite moves, between transmision of a pulse by the synthetic aperture radar and reception of the pulse by the synthetic aperture radar. 11. The method according to claim 10, wherein the primary coordinates (R, COL) comprises a primary column number (COL) of the primary pixels (P) in the primary mosaic image (MI), the primary function (δtbistatic) of bistatic correction of azimuth time (t) of the synthetic aperture radar varies in an affine manner with the primary column number (COL) of the primary pixels (P) in the primary mosaic image (MI) and with primary second parameters (Rank, τmid, IW2, EW3, PRI) given in the primary metadata or calculated from the primary metadata, and/or the secondary coordinates (R’, COL’) comprises a secondary column number (COL’) of the secondary pixels (P’) in the secondary mosaic image (MI’), the secondary function (δt’bistatic) of bistatic correction of azimuth time (t’) of the synthetic aperture radar varies in an affine manner with the secondary column number (COL’) of the secondary pixels (P’) in the secondary mosaic image (MI’) and with secondary second parameters (Rank’, τ’mid, IW2, EW3, PRI’) given in the secondary metadata or calculated from the secondary metadata. 12. The method according to any one of claims 7 to 11, wherein the primary corrections (C(Q(p), i, θ)) of the primary coordinates (R, COL) of the primary projected points (Q(p)) in the primary mosaic image (MI) comprise a primary function (δτintra_pulse,i) of intra pulse motion correction, and/or the secondary corrections (C(Q(s), i, θ)) of the secondary coordinates (R’, COL’) of the secondary projected points (Q(s)) in the secondary mosaic image (MI’) comprise a secondary function (δτ’intra_pulse,i) of intra pulse motion correction. 13. The method according to claim 12, wherein the primary coordinates (R, COL) comprises a primary column number (COLi) of the primary pixels (P) in the primary burst image (Ii) and a primary row number (Ri ) of the primary pixels (P) in the primary burst image (Ii), the primary function (δτintra_pulse,i) of intra pulse motion correction is negatively proportional to a primary Doppler centroid frequency (fDCi,ti)) calculated depending on primary third parameters (Kr, Nlburst, Fa, t0,i, fηc,ii), ka,i(τi), vs,i, fc, kψ,i, Nsswath) calculated from the primary metadata or given in the primary metadata and in an affine manner on the primary row number (Ri ) of the primary pixels (P) in the primary burst image (Ii), wherein at least one (fηc,i(τi), ka,i(τi)) of the primary third parameters depends on the primary column number (COLi) of the primary pixels (P) in the primary burst image (Ii) and/or the secondary coordinates (R’, COL’) comprises a secondary column number (COL’i) of the secondary pixels (P’) in the secondary burst image (I’i) and a secondary row number (R’i) of the secondary pixels (P’) in the secondary burst image (I’i), the secondary function (δτ’intra_pulse,i) of intra pulse motion correction is negatively proportional to a secondary Doppler centroid frequency (f’DCi,ti)) calculated depending on secondary third parameters (K’r, N’lburst, F’a, t’0,i, f’ηc,i(τi), k’a,i(τi), v’s,i, f’c, k’ψ,i, Ns’swath) calculated from the secondary metadata or given in the secondary metadata and in an affine manner on the secondary row number (R’i) of the secondary pixels (P’) in the secondary burst image (I’i), wherein at least one (f’ηci), k’ai)) of the secondary third parameters depends on the secondary column number (COL’i) of the secondary pixels (P’) in the secondary burst image (I’i). 14. The method according to any one of the preceding claims, wherein choosing (S2) points (Q) having each a longitude (Long), a latitude (Lat) and an altitude (h) in a digital elevation model (DEM) covering the primary mosaic image (MI) comprises reading from the primary metadata a geographical extent of the primary mosaic image (MI) and downloading (S2) from the primary metadata (Bi(p), hi, li) the digital elevation model (DEM) giving for each of the points (Q) of the digital elevation model (DEM) the longitude (Long) of the point (Q), the latitude (Lat) of the point (Q) and the altitude (h) of the point (Q). 15. The method according to any one of the preceding claims, wherein defining (S1) the primary mosaic image (MI) comprises calculating (S1) from the primary metadata (Bi(p), hi, li) for each primary couple (i, i+1) of primary burst images (Ii, Ii+1) which are partially overlapping from one to another, a primary overlapping zone (Zi,i+1) common to the primary burst images (Ii, Ii+1) of the primary couple (i, i+1) and a primary cut (Di,i+1), which is situated in the primary overlapping zone (Zi,i+1), wherein calculating (S4) to which primary burst image (Ii) of the primary mosaic image (MI) belongs each primary projected point (Q(p)) is carried out based on the primary burst images (Ii, Ii+1) delimited at least by the primary cuts (Di,i+1). 16. The method according to any one of the preceding claims, wherein defining (S1) the primary mosaic image (MI) comprises calculating (S1) from the primary metadata (Bi(p), hi, li) for each primary couple (i, i+1) of primary burst images (Ii, Ii+1) which are partially overlapping from one to another, a primary overlapping zone (Zi,i+1) common to the primary burst images (Ii, Ii+1) of the primary couple (i, i+1) and a primary cut (Di,i+1), which is situated in the primary overlapping zone (Zi,i+1), wherein assembling (S7) into the primary corrected mosaic image (MI2) the primary burst images (Ii) having been transformed by the primary warp functions (Ai(p)) having been calculated is carried out using the primary cuts (Di,i+1) delimiting the primary burst images (Ii) in the primary corrected mosaic image (MI2), and/or wherein assembling (S17) into the secondary corrected mosaic image (MI2’) the secondary burst images (I’i) having been transformed by the secondary warp functions (Ai(s)) having been calculated is carried out using the primary cuts (Di,i+1) delimiting the primary burst image s ( i) n the secondary corrected mosaic image (MI2’). 17. The method according to any one of the preceding claims, wherein the primary coordinates (R, COL) of the primary projected points (Q(p)) are defined by reference to the second primary origin (O(p)) calculated in the primary mosaic image (MI), the primary corrected coordinates (^^(^)) of the primary projected points (Q(p)) are defined by reference to the primary origin (O(p)) calculated in the primary mosaic image (MI), and/or the secondary coordinates (R’, COL’) of the secondary projected points (Q(s)) are defined by reference to the second secondary origin (O(s)) calculated in the secondary mosaic image (MI’), the secondary corrected coordinates (^^(^)) of the secondary projected points (Q(s)) are defined by reference to the secondary origin (O(s)) calculated in the secondary mosaic image (MI’). 18. A device for reconstructing a corrected primary mosaic image (MI2) from a set of primary burst images (Ii) of a ground surface area, wherein the device comprises a calculator configured to : input (S0) the set of primary burst images (Ii) of the ground surface area, having been acquired by at least one acquisition swath carried out by a synthetic aperture radar carried by a moving satellite, and primary metadata (Bi(p), hi, li), which have been provided by the at least one acquisition swath carried out by the synthetic aperture radar, which are associated to the primary burst images (Ii) and which comprise a first primary origin (Bi(p)) of each primary burst image (Ii), a height (hi) of each primary burst image (Ii) and a width (li) of each primary burst image (Ii), wherein the primary burst images (Ii) are partially overlapping from one to another and have primary pixels (P), define (S1) from the primary metadata (Bi(p), hi, li) and from the primary burst images (Ii) the primary mosaic image (MI) containing the primary burst images (Ii) and having primary coordinates (R, COL) of the primary pixels (P) defined relative to a second primary origin (O(p)) calculated in the primary mosaic image (MI), choose (S2) points (Q) having each a longitude (Long), a latitude (Lat) and an altitude (h) in a digital elevation model (DEM) covering the primary mosaic image (MI), project (S3) the points (Q) into prim ry projected points (Q(p)) having the primary coordinates (R, COL) of some of the primary pixels (P) in the primary mosaic image (MI), calculate (S4) to which primary burst image (Ii) of the mosaic image (MI) belongs each primary projected point (Q(p)), apply (S5) primary corrections (C(Q(p), i, θ)) to the primary coordinates (R, COL) of the primary projected points (Q(p)) in the primary mo mage (MI) for calculating primary corrected coordinates (^^(^)) of the primary projected points (Q(p)) in the primary mosaic image (MI), calculate (S6), for each primary burst image (Ii) to which at least one (Qi(p)) of the primary projected points (Q(p)) belongs, a primary warp function (Ai(p)), which is different from the primary corrections (C(Q(p), i, θ)) and which transforms the primary coordinates (R, COL) of said at least one (Qi(p)) of the primary projected points (Q(p)) into the primary corrected coordinates (^^(^)) of said at least one (Qi(p)) of the primary projected points (Q(p)) in the primary mosaic image (MI), transform (S7) each primary burst image (Ii) by the primary warp function (Ai(p)) having been calculated for the primary burst image (Ii) and assembling (S7) into a corrected primary mosaic image (MI2) the primary burst images (Ii) having been transformed by the primary warp functions (Ai(p)) having been calculated. 19. A non-transitory computer-readable storage medium having instructions that, when executed by a processor, cause the processor to execute operations for reconstructing a corrected primary mosaic image (MI2) from a set of primary burst images (Ii) of a ground surface area, wherein the operations comprise inputting (S0) the set of primary burst images (Ii) of the ground surface area, having been acquired by at least one acquisition swath carried out by a synthetic aperture radar carried by a moving satellite, and primary metadata (Bi(p), hi, li), which have been provided by the at least one acquisition swath carried out by the synthetic aperture radar, which are associated to the primary burst images (Ii) and which comprise a first primary origin (Bi(p)) of each primary burst image (Ii), a height (hi) of each primary burst image (Ii) and a width (li) of each primary burst image (Ii), wherein the primary burst images (Ii) are partially overlapping from one to another and have primary pixels (P), defining (S1) from the primary metadata (Bi(p), hi, li) and from the primary burst images (Ii) the primary mosaic image (MI) containing the primary burst images (Ii) and having primary coordinates (R, COL) of the primary pixels (P) defined relative to a second primary origin (O(p)) calculated in the primary mosaic image (MI), choosing (S2) points (Q) having each a longitude (Long), a latitude (Lat) and an altitude (h) in a digital elevation model (DEM) covering the primary mosaic image (MI), projecting (S3) the points (Q) into primary projected points (Q(p)) having the primary coordinates (R, COL) of some of the primary pixels (P) in the primary mosaic image (MI), calculating (S4) to which primary burst image (Ii) of the mosaic image (MI) belongs each primary projected point (Q(p)), applying (S5) primary corrections (C(Q(p), i, θ)) to the primary coordinates (R, COL) of the primary projected points (Q(p)) in the primary mosaic image (MI) for calculating primary corrected coordinates (^^(^)) of the i ary projected points (Q(p)) in the primary mosaic image (MI), calculating (S6), for each primary burst image (Ii) to which at least one (Qi(p)) of the primary projected points (Q(p)) belongs, a primary warp function (Ai(p)), which is different from the primary corrections (C(Q(p), i, θ)) and which transforms the primary coordinates (R, COL) of said at least one (Qi(p)) of the primary projected points (Q(p)) into the primary corrected coordinates (^^(^)) of said at least one (Qi(p)) of the primary projected points (Q(p)) in the primary mosaic image (MI), transforming (S7) each primary burst image (Ii) by the primary warp function (Ai(p)) having been calculated for the primary burst image (Ii) and assembling (S7) into a corrected primary mosaic image (MI2) the primary burst images (Ii) having been transformed by the primary warp functions (Ai(p)) having been calculated.

Description:
Method, device and storage medium for reconstructing a corrected mosaic image from burst images BACKGROUND OF THE INVENTION The invention concerns a computer-implemented method, a device and a non-transitory computer-readable storage medium for reconstructing a corrected primary mosaic image from a set of primary burst images of a ground surface area. A field of application of the invention is earth monitoring, for example for the detection of human activity or of any activity, by burst images of the ground surface area, having been acquired by a remote synthetic aperture radar (SAR) carried by a moving satellite. For example, the synthetic aperture radar may be carried by the satellite known as Sentinel-1, which would acquire the image in bursts when operating in the TOPSAR mode (meaning Terrain Observation by Progressive Scans of Synthetic Aperture Radar). A problem is that the burst images provided by this synthetic aperture radar are affected by defects of geometric precision, which have for consequences geolocation errors when using the burst images. The burst images can be acquired by the synthetic aperture radar in an interferometric wide (IW) swath mode providing a set of burst images along the trajectory (orbit) of the satellite, which is the azimuth coordinate. The synthetic aperture radar may obtain 3 subswaths during a same passage by orienting the synthetic aperture radar in different orientations transversally to the trajectory during the travel of the satellite along the trajectory, wherein each subswath acquires its own set of burst images. The added difficulty of this swath mode to the processing is the need to mosaic the burst images into a continuous image (mosaic image). Then, the problem arises when shifts depending on the burst images are present in the data, which would cause geometric discontinuities and inconsistencies at the burst boundaries in the mosaic image. A goal of the invention is to propose a computer-implemented method, a device and a non-transitory computer-readable storage medium for reconstructing a corrected primary mosaic image from a set of primary burst images of a ground surface area, which circumvent this difficulty and problem. SUMMARY OF THE INVENTION A first subject matter of the invention is a computer-implemented method A computer-implemented method for reconstructing a corrected primary mosaic image from a set of primary burst images of a ground surface area, wherein the method comprises : inputting the set of primary burst images of the ground surface area, having been acquired by at least one acquisition swath carried out by a synthetic aperture radar carried by a moving satellite, and primary metadata, which have been provided by the at least one acquisition swath carried out by the synthetic aperture radar, which are associated to the primary burst images and which comprise a first primary origin of each primary burst image, a height of each primary burst image and a width of each primary burst image, wherein the primary burst images are partially overlapping from one to another and have primary pixels, defining from the primary metadata and from the primary burst images a primary mosaic image containing the primary burst images and having primary coordinates of the primary pixels defined relative to a second primary origin calculated in the primary mosaic image, choosing points having each a longitude, a latitude and an altitude in a digital elevation model covering the primary mosaic image, projecting the points into primary projected points having the primary coordinates of some of the primary pixels in the primary mosaic image, calculating to which primary burst image of the mosaic image belongs each primary projected point, applying primary corrections to the primary coordinates of the primary projected points in the primary mosaic image for calculating primary corrected coordinates of the primary projected points in the primary mosaic image, calculating, for each primary burst image to which at least one of the primary projected points belongs, a primary warp function, which is different from the primary corrections and which transforms the primary coordinates of said at least one of the primary projected points into the primary corrected coordinates of said at least one of the primary projected points in the primary mosaic image, transforming each primary burst image by the primary warp function having been calculated for the primary burst image and assembling into a corrected primary mosaic image the primary burst images having been transformed by the primary warp functions having been calculated. Thanks to the method according to the invention, for reconstructing the corrected primary mosaic image from the set of primary burst images of the ground surface area are simple and generic and correct geometric errors in the burst images before assembling. The correction is based on the geolocation of a set of points sampled from a digital elevation model. Thanks to the invention, the difference in the shifts from one burst to another are corrected. According to an embodiment of the invention, inputting a secondary set of secondary burst images of the ground surface area, having been acquired by at least one acquisition swath carried out by the synthetic aperture radar carried by the moving satellite, and secondary metadata, which have been provided by the at least one acquisition swath carried out by the synthetic aperture radar, which are associated to the secondary burst images and which comprise a first secondary origin of each secondary burst image, a height of each secondary burst image and a width of each secondary burst image, wherein the secondary burst images are partially overlapping from one to another and have secondary pixels, defining from the secondary metadata and from the secondary burst images matching the primary burst images a secondary mosaic image containing the secondary burst images and having secondary coordinates for the secondary pixels defined relative to a second secondary origin calculated in the secondary mosaic image, projecting the points into secondary projected points having the secondary coordinates of some of the secondary pixels in the secondary mosaic image, calculating to which secondary burst image of the secondary mosaic image belongs each secondary projected point by finding the secondary burst images that match the primary burst images, applying secondary corrections to the secondary coordinates of the secondary projected points in the secondary mosaic image for calculating secondary corrected coordinates of the secondary projected points in the secondary mosaic image, calculating, for each secondary burst image to which at least one of the secondary projected points belongs, a secondary warp function, which transforms the primary coordinates of said at least one of the secondary projected points into the secondary corrected coordinates of said at least one of the secondary projected points in the secondary mosaic image, transforming each secondary burst image by the secondary warp function having been calculated for the secondary burst image and assembling into a secondary corrected mosaic image the secondary burst images having been transformed by the secondary warp functions having been calculated. According to an embodiment of the invention, the primary warp function is a primary affine resampling matrix for each primary burst image and/or the secondary warp function is a secondary affine resampling matrix for each secondary burst image. According to an embodiment of the invention, the primary warp function is defined as the primary affine resampling matrix transforming a first primary difference between the primary coordinates of said at least one of the primary projected points and the first primary origin of the primary burst image to which said at least one of the primary projected points belongs, into a second primary difference between the primary corrected coordinates of said at least one of the primary projected points and the first primary origin of the primary burst image to which said at least one of the primary projected points belongs. and/or the secondary warp function is defined as the secondary affine resampling matrix transforming a first secondary difference between the primary coordinates of said at least one of the primary projected points and the first primary origin of the primary burst image to which said at least one of the primary projected points belongs, into a second secondary difference between the secondary corrected coordinates of said at least one of the secondary projected points and the first secondary origin of the secondary burst image to which said at least one of the secondary projected points belongs. According to an embodiment of the invention, the primary warp function is a primary homography resampling matrix for each primary burst image and/or the secondary warp function is a secondary homography resampling matrix for each secondary burst image. According to an embodiment of the invention, the primary warp function is a primary polynomial function for each primary burst image and/or the secondary warp function is a secondary polynomial function for each secondary burst image. According to an embodiment of the invention, the primary corrections of the primary coordinates of the primary projected points in the primary mosaic image varies depending on at least one primary parameter calculated from the primary metadata or contained in the primary metadata, and/or the secondary corrections of the secondary coordinates of the secondary projected points in the secondary mosaic image varies depending on at least one secondary parameter calculated from the secondary metadata or contained in the secondary metadata. According to an embodiment of the invention, the primary corrections of the primary coordinates of the primary projected points in the primary mosaic image comprise a primary function of correction of an atmospheric path delay between the synthetic aperture radar and the ground surface area, wherein the primary function of correction of the atmospheric path delay varies at least with the altitude of the points corresponding to the primary projected points and with a first primary parameter calculated from the primary metadata, and/or the secondary corrections of the secondary coordinates of the secondary projected points in the secondary mosaic image comprise a secondary function of correction of an atmospheric path delay between the synthetic aperture radar and the ground surface area, wherein the secondary function of correction of the atmospheric path delay varies at least with the altitude of the points corresponding to the secondary projected points and with a first secondary parameter calculated from the secondary metadata. According to an embodiment of the invention, the primary function of correction of the atmospheric path delay comprises a primary polynomial function of the altitude of the points corresponding to the primary projected points, which has a degree equal to 2, the first primary parameter is cosα inc , wherein α inc is a primary incidence angle between a line of sight of the synthetic aperture radar and a local normal to the ground surface area, the primary function of correction of the atmospheric path delay is proportional to the primary polynomial function divided by cosαinc, and/or the secondary function of correction of the atmospheric path delay comprises a secondary polynomial function of the altitude of the points corresponding to the secondary projected points, which has a degree equal to 2, the first secondary parameter is cosα’inc, wherein α’inc is a secondary incidence angle between a line of sight of the synthetic aperture radar and a local normal to the ground surface area, the secondary function of correction of the atmospheric path delay is proportional to the secondary polynomial function divided by cosα’ inc . According to an embodiment of the invention, the primary corrections of the primary coordinates of the primary projected points in the primary mosaic image comprise a primary function of bistatic correction of an azimuth time of the synthetic aperture radar along an azimuth direction, along which the moving satellite moves, between transmision of a pulse by the synthetic aperture radar and reception of the pulse by the synthetic aperture radar, and/or the secondary corrections of the secondary coordinates of the secondary projected points in the secondary mosaic image comprise a secondary function of bistatic correction of an azimuth time of the synthetic aperture radar along an azimuth direction, along which the moving satellite moves, between transmision of a pulse by the synthetic aperture radar and reception of the pulse by the synthetic aperture radar. According to an embodiment of the invention, the primary coordinates comprises a primary column number of the primary pixels in the primary mosaic image, the primary function of bistatic correction of azimuth time of the synthetic aperture radar varies in an affine manner with the primary column number of the primary pixels in the primary mosaic image and with primary second parameters given in the primary metadata or calculated from the primary metadata, and/or the secondary coordinates comprises a secondary column number of the secondary pixels in the secondary mosaic image, the secondary function of bistatic correction of azimuth time of the synthetic aperture radar varies in an affine manner with the secondary column number of the secondary pixels in the secondary mosaic image and with secondary second parameters given in the secondary metadata or calculated from the primary metadata. According to an embodiment of the invention, the primary corrections of the primary coordinates of the primary projected points in the primary mosaic image comprise a primary function of intra pulse motion correction, and/or the secondary corrections of the secondary coordinates of the secondary projected points in the secondary mosaic image comprise a secondary function of intra pulse motion correction. According to an embodiment of the invention, the primary coordinates comprises a primary column number of the primary pixels in the primary burst image and a primary row number of the primary pixels in the primary burst image, the primary function of intra pulse motion correction is negatively proportional to a primary Doppler centroid frequency calculated depending on primary third parameters calculated from the primary metadata or given in the primary metadata and in an affine manner on the primary row number of the primary pixels in the primary burst image, wherein at least one of the primary third parameters depends on the primary column number of the primary pixels in the primary burst image, and/or the secondary coordinates comprises a secondary column number of the secondary pixels in the secondary burst image and a secondary row number of the secondary pixels in the secondary burst image, the secondary function of intra pulse motion correction is negatively proportional to a secondary Doppler centroid frequency calculated depending on secondary third parameters calculated from the secondary metadata or given in the secondary metadata and in an affine manner on the secondary row number of the secondary pixels in the secondary burst image, wherein at least one of the secondary third parameters depends on the secondary column number of the secondary pixels in the secondary burst image. According to an embodiment of the invention, choosing points having each a longitude, a latitude and an altitude in a digital elevation model covering the primary mosaic image comprises reading from the primary metadata a geographical extent of the primary mosaic image and downloading from the primary metadata the digital elevation model giving for each of the points of the digital elevation model the longitude of the point, the latitude of the point and the altitude of the point. According to an embodiment of the invention, defining the primary mosaic image comprises calculating from the primary metadata for each primary couple of primary burst images which are partially overlapping from one to another, a primary overlapping zone common to the primary burst images of the primary couple and a primary cut, which is situated in the primary overlapping zone, wherein calculating to which primary burst image of the primary mosaic image belongs each primary projected point is carried out based on the primary burst images delimited at least by the primary cuts. According to an embodiment of the invention, defining the primary mosaic image comprises calculating from the primary metadata for each primary couple of primary burst images which are partially overlapping from one to another, a primary overlapping zone common to the primary burst images of the primary couple and a primary cut, which is situated in the primary overlapping zone, wherein assembling into the primary corrected mosaic image the primary burst images having been transformed by the primary warp functions having been calculated is carried out using the primary cuts delimiting the primary burst images in the primary corrected mosaic image, and/or wherein assembling into the secondary corrected mosaic image the secondary burst images having been transformed by the secondary warp functions having been calculated is carried out using the primary cuts delimiting the primary burst images in the secondary corrected mosaic image. According to an embodiment of the invention, the primary coordinates of the primary projected points are defined by reference to the second primary origin calculated in the primary mosaic image, the primary corrected coordinates of the primary projected points are defined by reference to the primary origin calculated in the primary mosaic image, and/or the secondary coordinates of the secondary projected points are defined by reference to the second secondary origin calculated in the secondary mosaic image, the secondary corrected coordinates of the secondary projected points are defined by reference to the secondary origin calculated in the secondary mosaic image. A second subject matter of the invention is a device for reconstructing a corrected primary mosaic image from a set of primary burst images of a ground surface area, wherein the device comprises a calculator configured to : input the set of primary burst images of the ground surface area, having been acquired by at least one acquisition swath carried out by a synthetic aperture radar carried by a moving satellite, and primary metadata, which have been provided by the at least one acquisition swath carried out by the synthetic aperture radar, which are associated to the primary burst images and which comprise a first primary origin of each primary burst image, a height of each primary burst image and a width of each primary burst image, wherein the primary burst images are partially overlapping from one to another and have primary pixels, define from the primary metadata and from the primary burst images the primary mosaic image containing the primary burst images and having primary coordinates of the primary pixels defined relative to a second primary origin calculated in the primary mosaic image, choose points having each a longitude, a latitude and an altitude in a digital elevation model covering the primary mosaic image, project the points into primary projected points having the primary coordinates of some of the primary pixels in the primary mosaic image, calculate to which primary burst image of the mosaic image belongs each primary projected point, apply primary corrections to the primary coordinates of the primary projected points in the primary mosaic image for calculating primary corrected coordinates of the primary projected points in the primary mosaic image, calculate, for each primary burst image to which at least one of the primary projected points belongs, a primary warp function, which is different from the primary corrections and which transforms the primary coordinates of said at least one of the primary projected points into the primary corrected coordinates of said at least one of the primary projected points in the primary mosaic image, transform each primary burst image by the primary warp function having been calculated for the primary burst image and assembling into a corrected primary mosaic image the primary burst images having been transformed by the primary warp functions having been calculated. A third subject matter of the invention is a non-transitory computer- readable storage medium having instructions that, when executed by a processor, cause the processor to execute operations for reconstructing a corrected primary mosaic image from a set of primary burst images of a ground surface area, wherein the operations comprise inputting the set of primary burst images of the ground surface area, having been acquired by at least one acquisition swath carried out by a synthetic aperture radar carried by a moving satellite, and primary metadata, which have been provided by the at least one acquisition swath carried out by the synthetic aperture radar, which are associated to the primary burst images and which comprise a first primary origin of each primary burst image, a height of each primary burst image and a width of each primary burst image, wherein the primary burst images are partially overlapping from one to another and have primary pixels, defining from the primary metadata and from the primary burst images the primary mosaic image containing the primary burst images and having primary coordinates of the primary pixels defined relative to a second primary origin calculated in the primary mosaic image, choosing points having each a longitude, a latitude and an altitude in a digital elevation model covering the primary mosaic image, projecting the points into primary projected points having the primary coordinates of some of the primary pixels in the primary mosaic image, calculating to which primary burst image of the mosaic image belongs each primary projected point, applying primary corrections to the primary coordinates of the primary projected points in the primary mosaic image for calculating primary corrected coordinates of the primary projected points in the primary mosaic image, calculating, for each primary burst image to which at least one of the primary projected points belongs, a primary warp function, which is different from the primary corrections and which transforms the primary coordinates of said at least one of the primary projected points into the primary corrected coordinates of said at least one of the primary projected points in the primary mosaic image, transforming each primary burst image by the primary warp function having been calculated for the primary burst image and assembling into a corrected primary mosaic image the primary burst images having been transformed by the primary warp functions having been calculated. DESCRIPTION OF THE DRAWINGS The invention will be more clearly understood from the following description, given solely by way of non–limiting example in reference to the appended drawings, in which : - Figure 1 schematically illustrates burst images of the ground surface area, which have been acquired by a remote synthetic aperture radar carried by a moving satellite images and which can be processed by the method, the device and the non-transitory computer-readable storage medium for reconstructing a corrected mosaic image according to the invention, - Figure 2 schematically illustrates operations carried out by the method, the device and the non-transitory computer-readable storage medium for reconstructing a corrected primary mosaic image from primary burst images according to embodiments of the invention, - Figure 3 schematically illustrates operations carried out by the method, the device and the non-transitory computer-readable storage medium for reconstructing a corrected primary mosaic image from primary burst images and for reconstructing a corrected secondary mosaic image from secondary burst images according to embodiments of the invention, - Figure 4 schematically illustrates steps of the method for reconstructing a corrected mosaic image according to embodiments of the invention, - Figure 5 schematically illustrates a device a corrected mosaic image the method according to embodiments of the invention. DETAILED DESCRIPTION OF THE INVENTION Embodiments of the method and device 100 according to the invention for reconstructing a corrected mosaic image MI2 of a ground surface area of the earth or scene or terrestrial location from a set of N primary burst images Ii are described below in reference to Figures 1 to 5. The index i or identifier i uniquely identifies the primary burst image I i in the set of N primary burst images I i . The indexes i of the primary burst images I i may be an integer going from 1 to N. The integer N is higher than or equal to 2 or 3. The primary burst images are designated by Ii and Ii+1 and are partially overlapping for i going from 1 to N-1, as shown on Figures 2 and 3. The primary burst image Ii is also called respective primary burst image Ii and the primary burst image Ii+1 is also called next primary burst image Ii+1. Each primary burst image Ii has primary pixels P. The primary burst images I i are previously acquired by a remote synthetic aperture radar (SAR, as described above) carried by a moving satellite. The primary burst images Ii are satellite radar images. As illustrated in Figure 1, the satellite is in orbit around the earth and moves along a satellite flight path T, i.e. along the azimuth direction T relative to the earth. The synthetic aperture radar carried by the moving satellite acquires the set of primary burst image Ii during an acquisition swath which corresponds to the travel of the satellite along a certain distance (for example the ground surface area covered by the 9 primary burst images Ii shown in Figure 1 extends approximatively on a distance of 200 km on earth) and during a certain duration along the azimuth direction T relative to the earth. Consequently, the primary burst image I i are acquired at different positions of the satellite along the azimuth direction T. For example, the synthetic aperture radar may acquire several (for example 3 in Figure 1) subswaths IW1, IW2, IW3 during a same swath by orienting the synthetic aperture radar in different orientations transversally to the azimuth direction T during the travel of the satellite along the azimuth direction T. Each subswath IW1, IW2, IW3 acquires its own set of primary burst images I i . The ground surface areas covered by the subswaths IW1, IW2, IW3 of primary burst images Ii are different from each other and may be adjacent to each other or in the vicinity of each other, as illustrated in Figure 1. When, after a certain time, the satellite travels again over the same ground surface area, the synthetic aperture radar carried by the moving satellite may acquire the set of secondary burst image I’ i during another acquisition swath which corresponds to the travel of the satellite along a certain distance during a certain duration along the azimuth direction T relative to the earth and may acquire several (for example 3) other subswaths IW1, IW2, IW3 during this other swath. The secondary burst images I’i are acquired at different positions of the satellite along the azimuth direction T. Each secondary burst image I’ i has secondary pixels P’. The primary burst images I i and/or the secondary burst images I’ i may also be acquired in an EW mode (EW meaning Extra Wide swath). In a general manner, what is described above and below for the primary burst images Ii is also valid for the secondary burst images I’i, except otherwise indicated. The primary burst images Ii and/or the secondary burst images I’i are previously sent by the satellite to one or more stations located on earth and are stored in a media. The primary burst images I i and/or the secondary burst images I’i may be obtained by an image production module, which may be a computer or any automatic machine, able to download the images from the media. The media may be an external media, which may be a database, a distant server, a website, a distant computer or others. The media may be a local media, such as memory or storage media, for example of an image production module. Figures 2 and 3 show that the primary burst image Ii-1 is partially overlapping the primary burst image Ii and that the primary burst images Ii-1 and Ii cover an overlapping geographical zone Z i-1,i in common. The primary burst image I i is partially overlapping the primary burst image I i+1 and the primary burst images I i and Ii+1 cover an overlapping geographical zone Zi,i+1 in common. The primary burst images Ii partially overlap from each respective primary burst image to the next primary burst image in an overlap direction and the overlapping geographical zones Zi-1,i, Zi,i+1 are shifted in that overlap direction , so that they may form a strip of images extending in that overlap direction . For example, the overlap direction is perpendicular to the largest dimension (width li) of each primary burst images Ii, as shown on Figures 2 and 3. Of course, the overlap direction may be perpendicular to the smallest dimension (height hi) of each primary burst images Ii. For example, the overlap direction is along the azimuth direction T or is perpendicular to the azimuth direction T. A processing device 100 for processing the primary burst images I i and/or the secondary burst images I’ i to carry out the method(s) is also provided, as shown on Figure 5. The device 100 and method may be embodied by any means, such as for example computers, calculators, processors, microprocessors, permanent memories 103, servers, databases, computer programs, man-machine interfaces, user interfaces, network, framework, neural networks. The device 100 may comprise a calculator (one or more) 101 to carry out the steps and operations mentioned of the method of the invention, wherein this calculator 101 can comprise a computer (one or more), a processor (one or more), a microprocessor (one or more), a graphics processing unit (one or more), a memory 103 (one or more), a permanent or non-transitory computer-readable storage medium 103 (one or more), a permanent or non-transitory computer- readable storage memory 103 (one or more), a server (one or more), a database (one or more), a computer program (one or more), a network (one or more), a framework (one or more), a neural network (one or more), or others. The device 100 may comprise a physical input interface 102 (one or more), a physical output interface (one or more), a man-machine interface (one or more), a user interface (one or more) or others. Computer programs according to the invention may comprise instructions for executing the steps and operations of the method. The computer program may be recorded on any storage media, which may be permanent storage memory 103, i.e. a non-transitory computer-readable storage medium 103, or a non-permanent storage memory. The method according to the invention comprises the steps described below, as shown on Figure 4. In a step S0, the set of primary burst images Ii of the ground surface area, as well as primary metadata, which have been provided by the at least one acquisition swath carried out by the synthetic aperture radar in association with the primary burst images Ii are inputted to the input interface 102 of the device 100 and are stored in the permanent storage memory 103. The calculator 101 is connected to the input interface 102, to the permanent storage memory 103 and to the physical output interface 104. The physical output interface 104 may be for example a screen or a physical output of the device 100. The primary metadata comprise a first primary origin Bi (p) of each primary burst image Ii, a height hi of each primary burst image Ii and a width li of each primary burst image I i , which define the limits and burst size of each primary burst image I i , In a step S1, the calculator 101 defines from the primary metadata B i (p) , h i , l i and from the primary burst images Ii a primary mosaic image MI containing the primary burst images Ii. The calculator 101 calculates primary coordinates R (primary row number R), COL (primary column number COL) of the primary pixels P defined relative to a second primary origin O (p) in the primary mosaic image MI, such as for example primary coordinates R of rows of the pixels P and primary coordinates COL of columns of the pixels P. The primary coordinates R of rows of the pixels P may extend along the azimuth direction T. The primary coordinates COL of columns of the pixels P may extend along the range direction RD, which is perpendicular to the azimuth direction T in the plane of each primary burst image Ii. The calculator 101 thus determines the limits of the primary burst images Ii. In a step S2, the calculator 101 chooses points Q having each a longitude Long, a latitude Lat and an altitude h in a digital elevation model DEM covering the primary mosaic image MI. The calculator 101 may determine from the primary metadata Bi (p) , hi, li a geographical extent (minimum/maximum longitude, minimum/maximum latitude) of the primary mosaic image MI containing the primary burst images Ii (the geographic extent of the primary mosaic image may be deduced from other fields in the metadata). The calculator 101 may download the digital elevation model DEM covering the geographical extent of the primary mosaic image MI containing the primary burst images Ii. The calculator 101 may download the digital elevation model DEM using external tools. The digital elevation model DEM gives for each of the points Q of the digital elevation model DEM the longitude Long of the point Q, the latitude Lat of the point Q and the altitude h of the point Q. The primary coordinates R, COL of one of the primary pixels P in the primary mosaic image MI can be computed by a geolocation operation (in this case, a projection operation). The digital elevation model DEM gives the longitude, latitude and altitude of some or all points of a surface of the ground surface area covered by the primary mosaic image MI containing the primary burst images Ii. The points Q are chosen among the points of this surface of the digital elevation model DEM. In a step S3, the calculator 101 projects the points Q into primary projected points Q (p) having the primary coordinates R, COL of some of the primary pixels P in the primary mosaic image MI (The coordinates R, COL may be real numbers and not necessarily integers). The calculator 101 calculates for the set of the points Q of the digital elevation model DEM the primary coordinates R, COL of the primary pixels P of the primary projected points Q (p) in the mosaic image MI, which are corresponding according to the digital elevation model DEM to the longitude Long of the points Q, the latitude Lat of the points Q and the altitude h of the points Q. In a step S4, the calculator 101 calculates to which one of the primary burst images I i of the mosaic image MI belongs each primary projected point Q (p) . The calculator 101 determines the index i of the primary burst image Ii to which the primary projected point Q (p) belongs. This can be done by calculator 101 using primary cuts Di,i+1, Di-1,i which are determined by the calculator 101 between the primary burst images Ii and Ii+1 which are partially overlapping from one to another, as will be described below. In a step S5, the calculator 101 applies primary corrections C or C(Q (p) , i, θ) to the primary coordinates R, COL of the primary projected points Q (p) in the primary mosaic image MI and calculates primary corrected coordinates of the primary projected points Q (p) in the primary mosaic image MI. In a step S6, the calculator 101 calculates a primary warp function Ai (p) , which transforms the primary coordinates R, COL of the primary projected points Q (p) into the primary corrected coordinates of the primary projected points Q (p) in the primary mosaic image MI. This is done by the calculator 101 for each primary burst image Ii to which at least one Qi (p) of the primary projected points Q (p) belongs, as detemined by the step S4 by the calculator 101. In a step S7, the calculator 101 transforms each primary burst image Ii by the primary warp function A i (p) . The calculator 101 assembles into a corrected primary mosaic image MI2 the primary burst images I i having been transformed by the primary warp functions A i (p) . The calculator 101 may output the corrected primary mosaic image MI2 on the physical output interface 104. This assembling may be done as described below. In an embodiment of the invention, shown in Figures 3 and 4, steps S10, S11, S13, S15, S16 and S17 similar respectively to steps S0, S1, S3, S5, S6 and S7 may be carried out for the secondary burst images I’ i described above, in addition to steps S0, S1, S2, S3, S4, S5, S6 and S7 carried out for the primary burst images I i . In steps S10, S11, S13, S14, S15, S16 and S17, the features having the term « primary » are replaced by features having the term « secondary », such as : secondary set of secondary burst images I’i ; secondary metadata comprising a first secondary origin Bi (s) of each secondary burst image I’i, height h’i of each secondary burst image I’ i , width l’ i of each secondary burst image I’ i ; secondary mosaic image MI’; secondary coordinates R’, COL’ ; second secondary origin O (s) ; secondary projected points Q (s) ; secondary corrections C or C(Q (s) , i, θ); secondary corrected coordinates secondary warp function Ai (s) ; secondary corrected coordinates ; secondary corrected mosaic image MI2’. Of course, one or more secondary corrected mosaic image MI2’ may be calculated by the calculator 101 base on one or more sets of secondary burst images I’i, having been acquired at different times by the radar. In step S14, the calculator 101 calculates to which one of the secondary burst images I’ i of the secondary mosaic image MI’ belongs each secondary projected point Q (s) by finding the secondary burst image I’i that match the primary burst images Ii. The calculator 101 determines the index i of secondary burst image I’i to which the primary projected point Q (p) belongs. This can be done by calculator 101 using the primary cuts Di,i+1 , Di-1,i which are determined by the calculator 101 between the primary burst images Ii and Ii+1 which are partially overlapping from one to another, as will be described below, and by finding the secondary burst image I’ i that match the primary burst images I i . In an embodiment of the invention, the primary warp function Ai (p) is a primary affine resampling matrix for each primary burst image Ii. The secondary warp function Ai (s) may be a secondary affine resampling matrix for each secondary burst image I’i. In an embodiment of the invention, in step S6, the calculator 101 calculates the primary warp function A i (p) as the primary affine resampling matrix transforming a first primary difference Q (p) -Bi (p) between the primary coordinates R, COL of the primary projected points Q (p) and the first primary origin Bi (p) , into a second primary difference ^^ (^) -B i (p) between the primary corrected coordinates of the primary projected points and the first primary origin B i (p) of the This is done by the calculator 101 for each primary burst image I i to which at least one Q i (p) of the primary projected points Q (p) belongs, as detemined by the step S4 by the calculator 101. In step S16, the calculator 101 may calculate the secondary warp function Ai (s) as the secondary affine resampling matrix transforming a first primary difference Q (p) -Bi (p) between the primary coordinates R, COL of the primary projected points Q (p) and the first primary origin B i (p) of the secondary burst image Ii, into a second secondary difference ^^ (^) -Bi (s) between the secondary of the secondary projected points and the first secondary origin Bi (s) of the secondary burst image Ii, i. e. : This is done by the calculator 101 for each secondary burst image I’ i to which at least one Qi (s) of the secondary projected points Q (s) belongs, as detemined by the step S14 by the calculator 101. In another embodiment of the invention, the primary warp function Ai (p) is a primary homography resampling matrix for each primary burst image I i . The secondary warp function Ai (s) may be a secondary homography resampling matrix for each secondary burst image I’i). In another embodiment of the invention, the primary warp function A i (p is a primary polynomial function for each primary burst image I i . The secondary warp function Ai (s) may be a secondary polynomial function for each secondary burst image I’i. In an embodiment of the invention, the primary corrections C(Q (p) , i, θ) of the primary coordinates R, COL of the primary projected points Q (p) in the primary mosaic image MI vary depending on at least one primary parameter calculated from the primary metadata or contained in the primary metadata. Embodiments of primary parameters cosαinc; Rank, τmid, IW2, EW3, PRI ; cosαinc; Kr, Nlburst, Fa, t0,i, fηc,i(τi), ka,i(τi), vs,i, fc, kψ,i, Nsswath are described below. The secondary corrections C(Q (s) , i, θ) of the secondary coordinates R’, COL’ of the secondary projected points Q (s) in the secondary mosaic image MI’ may vary depending on at least one secondary parameter calculated from the secondary metadata or contained in the secondary metadata. Embodiments of secondary parameters cosα’ inc ; Rank’, τ’ mid , IW2, EW3, PRI’ ; K’ r , N’ lburst , F’ a , t’ 0,i , f’ ηc,i i ), k’ a,i i ), v’ s,i , f’ c , k’ ψ,i , Ns’ swath are described below. In an embodiment of the invention, in step S1, a burst azimuth time ti may be defined by the calculator 101 as being an affine function of the primary burst coordinate Ri of rows along the azimuth direction T in each burst image Ii, as follows : where ^ ^ is an azimuth sampling frequency contained in the primary metadata, t0,i is a burst reference azimuth time given by the primary metadata. In an embodiment of the invention, in step S1, a burst range time τ i may be defined by the calculator 101 as being an affine function of the primary burst coordinate COLi of columns along the range direction RD in each burst image Ii, as follows : where ^ ^ is a range sampling frequency contained in the primary metadata, τ 0,i is a burst reference range time given by the primary metadata. In the case of projection (specific case of geolocation), knowing the tridimensional coordinates of a point ^ = ( ^, ^, ^ ) , and the satellite position ^^^ ( ^ ) = ^^^^ ^ , ^^^ ^ , ^^^ ^ ^ with azimuth time, it is possible to find the azimuth time ^ ^ of the satellite when it "captured" the point (in many cases, it is simply the time of closest distance between the satellite and the point, i.e. ^ ^ = ^^^^^^ ^ ∥ ^ − ^^^(^) ∥^ ^ , which is a problem that may be solved used various optimization techniques). The distance between the satellite and the point can thus be computed as ^^^^ =∥ ^ − ^^^ ^ ^^ ^ ∥^ and converted to time ^^ using the speed of light in vacuum ^ : 2 × ^^^ ^ ^ ^ = ^ Knowing ^^ and ^^ , the calculator 101 calculates the primary burst coordinates Ri, COLi of each point Q in each burst image Ii from the above equations. The sampling frequencies ^ ^ , ^ ^ and samples of the satellite position as function of azimuth time (that can be used to fit a polynomial model) are given in the metadata as "rangeSamplingRate", "azimuthFrequency", "orbit" (more precise orbits can also be downloaded). ^ ^,^ , ^ ^,^ are given in the metadata for every burst image Ii (to be precise, the burst images Ii are distributed with a number of invalid samples in the border, the position of these samples are specified in the metadata and the burst image I i can be redefined at its valid boundaries). Using this, along with the burst size (also calculated from the metadata), the minimal continuous mosaic image MI containing all the burst image Ii is calculated along with the position of each burst image Ii in it. The calculator calculates the primary coordinates R and COL from the burst azimuth time ti, from theburst range time τi, from the primary burst coordinates COLi and R i , from the azimuth sampling frequency F a , from the range sampling frequency Fr, from the burst reference azimuth time t0,i and from the burst reference range time τ0,i. Thus, an azimuth time t may be defined by the calculator 101 as being an affine function of the primary coordinate R of rows along the azimuth direction T in the primary mosaic image MI, as follows : ^ ^ = + ^ ^ ^ ^ where ^ ^ is an azimuth sampling frequency contained in the primary metadata, wherein t 0 is a reference azimuth time calculated by the calculator 101 from the azimuth sampling frequency F a and from the burst reference azimuth time t 0,i . The range time τ may be defined by the calculator 101 as being an affine function of the primary coordinate COL of columns along the range direction RD in the primary mosaic image MI, as follows : ^^^ ^ = + ^ ^ ^ ^ where ^ ^ is a range sampling frequency contained in the primary metadata, wherein τ0 is a reference range time calculated by the calculator 101 from the range sampling frequency Fr and from the burst reference range time τ0,i. Thus, a SAR image’s reference system (origin) can be defined by the calculator 101 as times ^ ^ , ^ ^ , where ^ ^ is the reference azimuth time (time along the satellite flight path, i.e. along the azimuth direction T) and ^ ^ is the reference range time (time of flight transversally to the flight path, i.e. along the range direction RD), along the radar wave propagation direction. It is two- way because the wave is reflected on the ground, i.e. it is associated with twice the distance from satellite to ground). Then, knowing the azimuth sampling frequency F a and range sampling frequency F r , any primary coordinates R, COL of pixels P in the primary mosaic image MI can be converted by the calculator 101 to the azimuth time t and to the range time τ and vice-versa. Correction of an atmospheric path delay In an embodiment of the invention, in step S5, the calculator 101 calculates the primary corrections C(Q (p) , i, θ) of the primary coordinates R, COL of the primary projected points Q (p) in the primary mosaic image MI, which comprise a primary function δτAPD of correction of an atmospheric path delay between the synthetic aperture radar and the ground surface area. The calculator 101 may calculate the primary function δτ APD of correction of the atmospheric path delay as varying at least with the altitude h of the points Q corresponding to the primary projected points Q (p) and with a first primary parameter cosαinc calculated from the primary metadata (or from the geometry, by applying the cosine rule in the triangle between the center of the earth, the tridimensional point Q and the satellite position that corresponds to Q, ^^^^^ ^ ^). The calculator 101 may calculate the primary function δτAPD of correction of the atmospheric path delay as being proportional to a primary polynomial function ψ tropo,zenith of the altitude h of the points Q corresponding to the primary projected points Q (p) . The primary polynomial function ψtropo,zenith of the altitude h of the points Q may has a degree equal to 2. The calculator 101 may calculate the first primary parameter as cosα inc , wherein α inc is a primary incidence angle between a line of sight of the synthetic aperture radar and a local normal to the ground surface area. The calculator 101 may calculate the primary polynomial function ψtropo,zenith, expressed in meters for h expressed in meters, according to the following . The calculator 101 may calculate the primary function δτ APD of correction of the atmospheric path delay as When the electromagnetic wave emitted by the radar travels in the earth’s atmosphere, its velocity is actually lower than the speed of light in vacuum. Therefore, there is a travel delay due to the fact that the wave is slower than actually assumed. The atmospheric delay can be decomposed into two parts for radio signals of up to 30 GHz: the tropospheric delay of the order of 2.5-4 m, and the ionospheric delay that may reach up to 1m for C-band (5.405 GHz, Sentinel-1 satellite). The primary function δτ APD of correction of the atmospheric path delay enables to correct the tropospheric delay with a simple model. Using the altitude h, the delay in the zenith direction can be approximated by the primary polynomial function ψtropo,zenith. In an embodiment of the invention, in step S15, the calculator 101 calculates the secondary corrections C(Q (s) , i, θ) of the secondary coordinates R’, COL’ of the secondary projected points Q (s) in the secondary mosaic image MI’, which comprise a secondary function δτ’ APD of correction of an atmospheric path delay between the synthetic aperture radar and the ground surface area. The calculator 101 may calculate the secondary function δτ’APD of correction of the atmospheric path delay as varying at least with the altitude h of the points Q corresponding to the secondary projected points Q (s) and with a first secondary parameter cosα’ inc calculated from the secondary metadata (in the same manner as described above). The calculator 101 may calculate the secondary function δτ’APD of correction of the atmospheric path delay as being proportional to a secondary polynomial function ψ’tropo,zenith of the altitude h of the points Q corresponding to the secondary projected points Q (s) . The secondary polynomial function ψ’tropo,zenith of the altitude h of the points Q may has a degree equal to 2. The calculator 101 may calculate the first secondary parameter as cosα’ inc , wherein α’ inc is a secondary incidence angle between a line of sight of the synthetic aperture radar and a local normal to the ground surface area. The calculator 101 may calculate the secondary polynomial function ψ’tropo,zenith, expressed in meters for h expressed in meters, according to the following equation : The calculator 101 may calculate the secondary function δτ’ APD of correction of the atmospheric path delay as with ′ tropo, zenith ^^^^′ inc Bistatic correction of azimuth time In an embodiment of the invention, in step S5, the calculator 101 calculates the primary corrections C(Q (p) , i, θ) of the primary coordinates R, COL of the primary projected points Q (p) in the primary mosaic image MI, which comprise a primary function δtbistatic of bistatic correction of the azimuth time t of the synthetic aperture radar along the azimuth direction T. The moving satellite moves, between transmision of a pulse by the synthetic aperture radar and reception of the pulse by the synthetic aperture radar along the azimuth direction T. The calculator 101 may calculate the primary function δtbistatic of bistatic correction of azimuth time t of the synthetic aperture radar as varying in an affine manner with the primary column number COL of the primary pixels P in the primary mosaic image MI and with primary second parameters Rank, PRI given in the primary metadata. PRI is a Pulse Repetition Interval, Rank is reception rank between reception and emission of each pulse. The calculator 101 may calculate the primary function δtbistatic of bistatic correction of azimuth time t of the synthetic aperture radar as: ^ ^ where ^ is obtained by the above equations, ^ ^^^ is the range time τ computed at the middle column of the middle subswath (IW2 in IW mode or EW3 in EW mode). When forming each SAR burst image Ii, a pulse is transmitted then echos are received. The satellite does not stop at a fixed position in orbit and wait for the pulse to reflect on the ground until reception (stop-and-go approximation). The satellite will keep moving and transmitting pulses each Pulse Repetition Interval (PRI). Within a PRI, after each transmission window, there is a reception window. A pulse transmitted in the ^ ^^ PRI is received in the ( ^ + ^^^^ )^^ PRI received window. The monodirectional echo packets need to be reorganized as a bidimensional image, by associating to each echo the azimuth time t. However, this is not as straightforward as the stop-and-go approximation. The time elapsed between transmission and reception needs to be taken into account. Ideally, the time that should be associated to the echo is in the middle between transmission and reception. However, this is not the case for the time associated to the echo in Sentinel- 1 satelite. A correction was performed by the Sentinel-1 Instrument Processing Facility (IPF), but it is a rough correction that does not achieve the highest possible accuracy. The formula for the most accurate bistatic correction in Sentinel-1 satellite is : where ^ is obtained by the above equations, ^ ^^^ is the range time τ computed at the middle column of IW2 or EW3 (i.e. the column at the center of the whole swath). ^ ^^^ is used to reverse the IPF correction. In an embodiment of the invention, in step S15, the calculator 101 calculates the secondary corrections C(Q (s) , i, θ) of the secondary coordinates R’, COL’ of the secondary projected points Q (s) in the secondary mosaic image MI, which comprise a secondary function δt’bistatic of bistatic correction of the azimuth time t’ of the synthetic aperture radar along the azimuth direction T. The moving satellite moves, between transmission of a pulse by the synthetic aperture radar and reception of the pulse by the synthetic aperture radar along the azimuth direction T. The calculator 101 may calculate the secondary function δt’bistatic of bistatic correction of azimuth time t’ of the synthetic aperture radar as varying in an affine manner with the secondary column number COL’ of the secondary pixels P in the secondary mosaic image MI’ and with secondary second parameters Rank’, PRI’ given in the secondary metadata. PRI’ is a Pulse Repetition Interval, Rank’ is reception rank between reception and emission of each pulse. The calculator 101 may calculate the primary function δt’bistatic of bistatic correction of azimuth time t’ of the synthetic aperture radar as: where ^′ is obtained by the above equations, ^′ ^^^ is the range time τ computed at the middle column of the middle subswath (IW2 in IW mode or EW3 in EW mode). Intra pulse motion correction In an embodiment of the invention, in step S5, the calculator 101 calculates the primary corrections C(Q (p) , i, θ) of the primary coordinates (R, COL) of the primary projected points Q (p) in the primary mosaic image MI as comprising a primary function δτ intra_pulse,i of intra pulse motion correction. The radar provides Each burst image Ii is affected differently from the other by the intra pulse shift. The primary function δτintra_pulse,i of intra pulse motion correction corrects these shifts. In an embodiment of the invention, in step S5, the calculator 101 calculates the primary function δτ intra_pulse,i of intra pulse motion correction as being negatively proportional to a primary Doppler centroid frequency f DC,i i ,t i ). The calculator 101 calculates the primary Doppler centroid frequency fDC,i(τi,ti) depending on primary third parameters calculated from the primary metadata or given in the primary metadata. The primary third parameters may comprise Kr, Nlburst, Fa, t0,i, fηc,i(τi), ka,i(τi), vs,i, fc, kψ,i, Nsswath as described below. The calculator 101 calculates the primary Doppler centroid frequency fDC,i(τi,ti) as depending in an affine manner on the primary row number Ri of the primary pixels P in the primary burst image Ii. The calculator 101 calculates at least one (such as fηc,i(τi) and ka,i(τi) as described below) of the primary third parameters as depending on the primary column number COL i of the primary pixels P in the primary mosaic image MI. In an embodiment of the invention, in step S5, the calculator 101 calculates the primary function δτ intra_pulse,i of intra pulse motion correction as where ^ ^ is a (positive, not null) range chirp modulation rate given in the primary metadata (for example as "txPulseRampRate"). This correction can be modeled as a removal of a range shift dependent on the Doppler centroid frequency, where the primary function δτ intra_pulse,i of intra pulse motion correction is the two-way range time shift present in the data in seconds, ^ ^^,^ (^ ^ , ^ ^ ) is the Doppler centroid (DC) in Hz and ^ ^ is the chirp rate in Hz/second. ^ ^ is constant on the swath, while ^ ^^,^ ( ^ ^ , ^ ^ ) varies linearly inside a burst image Ii with respect to the burst azimuth time ti referenced to the middle of the burst. The Doppler centroid frequency ^ ^^,^ ( ^ ^ , ^ ^ ) comprises the geometric Doppler and the contributions of the radar antenna steering. The Doppler frequency of a target varies along track (azimuth) even within a burst image I i , whereas the burst-level processor assumes a single frequency for the duration of the burst image Ii. As a result, a range bias will be introduced for a given target that depends on its along-track offset from the burst image center. For IW products, this effect may cause azimuth geolocation errors up to 0.5 m in either direction. The primary function δτintra_pulse,i of intra pulse motion correction enables to accurately predict the position of Q within the mosaic image MI. In an embodiment of the invention, in step S5, the calculator 101 calculates the primary Doppler centroid frequency f DC i ,t i ) according to the following equation: where The burst azimuth time t i is defined by the calculator 101 as being an affine function of the primary burst coordinate Ri of rows along the azimuth direction T in the primary burst image Ii, as indicated above. The calculator 101 calculates ^ ^^^,^ according to the following equation: where ^ ^^^^^^ is the number of rows per image burst Ii and is given in the primary metadata (for example as "linesPerBurst"), where ^ ^ is the azimuth sampling frequency contained in the primary metadata, where t0,i is the burst reference azimuth time given by the primary metadata. ^ ^ ( ^ ^ ) is therefore the burst azimuth time ti referenced to the middle of the image burst Ii. is given in the primary metadata. The calculator 101 calculates according to the following equation: where ^ ^,^ ( ^ ^ ) is the azimuth frequence modulation rate given in the primary metadata and depends on the burst range time τi. are given in the metadata as sequence of range polynomials regularly updated with the burst azimuth time ti. For the burst image Ii, it is recommended to use the closest polynomial ^ ^^^,^ of the burst image Ii. The calculator 101 calculates ^ ^,^ according to the following equation: where ^ ^,^ is the satellite velocity at the middle of the burst image Ii ^^^^ ^ ^,^ = ^ ^^ ^ ^^^,^ ^ ^ is the carrier frequency given in the primary metadata (for example as "radarFrequency"), ^ ^,^ is the steering rate given in the primary metadata (for example as "azimuthSteeringRate"). The calculator 101 calculates the satellite velocity ^ ^,^ from the primary metadata. The calculator 101 calculates ^ ref,i ( ^ ) depending on the burst range time τ i according to the following equations: with where ^^ swath is the number of range samples in the swath given in the primary metadata (for example as "samplesPerBurst"), where ^ ^,^ ^ ^^ swath ^ ^ is ^ ^,^ ( ^ ^ ) calculated by the calculator 101 for the burst range time τi corresponding to half of ^^ swath . Of course, all the 3 corrections : primary function δτAPD of correction of an atmospheric path delay, primary function δtbistatic of bistatic correction and primary function δτintra_pulse,i of intra pulse motion correction mentionned above may be made during step S5 and/or all the 3 corrections : secondary function δτ’ APD of correction of an atmospheric path delay, secondary function δt’ bistatic of bistatic correction and secondary function δτ’ intra_pulse,i of intra pulse motion correction mentionned above may be made during step S15. Of course, one or more secondary burst images may be processed by the method. Of course, the aspects, embodiments, features, possibilities, steps, operations and examples of the disclosure mentioned above may be combined one with another or may be selected independently one from another.