Satellite-image-based surface reconstruction in `MicMac`
A tutorial
- Goal
- Projet set-up
- 1. Extract SIFT tie-points
- 2. RPC-bundle adjustment
- 3. Surface reconstruction
Goal
In this tutorial we will introduce you to satellite image processing in MicMac. The goal of the exercise is to compute the surface and generate an orthohoto given a set of modern satellite images and their RPC geolocalisations. After setting-up MicMac and downloading the dataset, the pipeline is as follows:
- Tie-points extraction
- RPC-bundle adjustement
-
Surface computation
3.1. Method1: Matching in object geometry
3.3. Method2: Matching in image geometry and fusion
Contact: ewelina.rupnik(at)ign.fr
import os
from os.path import exists, join, basename, splitext
import numpy as np
import cv2
import matplotlib.pyplot as plt
Dependencies_install = True
MicMac_clone = True
MicMac_cmake = True
MicMac_build = True
YOUR_PATH = '/content/'#MyDrive/micmac/satellites/'
!cd $YOUR_PATH
!pwd
if Dependencies_install:
!apt update
!apt install -y cmake
!pip install dlib
!apt-get install imagemagick proj-bin exiv2
!pip install wget gdown
if MicMac_clone:
if not exists(YOUR_PATH+'micmac/'):
git_repo_url = 'https://github.com/micmacIGN/micmac.git'
!git clone $git_repo_url
if MicMac_cmake:
!cd micmac
if not exists(YOUR_PATH+'micmac/build'):
!mkdir $YOUR_PATH"micmac/build"
!cd $YOUR_PATH"micmac/build"
!cmake $YOUR_PATH"micmac" -DBUILD_POISSON=OFF
if MicMac_build:
!make install -j28
import os
os.environ['PATH'] += ":/content/micmac/bin/"
!echo $PATH
# if you can see the commands printed to the screen, everything is OK
!mm3d
dataset_url = 'https://drive.google.com/uc?id=18hmQL5kIqhcnR5ahp8IUsMZxLv7jvjgB'
!gdown $dataset_url -O "satellite_data.tar.gz"
# unpack
if not exists(YOUR_PATH+'satellite_data'):
!mkdir $YOUR_PATH'satellite_data'
!tar -xf satellite_data.tar.gz -C $YOUR_PATH'satellite_data'
%cd $YOUR_PATH'satellite_data'
# utility functions to visualise tie-points
utils_url='https://drive.google.com/uc?id=1ATO1Nz_aXApxVnm6l7x1xappGXtcjuvp'
!gdown $utils_url -O "mm3d_utils.py"
1. Extract SIFT tie-points
-
computation strategy: there exist several predefined strategies to compute tie-points:
Line,All,MulScale,File. We will use theAllstrategy where tie-points are searched between all possible pairs. Refer to MicMac documentation for the other modes. -
image resolution: tie-points extraction is very costly, and to limit the computation time we usually downsampled the images; in this example, indicate resolution of
-1which means full-resolution images; otherwise, if set to, e.g.,2000, the images will be downsampled such that the larger image dimension (typically the width) will have2000pixels; the other dimension will have a size that is proportionally smaller; - ExpTxt=1: the extracted tie-points will be saved in a text format (as opossed to the default dat format)
-
results: tie-points are stored in the
Homoldirectory. For instance, tie-points correponding to imageIm1.tifwill be stored inHomol/PastisIm1.tif/directory. IfIm1.tifoverlaps withIm2.tifandIm3.tif, their tie-points will be stored inHomol/PastisIm1.tif/Im2.tif.datandHomol/PastisIm1.tif/Im3.tif.dat, respectively. If you chose to export in the text format, thedatextension will be replaced withtxt.
Note: Intermediary results are stored in the Pastis directory. It takesa significant amount of space and is not used at later processing stages, therefore you may delete it.
!mm3d Tapioca All .*tif -1 ExpTxt=1 @ExitOnBrkp
import mm3d_utils
aIm1 = cv2.imread('TPMM_0435.tif',cv2.IMREAD_IGNORE_ORIENTATION)
aIm2 = cv2.imread('TPMM_0566.tif',cv2.IMREAD_IGNORE_ORIENTATION)
TPtsVec = mm3d_utils.ImportHom("Homol/PastisTPMM_0435.tif/TPMM_0566.tif.txt")
mm3d_utils.plot_images([np.asarray(aIm1),np.asarray(aIm2)])
mm3d_utils.plot_tiepts2(np.asarray(TPtsVec,dtype=float))
Read the RPCs in DIMAP format
This function reads the DIMAP format RPCs and converts it to a MicMac format. Several parameters are specified here:
-
(.*).tifthis is the pattern of input images (note the dot preceding the star which is the posix convention) -
\$1.xmlis the corresponding pattern of RPC files; I use here a regular expression that associates the image name with its corresponding RPC file name; you may also run the command independently for each image if you're not familiar with regular expressions; -
RPC-d0is the directory name where the converted files will be stored; it will serve as input in the following step, i.e., the bundle adjustment; -
Degre=0, the degree of the polynomial correction;By choosing a zero-degree polynomial we will correct the satellite's geolocalisation by modelling a 3D image shift; please refer to [Rupnik et al., 2016] for more on the method.
-
ChSys=WGS84toUTM.xmldefinition of the projection coordinate sytem; MicMac expects that the processing coordinate frame is euclidean and all three coordinates have the same unit. The RPCs are expressed in geographical coordinates which are neither euclidean, nor unique in terms of units. To overcome that, MicMac will transfer, on the fly, the RPCs to a user-defined coordinate system, in this exemple defined in theWGS84toUTM.xmlfile. The definition of the coordinate system follows theproj4library convention. You can retrieve the code corresponding to the coordinate frame of your interest fromhttps://spatialreference.org/
Rupnik, E., Deseilligny, M.P., Delorme, A. and Klinger, Y., 2016. Refined satellite image orientation in the free open-source photogrammetric tools Apero/Micmac. ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 3, p.83.</p>
</div>
</div>
</div>
Within the bundle adjustemnt, MicMac will estimate the parameters of $D_x$ and $D_y$ functions, that are the polynomials you have defined in $\begin{equation}
\begin{split}
y = g (\phi, \lambda, h) + {D_y (x,y)}\\
x = h (\phi, \lambda, h) + {D_x (x,y)}
\end{split}
\end{equation}$ where $g$ and $h$ are the RPC functions, and $\begin{equation}
\begin{split}
D_y (x,y) = \sum_{i=0}^m \sum_{j=0}^n a_{ij} \cdot x^i y^j \\
D_x (x,y) = \sum_{i=0}^m \sum_{j=0}^n b_{ij} \cdot x^i y^j
\end{split}
\end{equation}$ The input parameters: One way to asses the quality of the adjustment is to look at the tie-points residual (for more sophisticated quality estimates see The bundle adjustment is carried out in several iterations. Let's look at image We will now calculate the surface with the semi-global dense image matching [Deseilligny \& Paparoditis, 2006]. Deseilligny, M. and Paparoditis, N., 2006. A multiresolution and optimization-based image matching approach: An application to surface reconstruction from SPOT5-HRS stereo imagery. Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences, 36(1/W41), pp.1-5.</p>
</div>
</div>
</div>
The computation will be carried out in the so-called terrain geometry, where the optimization is defined in the The input parameters are: The matching is carried out at multi-resolutions, i.e., we first calculate the surface using low resolution images (top-most level of the image pyramid), then we propagate the solution to lower levels and refine it, so long we have not reached the bottom of the image pyramid. The surface reconstructions at each level are stored inside the The computation will be carried out in the image geometry, where the optimization is defined in the The multiview pipeline is as follows [Rupnik et al., 2018]: Rupnik, E., Pierrot-Deseilligny, M. and Delorme, A., 2018. 3D reconstruction from multi-view VHR-satellite images in MicMac. ISPRS Journal of Photogrammetry and Remote Sensing, 139, pp.201-211.</p>
</div>
</div>
</div>
We will compute two surfaces using two different subsets of the images. Input parameters: The $1^{st}$ command The $2^{nd}$ command The fusion takes all the surfaces specified by the regular expression and merges it. The fusion takes into account the correlation images and treats it as confidence maps. Input parameters: Tha result is saved to !mm3d Convert2GenBundle -help @ExitOnBrkp
!mm3d Convert2GenBundle "(.*).tif" "\$1.xml" RPC-d0 ChSys=WGS84toUTM.xml Degre=0 @ExitOnBrkp
Run the adjustment
Convert2GenBundle method. In this example our observations are the tie-points. It is also possible to include ground control points.
RPC-d0 is the folder with the initial geolocalisationRPC-d0_adj is the folder where the adjusted geoloc is savedExpTxt=1 indicates that tie-points are stored in text format !mm3d Campari ".*tif" RPC-d0 RPC-d0-adj ExpTxt=1 @ExitOnBrkp
Interpreting the results
MMTestOrient in MicMac documentation).TPMM_0435.tif in the last iteration:
RES:[TPMM_0435.tif][g] ER2 0.24636 Nn 100 Of 11753 Mul 5171 Mul-NN 5171 Time 1.15821
0.24636 pixels is the mean residual calculated over all tie-points (i.e., $\sigma$ of the bundle)
Nn 100 means that 100$\%$ of tie-points were considered as inliers11753 there were as many tie-points found5171 there were as many multiple tie-points found (out of the 11753), i.e., tie-points observed in at least 3 images;
3. Surface reconstruction
Method1: Matching in object geometry
(X,Y) of the object space, and we are looking for the most optimal Z-coordinates (see Figure below). This geometry is well adapted to 2.5D surface computation.
UrbanMNE is a predefined term and it defines a number of processing parameters (e.g., low regularization, small matching wodows, terrain geometry).*tif is the image set that will be used in the processingRPC-d0-adj is the name of the directory containing the geoloclisationSzW=2 defines the matching window size, i.e., with size set to 2, the window size is 5x5Regul is the regularization term $\alpha$; in UrbanMNE it is by default set to 0.02 becase in urban zones we're typically interested in reconstructing fine details; the dataset used in this example, however, represents a smooth surface so we're just fine with a more agressive regularization; moreover, our images are quite noisy and by adding more regularization we will avoid noisy surface reconstructions;DoOrtho=1, this parameters wil force MicMac to create individual orthomosaic, i.e., rectify each image; the rectified images are stored in Ortho-MEC-Malt/Ort_*.tif; to create the final orthophotomap we will still need to do mosaicing with Tawny (later in this tutorial);NbVI=2 sets the necessary minimum number of images for MicMac to compute the surface; by default the value is set to 3 which means that in areas with only two overlapping images, the surface will not be computed;EZA=1, this parameter will force the output surface raster to save the Z-coordinates in their absolute values; without explicitly forcing MicMac to do that, to avoid having to store large values inside the raster tiff, it will apply a normalisation (normalisation parameters OrigineAlti and ResolutionAlti are to be found in the accompanying xml files);!mm3d Malt UrbanMNE ".*tif" RPC-d0-adj SzW=2 Regul=0.2 DoOrtho=1 NbVI=2 EZA=1 @ExitOnBrkp
Reading the output files
MEC-Malt directory. Here's how to decipher the files:
MEC-Malt/Z_Num8_DeZoom1_STD-MALT.tif represents the surface raster at the highest resolution; MEC-Malt/Z_Num8_DeZoom1_STD-MALT.xml is its metadata file that encodes the georeferencing; let's assume you'd like to convert a pixel $(i,j)$ from its image coordinates to its georeferenced coordinates (i.e., object coordinates):
Z-coordinate normalised :$Z^{img}_{i,j} = Z\_Num\_File^{img}(i,j)$, > Z-coordinate :$Z^{obj}_{i,j} = OrigineAlti + ResolutionAlti \cdot Z^{img}_{i,j}$> XY-coordinates :$(X,Y) = OriginePlani + ResolutionAlti \cdot (i,j)$MEC-Malt/Masq_STD-MALT_DeZoomX.tif is a binary mask file that is a result of your input mask (if you used one) and a mask that is automatically calculated in the matching optimisation phase;MEC-Malt/Correl_STD-MALT_Num_X.tif are the images with storing the correlation scores (it is not pure correlation, it is the correlation store after the aggregation step)!mm3d GrShade MEC-Malt/Z_Num8_DeZoom1_STD-MALT.tif ModeOmbre=IgnE Mask=MEC-Malt/Masq_STD-MALT_DeZoom1.tif @ExitOnBrkp
surface_shade_im = cv2.imread("MEC-Malt/Z_Num8_DeZoom1_STD-MALTShade.tif",cv2.IMREAD_IGNORE_ORIENTATION)
fig, ax = plt.subplots(figsize=(30, 10))
ax.imshow(surface_shade_im,cmap="gray")
plt.tight_layout()
Generate an orthophotmap
Tawny will mosaic the per-image orthopĥotomosaics created in Malt (i.e., Ortho-MEC-Malt/Ort_*.tif), during matching. It will additionally perform some basic radiometry equalization. The output orthoimage is stored in Ortho-MEC-Malt/Orthophotomosaic.tif, and its georeferencing is stored in Ortho-MEC-Malt/Orthophotomosaic.twf.!mm3d Tawny Ortho-MEC-Malt/ @ExitOnBrkp
ortho_im = cv2.imread("Ortho-MEC-Malt/Orthophotomosaic.tif",cv2.IMREAD_IGNORE_ORIENTATION)
fig, ax = plt.subplots(figsize=(30, 10))
ax.imshow(ortho_im,cmap='gray')
plt.tight_layout()
Method2: Multiview matching in image geometry and fusion
(x,y) of the image space, and we are looking for the most optimal depths (see Figure below). This geometry is well adapted to true 3D surface reconstruction. Because the individual reconstructions are computed in image coordinate frames, a fusion will but carried out at the end.
Malt GeomImage)NuageBascule)SMDM)
Do two per-triplet image matching
GeomImage is a predefined term and it defines a number of processing parameters (e.g., low regularization, small matching wodows, image geometry)TPMM_(0435|0566|1088).*tif is the image set that will be used in the processing; in this example it is an image triplet; we take 3 consecutive images to make sure that the $\frac{B}{H}$ ratios in the set are relatively smallRPC-d0-adj is the name of the directory containing the geolocalisationMaster=TPMM_0566.tif is the master image, i.e., the optimization is defined over each pixel of that imageSzW=2 and Regul=0.2, similarily to Method1, we add regularization and use bigger correlation windows because: (1) the geometry of the surface is smooth (i.e., no discontinuities), and (2) the images are quiet noisy;NbVI=2 sets the necessary minimum number of images for MicMac to compute the surface; by default the value is set to 3 which means that in areas with only two overlapping images, the surface will not be computed;!mm3d Malt GeomImage "TPMM_(0435|0566|1088).*tif" RPC-d0-adj Master=TPMM_0566.tif SzW=1 Regul=0.1 NbVI=2 ZPas=1 @ExitOnBrkp
!mm3d Malt GeomImage "TPMM_(0566|1088|1216).*tif" RPC-d0-adj Master=TPMM_1088.tif SzW=1 Regul=0.1 NbVI=2 ZPas=1 @ExitOnBrkp
Transform the depth maps to a common coordinate frame
Malt UrbanMNE does image matching in ground geometry (we've used it before), but since we indicate DoMEC=0 it will not calculate the matching, all it will do is to create metadata defining the coordinate frame of the ground geometry. The output, as before, is stored in MEC-Malt folder;NuageBascule will apply the transformation; the parameters are:
MM-Malt-Img-TPMM_0566/NuageImProf_STD-MALT_Etape_8.xml is the metadata file defining the input coordinate frameMEC-Malt/NuageImProf_STD-MALT_Etape_8.xml is the metadata file defining the target coordinate frame (i.e., it is the terrain geometry)Fusion/DSM_Tri1.xml is the output metadata file, i.e., the input file transformed to the target coordinate frame; it will be accompagned by several other files containing the surface itself, the mask and the correlation image (see the inside of the Fusion folder);!mm3d Malt UrbanMNE ".*tif" RPC-d0-adj DoMEC=0 @ExitOnBrkp
# create a directory that will store the fused surface
!mkdir Fusion
# do 3D spatial similarity of the first triplet depth map
!mm3d NuageBascule MM-Malt-Img-TPMM_0566/NuageImProf_STD-MALT_Etape_8.xml MEC-Malt/NuageImProf_STD-MALT_Etape_8.xml Fusion/DSM_Tri1.xml @ExitOnBrkp
# do 3D spatial similarity of the second triplet depth map
!mm3d NuageBascule MM-Malt-Img-TPMM_1088/NuageImProf_STD-MALT_Etape_8.xml MEC-Malt/NuageImProf_STD-MALT_Etape_8.xml Fusion/DSM_Tri2.xml @ExitOnBrkp
Fuse the individual depth maps
Fusion/DSM_Tri.*xml the subset of surfaces that will be merged; Fusion/Fusion_Prof.tif, there is a corresponding mask and a correlation map named with _Mask and Correl postfixes, respectively.!mm3d SMDM Fusion/DSM_Tri.*xml @ExitOnBrkp
!mm3d GrShade Fusion/Fusion_Prof.tif Out=Fusion/Fusion_GShade.tif ModeOmbre=IgnE @ExitOnBrkp
surface_fused_shade_im = cv2.imread("Fusion/Fusion_GShade.tif",cv2.IMREAD_IGNORE_ORIENTATION)
fig, ax = plt.subplots(1,2,figsize=(15, 15))
ax[0].imshow(surface_fused_shade_im,cmap="gray")
ax[1].imshow(surface_shade_im,cmap="gray")
plt.tight_layout()
# export to ply
#!mm3d Nuage2Ply Fusion/Fusion.xml Out=Fusion.ply
