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:

  1. Tie-points extraction
  2. RPC-bundle adjustement
  3. Surface computation

    3.1. Method1: Matching in object geometry

    3.3. Method2: Matching in image geometry and fusion


Contact: ewelina.rupnik(at)ign.fr

Projet set-up

Download and compile MicMac & dependencies

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

Add environmental variable & download the dataset

The dataset consists of:

  • 4 images (tif)
  • 4 corresponding RPCs (xml)
  • WGS84toUTM.xml with the definition of a projection coordinate system (proj4 format)
import os
os.environ['PATH'] += ":/content/micmac/bin/"
!echo $PATH

# if you can see the commands printed to the screen, everything is OK
!mm3d
/usr/local/nvidia/bin:/usr/local/cuda/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/tools/node/bin:/tools/google-cloud-sdk/bin:/opt/bin:/content/micmac/bin/

************************************************************************ 
**                                                                    ** 
**    MicMac: a  free open source project  for photogrammetry         ** 
**     hosted at Ecole Nationale des Sciences Geographiques           ** 
**               in Marne-la-Vallee, for IGN-France                   ** 
**                                                                    ** 
**                                                                    ** 
**  The project is funded by :                                        ** 
**                                                                    ** 
**  - Institut National de l'Information Geographique et Forestiere   ** 
**    (IGN main funder since 2003)                                    ** 
**  - the french FUI Project "Culture 3D Cloud" (and MAP-CNRS)        ** 
**  - the french ANR Project "MONUMENTUM"(collaborating with MAP-CNRS)** 
**                                                                    ** 
**  Research also currently supported by :                            ** 
**  - CNES (French Space Agency) - via TOSCA Committee (and IPGP)     ** 
**  - Compagnie Nationale du Rhone                                    ** 
**  - Vinci-Construction-Terrassement                                 ** 
**  - ERC Advanced Grant A.Kaeaeb "ICEMASS" (University of Oslo)      ** 
**                                                                    ** 
**                                                                    ** 
**  Current Team: MP Deseilligny, D Jouin, J Belvaux, G Maillet,      ** 
**    L Girod, E Rupnik, JM Muller, M Daakir, TG Nguyen               ** 
**                                                                    ** 
**    Contact for participating : Marc.Pierrot-Deseilligny@ensg.eu    ** 
**                                                                    ** 
**    Hope you enjoy, todo list in case of any problem using MicMac : ** 
**      (0) Take a Pastis                                             ** 
**      (1) Switch your computer off and on again                     ** 
**      (2) Install it on Gnu-Linux (work better on)                  ** 
**      (3) See the wiki at http://micmac.ensg.eu/                    ** 
**      (4) Contact the forum http://forum-micmac.forumprod.com/      ** 
**                                                                    ** 
************************************************************************ 
mm3d : Allowed commands 
 AimeApprent	Stat on homologous point using orientation of 3D Model
 AllDev	Force development of all tif/xif file
 AlphaGet27	 Tool for relative positioning of objects on images
 AnalyseTrajStereopolis	Analyse trajectory of Stereopolis-like acquisition
 Ann	 matches points of interest of two images
 AperiCloud	 Visualization of camera in ply file
 Apero	 Compute external and internal orientations
 Apero2Meshlab	Convert Orientation from Apero-Micmac workflow to a meshlab-compatible format
 Apero2NVM	Matthieu Moneyrond's convertor to VSfM, MVE, SURE, MeshRecon 
 Apero2PMVS	 Convert Orientation from Apero-Micmac workflow to PMVS format
 AperoChImSecMM	 Select secondary images for MicMac
 ApplyHomog	Homographie application on images 
 Arsenic	 IN DEV : Radiometric equalization from tie points
 Aspro	 Init  External orientation of calibrated camera from GCP 
 BAR	Bascule robutse 
 Bascule	 Generate orientations coherent with some physical information on the scene
 BatchFDC	 Tool for batching a set of commands
 Blinis	Block Initialisation 
 C3DC	Automatic Matching from Culture 3D Cloud project
 CalcImScale	Calculate scale of image
 CalcMapAnalytik	Compute map2d between images using various model 
 CalcMapOfT	Compute value of map evol for a given T 
 CalcMapXYT	Compute map2d evol of T 
 CalibFinale	 Compute Final Radial distortion model
 CalibInit	 Compute Initial Radial distortion model
 Campari	 Interface to Apero, for compensation of heterogeneous measures
 CASA	 Analytic Surface Estimation
 CatImSaisie	 Do some stuff
 CenterBascule	 Relative to absolute using embedded GPS
 ChantierClip	 Clip Chantier
 CheckDependencies	 check dependencies to third-party tools
 CheckGCPStereopolis	Check GCP with strategy optimized for Stereopolis-like acquisition
 ChgSysCo	 Change coordinate system of orientation
 CleanPatByOri	Clean a pattern of image by Ori-XXX folder
 ClipIm	 Clip Chantier
 CmpCalib	 Compare two  calibrations
 CmpDenseMap	comparison of dense map 
 CmpIm	 Basic tool for images comparison
 CmpOri	 Compare two sets of orientation
 cod	 Do some stuff
 CoherEpip	 Test coherence between conjugate epipolar depth-map
 Compens	 Do some stuff
 ContrastFilter	Some contrast filtering 
 Convert2GenBundle	Import RPC or other to MicMac format, for adjustment, matching ...
 ConvertCalib	 Conversion of calibration from one model 2 the other
 ConvertIm	 Tool for convertion inside tiff-format
 ConvertOriCalib	Convert external orientation with new internal orientation 
 ConvertPolygone	 Do some stuff
 CoronaCorrect	 Correct some geometric defautl of corona images
 CreateEpip	 Create epipolar images
 DebugAI4GeoMasq	 For debuging masq problem appeard with AI4Geo
 Dequant	 Tool for dequantifying an image
 Devlop	 Do some stuff
 Digeo	 In development- Will compute tie points 
 DistPxFrac	Compute distribution of fractional part of paralax 
 DIV	Videos development (require ffmpeg)
 DivFilter	 Test divers filters
 DMatch2Hom	Dense matching 2 homologues 
 Donuts	Cyl to Torus (Donuts like)
 DroneFootPrint	Draw footprint from image + orientation (drone) in PLY and QGIS format
 Drunk	 Images distortion removing tool
 EditSet	Edition creation of a set of images/files
 ElDcraw	 Do some stuff
 EstimHomog	Homographie estimation from GCPs and image measurements 
 ExtractAppui3D	 Extract points from a 3D appui points xml file
 ExtractMesure2D	 Extract points from a 2D measures xml file
 ExtractRaw	Convert raw image with XML descriptor to tiff 
 ExtractStdRaw	Convert raw image with predefined XML descriptor (in XML_MicMac/DataBaseCameraRaw) to tiff 
 FermDenseMap	Consistency by closing on dense maps 
 FFTKugelhupf	 Version of Kugelhupf using FFT, expecetd faster when it works (if ever ...)
 FieldDep3d	 To export results of matching as 3D shifting
 FitsMatch	Test Match Images NewPHom 
 FusionDepl	Fusion carte de deplacement 
 GCP2MeasuresL2D	 Convert a set of GCP in measure of 2D lines using convention NameLine_x with x={1,2}
 GCPBascule	 Relative to absolute using GCP
 GCPConvert	Convert GCP from Txt 2 XML
 GCPCtrl	 Control accuracy with GCP
 GCPMerge	 Merging of different GCP files
 GCPVisib	 Print a list of GCP visibility in images
 GenCode	 Do some stuff
 Genepi	 Generate 3D/2d synthetical points from orientation
 GenerateBorderCam	 Generate the polygone of image contour undistorded
 Genere_Header_TiffFile	 Generate Header for internal tiling format 
 genmail	 Do some stuff
 GenPairsFile	Generate pairs files between one image and a pattern
 GenPrime	Generate prime 
 GenXML2Cpp	 Do some stuff
 GrapheHom	Compute XML-Visibility graph from approximate orientation
 GrapheStereopolis	Compute Pair of Image for Stereopolis
 Gri2Bin	 Do some stuff
 GrShade	 Compute shading from depth image
 HackToF	Hack ToF format 
 Help	Help on existing MicMac commands 
 Homol2GND	 Creates fake ground points for aerotriangulation wedge
 HomolFilterMasq	 Tool for filter homologous points according to masq
 HomolMergePDVUnik	 Tool for merge homologous point from unik point of view
 HomProfPx	 Export pixel correspondences from Px1 et Px2
 Im2XYZ	 tool to transform a 2D point (text file) to their 3D cloud homologous
 ImMire	 For generation of some synthetic calibration image
 Impaint	Basic Impainting
 ImRandGray	 Generate Random Gray Textured Images
 Init11P	 Init Internal & External from GCP using 11-parameters algo
 InitOriLinear	 Initialize orientation for linear acquisition
 InvHomolHomog	Homographie application on images 
 Kugelhupf	 Semi-automatic fiducial points determination
 L2L	 Project row/column in one image to another
 L3D2Ply	 Convert a set of 3D lines in space to a ply file
 Liquor	Orientation specialized for linear acquisition
 LumRas	 Compute image mixing with raking light
 Luxor	Orientation specialized for linear acquisition using a sliding window
 MakeGrid	 Generate orientations in a grid format
 Malt	 Simplified matching (interface to MicMac)
 MapCmd	 Transforms a command working on a single file in a command working on a set of files
 Martini	 New orientation initialisation (uncomplete, still in dev...) 
 MartiniGin	 New orientation initialisation (uncomplete, still in dev...) 
 MasqMaker	 Create Mask form image values
 MeasuresL2D2L3D	 Convert a set of images measures of 2D lines to 3D lines in space
 MergeDepthMap	 Merging of individual, stackable, depth maps 
 MergeHomol	 Merge Homol dir
 MergePly	 Merge ply files
 MergeSOMAF	 Tool for merging SetOfMesureAppuisFlottants XMLs
 MeshProjOnImg	 Reproject mesh on image
 MICMAC	 Computes image matching from oriented images
 MICMACSaisieLiaisons	 Low level version of SEL, not recommended
 MM1P	 Matching One Pair of images
 MM2DPosSism	 Simplified interface for post 2D post sismic deformation 
 MMAI4Geo	 Basic Matching for AI4Geo Satellite
 MMByP	 Matching By Pair of images
 MMCalcSzWCor	 Compute Image of Size of correlation windows (Atomic tool, for adaptive window in geom image)
 MMHomCorOri	 Tool to compute homologues for correcting orientation in epip matching
 MMInitialModel	 Initial Model for MicMac 
 MMMergeCloud	 Merging of low resol cloud, in preparation 2 MicMac 
 MMPyram	 Computes pyram for micmac (internal use)
 MMRename	 Renaming a MicMac dataset respecting MicMac convention 
 MMTestAllAuto	 Full automatic version for 1 view point, test mode 
 MMTestMMVII	Basic Matching for insert in testing MMVII 
 MMTestOrient	 Tool for testing quality of orientation
 MMXmlXif	 Generate Xml from Xif (internal use mainly)
 mmxv	 Interface to xv (due to problem in tiff lib)
 Morito	Merge set of Orientations with common values
 MpDcraw	 Interface to dcraw
 MPDtest	 My own test
 MSD	 In development- Will compute tie points 
 MyRename	 File renaming using posix regular expression 
 NewTapas	Replace OldTapas - now same as Tapas
 Nikrup	Generik image filter, using invert polish like notation ;-) 
 Nuage2Homol	 Create Tie Points from a depth map
 Nuage2Ply	 Convert depth map into point cloud
 NuageBascule	 To Change geometry of depth map 
 OldTapas	 Old Tapas
 Ori2Xml	Convert "historical" Matis'Ori format to xml 
 OriConvert	Convert Orientation from Txt 2 XML
 OriExport	Export orientation from XML to XML or TXT with specified convention
 OriFromBlock	Use Rigid Block to complete orientation 
 OriRedTieP	Tie points filtering, using Martini results 
 PanelIm	Tool for creating a panel of images 
 Pasta	 Compute external calibration and radial basic internal calibration
 PastDevlop	 Do some stuff
 Pastis	 Tie points detection
 PatFromOri	Get pattern of images from Ori folder
 PHO_MI	 Filter homologue points from initial orientation to reduce number of observations
 PHom_ApBin	Test Binary 
 PHom_RenameRef	Rename Ref for PHom
 PIMs	Per Image Matchings
 PIMs2Mnt	Generate Mnt from Per Image Matchings
 PIMs2Ply	Generate Ply from Per Image Matchings
 PointeInitPolyg	 Do some stuff
 PolynOfImage	Approximate image by polynom
 PolynOfImageV2	Approximate image by polynom ver2
 Porto	 Generates a global ortho-photo
 PPMD_MatEss2Orient	transform essential matrix as list of orient 
 Prep4masq	 Generates files for making Masks (if SaisieMasq unavailable)
 ProfilIm	Image profiling  2D->1D 
 Ratafia	Tie point reduction
 RechCibleDRad	 Do some stuff
 RechCibleInit	 Do some stuff
 Recover	 Basic tool for recover files
 RedTieP	Test tie points filtering 
 Reduc2MM	 Do some stuff
 ReducHom	 Do some stuff
 ReechImMap	Resample image using 2d map 
 RepLocBascule	 Tool to define a local repair without changing the orientation
 ReprojImg	 Reproject an image into geometry of another
 ReSampFid	Resampling using one fiducial mark
 SaisieAppuisInit	 Interactive tool for initial capture of GCP
 SaisieAppuisPredic	 Interactive tool for assisted capture of GCP
 SaisieBasc	 Interactive tool to capture information on the scene
 SaisieCyl	 Interactive tool to capture information on the scene for cylinders
 SaisieMasq	 Interactive tool to capture masq
 SaisiePts	 Tool to capture GCP (low level, not recommended)
 Sake	 Simplified MicMac interface for satellite images
 SampleMap	Test values of maps on few points
 SAT4GEO	Satellite 3D pipeline
 SateLib	 Library of satellite images meta-data handling - early work in progress!
 SBGlobBascule	 Tool for 'scene based global' bascule
 ScaleIm	 Tool for image scaling
 ScaleNuage	 Tool for scaling internal representation of point cloud
 ScalePat	 Tool for pattern scaling
 Schnaps	 Reduction of homologue points in image geometry
 ScriptCalib	 Do some stuff
 SEL	 Tool to visualize tie points
 SetExif	Modification of exif file (requires exiv2)
 SetGpsExif	Add GPS infos in images exif meta-data (requires exiv2)
 Sift	 Tool for extracting points of interest using Lowe's SIFT method
 SimplePredict	 Project ground points on oriented cameras
 SimuLib	 Library (almost empty now)  for simulating
 SMDM	 Simplified Merging of individual, stackable, depth maps 
 SplitMPO	tool to develop MPO stereo format in pair of images
 StackFlatField	Basic Flat Field estimation by image stacking
 StatIm	 Tool for basic stat on an image
 StatPHom	Stat on homologous point using orientation of 3D Model
 SupMntIm	 Tool for superposition of Mnt Im & level curve
 SysCoordPolyn	 Tool for creating a polynomial coordinate system from a set of known pair of coordinate
 Tapas	Interface to Apero to compute external and internal orientations
 Tapioca	 Interface to Pastis for tie point detection and matching
 Tarama	 Compute a rectified image
 Tawny	 Interface to Porto to generate ortho-image
 Tequila	 Texture mesh
 TestBundleInter	Block Initialisation 
 TestCam	 Test camera orientation convention
 TestChantier	 Test global acquisition
 TestCmds	 Test MM3D commands on micmac_data sets
 TestDistM2C	 Basic Test for problematic camera 
 TestDistortion	 Basic Test of distortion formula 
 TestKey	 Test Keys for Sets and Assoc
 TestLib	 To call the program illustrating the library
 TestMTD	 Test meta data of image
 TestNameCalib	 Test Name of calibration
 TestPbRPC	Test possible Problems on RPC 
 TestRegEx	 Test regular expression
 TiePAll	 matches points of interest of two images
 TiePByMesh	 Raffiner pts homologue par mesh
 TiePHistoP	Inter-date features extraction => historical images pipeline
 TiePLine	 matches points of interest of two images
 TiePMS	 matches points of interest of two images
 TifDev	 Develop raw-jpg-tif, in suitable tiff file
 tiff_info	 Tool for giving information about a tiff file
 TiPunch	 Compute mesh
 to8Bits	 Tool for converting 16 or 32 bit image in a 8 bit image.
 TripleSec	Test Non Regression
 Turn90Im	Turn image of 90 degre
 Txt2Dat	 Convert an ascii tie point file to binary
 Undist	 Tool for removing images distortion
 vic	 Do some stuff
 Vino	Image Viewer
 VisuRedHom	 Create a visualisation of residual on tie points
 Vodka	 IN DEV : Compute the vignette correction parameters from tie points
 VV	 A very simplified tool for 3D model of visage out of video, just for fun
 XifDate2Txt	Export embedded EXIF Date data 2 Txt
 XifGps2Txt	Export embedded EXIF GPS data 2 Txt
 XifGps2Xml	Create MicMac-Xml struct from GPS embedded in EXIF
 XLib	 Xeres Lib - early work in progress!
 XYZ2Im	 tool to transform a 3D point (text file) to their 2D proj in cam or cloud
 Zlimit	 Crop Depth image (or DEM) in Z
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"   
Downloading...
From: https://drive.google.com/uc?id=18hmQL5kIqhcnR5ahp8IUsMZxLv7jvjgB
To: /content/satellite_data.tar.gz
13.5MB [00:00, 119MB/s]
/content/satellite_data
Downloading...
From: https://drive.google.com/uc?id=1ATO1Nz_aXApxVnm6l7x1xappGXtcjuvp
To: /content/satellite_data/mm3d_utils.py
100% 2.90k/2.90k [00:00<00:00, 5.59MB/s]

1. Extract SIFT tie-points

  • computation strategy: there exist several predefined strategies to compute tie-points: Line, All, MulScale, File. We will use the All strategy 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 -1 which 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 have 2000 pixels; 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 Homol directory. For instance, tie-points correponding to image Im1.tif will be stored in Homol/PastisIm1.tif/ directory. If Im1.tif overlaps with Im2.tif and Im3.tif, their tie-points will be stored in Homol/PastisIm1.tif/Im2.tif.dat and Homol/PastisIm1.tif/Im3.tif.dat, respectively. If you chose to export in the text format, the dat extension will be replaced with txt.

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

Visualise tie-points

Read any pair of images and visalise their tie-points.

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))

2. RPC-bundle adjustment

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:

  • (.*).tif this is the pattern of input images (note the dot preceding the star which is the posix convention)

  • \$1.xml is 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-d0 is 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.xml definition 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 the WGS84toUTM.xml file. The definition of the coordinate system follows the proj4 library convention. You can retrieve the code corresponding to the coordinate frame of your interest from https://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>

!mm3d Convert2GenBundle -help @ExitOnBrkp
!mm3d Convert2GenBundle "(.*).tif" "\$1.xml" RPC-d0 ChSys=WGS84toUTM.xml Degre=0  @ExitOnBrkp
*****************************
*  Help for Elise Arg main  *
*****************************
Mandatory unnamed args : 
  * string :: {Name of Image}
  * string :: {Name of input Orientation File}
  * string :: {Directory of output Orientation (MyDir -> Oi-MyDir)}
Named args : 
  * [Name=ChSys] string :: {Change coordinate file (MicMac XML convention)}
  * [Name=Degre] INT :: {Degre of polynomial correction (Def=2)}
  * [Name=Type] string :: {Type of sensor (see eTypeImporGenBundle)}
  * [Name=PertubAng] REAL :: {Type of sensor (see eTypeImporGenBundle)}
"(.*).tif": 4 matches.

Run the adjustment

Within the bundle adjustemnt, MicMac will estimate the parameters of $D_x$ and $D_y$ functions, that are the polynomials you have defined in Convert2GenBundle method. In this example our observations are the tie-points. It is also possible to include ground control points.

$\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:

  • RPC-d0 is the folder with the initial geolocalisation
  • RPC-d0_adj is the folder where the adjusted geoloc is saved
  • ExpTxt=1 indicates that tie-points are stored in text format
!mm3d Campari ".*tif" RPC-d0 RPC-d0-adj ExpTxt=1 @ExitOnBrkp
"/content/micmac/bin/mm3d" Apero "/content/micmac/include/XML_MicMac/Apero-Compense.xml" DirectoryChantier=./  +SetIm="NKS-Set-OfPattern@[[.*tif]]"  +PatterIm0=".*tif"  +AeroIn=-RPC-d0  +AeroOut=-RPC-d0-adj  +NbMinIterFin=4  +NbMaxIterFin=4  +NbLiais=100  +PdsGBRot=0.002000  +PdsGBId=0.000000  +PdsGBIter=0.000001  +Ext=txt  +DSElimB=1 
BEGIN Pre-compile
"NKS-Set-OfPattern@[[.*tif]]": 4 matches.
NB 20 21 [10,10]
NB 19 22 [9,11]
NB 18 23 [9,11]
NB 20 21 [10,10]
"NKS-Set-PMulHom@@dat": 0 matches.
"NKS-Set-PMulHom@@txt": 0 matches.
BEGIN Load Observation
"NKS-Set-Homol@@txt": 12 matches.
NB PACK PTS 12
"NKS-Set-Orient@-RPC-d0": 0 matches.
Pack Obs NKS-Set-Orient@-RPC-d0 NB 0
BEGIN Init Inconnues
"NKS-Set-OfPattern@[[.*tif]]": 4 matches.
BEGIN Compensation
BEGIN AMD 
END AMD 
APPLI APERO, NbUnknown = 11
WARN No Math for ContraintesCamerasInc .*
RES:[TPMM_0435.tif][g] ER2 1.33835 Nn 100 Of 11753 Mul 5171 Mul-NN 5171 Time 1.2886
RES:[TPMM_0566.tif][g] ER2 1.21445 Nn 100 Of 13220 Mul 6415 Mul-NN 6415 Time 1.44213
RES:[TPMM_1088.tif][g] ER2 1.21915 Nn 99.9919 Of 12293 Mul 6180 Mul-NN 6180 Time 1.3787
RES:[TPMM_1216.tif][g] ER2 1.34771 Nn 99.9718 Of 10657 Mul 5165 Mul-NN 5165 Time 1.17502
----- Stat on type of point (ok/elim) ----
     *   Perc=99.9917% ;  Nb=47919 for Ok
     *   Perc=0.00834672% ;  Nb=4 for PdsResNull
---------------------------------------
| |  Residual = 1.28148
| |  Worst, Res 1.34771 for TPMM_1216.tif,  Perc 99.9718 for TPMM_1216.tif
| |  Cond , Aver 18.6772 Max 39.3715 Prop>100 0
--- End Iter 0 STEP 0

RES:[TPMM_0435.tif][g] ER2 0.25059 Nn 100 Of 11753 Mul 5171 Mul-NN 5171 Time 1.28159
RES:[TPMM_0566.tif][g] ER2 0.250339 Nn 99.9924 Of 13220 Mul 6415 Mul-NN 6415 Time 1.44172
RES:[TPMM_1088.tif][g] ER2 0.245999 Nn 99.9837 Of 12293 Mul 6180 Mul-NN 6180 Time 1.37154
RES:[TPMM_1216.tif][g] ER2 0.247707 Nn 99.9812 Of 10657 Mul 5165 Mul-NN 5165 Time 1.14543
----- Stat on type of point (ok/elim) ----
     *   Perc=99.9896% ;  Nb=47918 for Ok
     *   Perc=0.00626004% ;  Nb=3 for PdsResNull
     *   Perc=0.00417336% ;  Nb=2 for OutIm
---------------------------------------
| |  Residual = 0.248666 ;; Evol, Moy=1.46206e-06 ,Max=4.02777e-06
| |  Worst, Res 0.25059 for TPMM_0435.tif,  Perc 99.9812 for TPMM_1216.tif
| |  Cond , Aver 25.0287 Max 172242 Prop>100 4.17358e-05
--- End Iter 1 STEP 0

RES:[TPMM_0435.tif][g] ER2 0.246393 Nn 100 Of 11753 Mul 5171 Mul-NN 5171 Time 1.25926
RES:[TPMM_0566.tif][g] ER2 0.247016 Nn 99.9924 Of 13220 Mul 6415 Mul-NN 6415 Time 1.44963
RES:[TPMM_1088.tif][g] ER2 0.242779 Nn 99.9837 Of 12293 Mul 6180 Mul-NN 6180 Time 1.36234
RES:[TPMM_1216.tif][g] ER2 0.244072 Nn 99.9812 Of 10657 Mul 5165 Mul-NN 5165 Time 1.16046
----- Stat on type of point (ok/elim) ----
     *   Perc=99.9896% ;  Nb=47918 for Ok
     *   Perc=0.00626004% ;  Nb=3 for PdsResNull
     *   Perc=0.00417336% ;  Nb=2 for OutIm
---------------------------------------
| |  Residual = 0.245071 ;; Evol, Moy=2.76758e-08 ,Max=7.38043e-08
| |  Worst, Res 0.247016 for TPMM_0566.tif,  Perc 99.9812 for TPMM_1216.tif
| |  Cond , Aver 27.1458 Max 172242 Prop>100 5.56471e-05
--- End Iter 2 STEP 0

RES:[TPMM_0435.tif][g] ER2 0.246361 Nn 100 Of 11753 Mul 5171 Mul-NN 5171 Time 1.25612
RES:[TPMM_0566.tif][g] ER2 0.246988 Nn 99.9924 Of 13220 Mul 6415 Mul-NN 6415 Time 1.45497
RES:[TPMM_1088.tif][g] ER2 0.242761 Nn 99.9837 Of 12293 Mul 6180 Mul-NN 6180 Time 1.35371
RES:[TPMM_1216.tif][g] ER2 0.244071 Nn 99.9812 Of 10657 Mul 5165 Mul-NN 5165 Time 1.16591
----- Stat on type of point (ok/elim) ----
     *   Perc=99.9896% ;  Nb=47918 for Ok
     *   Perc=0.00626004% ;  Nb=3 for PdsResNull
     *   Perc=0.00417336% ;  Nb=2 for OutIm
---------------------------------------
| |  Residual = 0.245052 ;; Evol, Moy=1.38882e-09 ,Max=2.37967e-09
| |  Worst, Res 0.246988 for TPMM_0566.tif,  Perc 99.9812 for TPMM_1216.tif
| |  Cond , Aver 28.2043 Max 172242 Prop>100 6.26027e-05
--- End Iter 3 STEP 0

RES:[TPMM_0435.tif][g] ER2 0.24636 Nn 100 Of 11753 Mul 5171 Mul-NN 5171 Time 1.2708
RES:[TPMM_0566.tif][g] ER2 0.246987 Nn 99.9924 Of 13220 Mul 6415 Mul-NN 6415 Time 1.46042
RES:[TPMM_1088.tif][g] ER2 0.242761 Nn 99.9837 Of 12293 Mul 6180 Mul-NN 6180 Time 1.34504
RES:[TPMM_1216.tif][g] ER2 0.244072 Nn 99.9812 Of 10657 Mul 5165 Mul-NN 5165 Time 1.14745
----- Stat on type of point (ok/elim) ----
     *   Perc=99.9896% ;  Nb=47918 for Ok
     *   Perc=0.00626004% ;  Nb=3 for PdsResNull
     *   Perc=0.00417336% ;  Nb=2 for OutIm
---------------------------------------
| |  Residual = 0.245051 ;; Evol, Moy=7.57123e-10 ,Max=8.16564e-10
| |  Worst, Res 0.246987 for TPMM_0566.tif,  Perc 99.9812 for TPMM_1216.tif
| |  Cond , Aver 28.8394 Max 172242 Prop>100 6.6776e-05
--- End Iter 4 STEP 0

RES:[TPMM_0435.tif][g] ER2 0.24636 Nn 100 Of 11753 Mul 5171 Mul-NN 5171 Time 1.25588
RES:[TPMM_0566.tif][g] ER2 0.246987 Nn 99.9924 Of 13220 Mul 6415 Mul-NN 6415 Time 1.44709
RES:[TPMM_1088.tif][g] ER2 0.242761 Nn 99.9837 Of 12293 Mul 6180 Mul-NN 6180 Time 1.33718
RES:[TPMM_1216.tif][g] ER2 0.244072 Nn 99.9812 Of 10657 Mul 5165 Mul-NN 5165 Time 1.14777
----- Stat on type of point (ok/elim) ----
     *   Perc=99.9896% ;  Nb=47918 for Ok
     *   Perc=0.00626004% ;  Nb=3 for PdsResNull
     *   Perc=0.00417336% ;  Nb=2 for OutIm
---------------------------------------
| |  Residual = 0.245051 ;; Evol, Moy=6.14958e-10 ,Max=6.31162e-10
| |  Worst, Res 0.246987 for TPMM_0566.tif,  Perc 99.9812 for TPMM_1216.tif
| |  Cond , Aver 29.2628 Max 172242 Prop>100 6.95582e-05
--- End Iter 5 STEP 0


 *********************************
 *     A-erotriangulation        *
 *     P-hotogrammetrique        *
 *     E-xperimentale            *
 *     R-elativement             *
 *     O-perationelle            *
 *********************************


 *********************************************
 *     C-ompensation of                      *
 *     A-lter                                *
 *     M-easurements for                     *
 *     P-hotomatric                          *
 *     A-djustment after                     *
 *     R-otation (and position and etc...)   *
 *     I-nitialisation                       *
 *********************************************


************************************************************************ 
**                                                                    ** 
**    MicMac: a  free open source project  for photogrammetry         ** 
**     hosted at Ecole Nationale des Sciences Geographiques           ** 
**               in Marne-la-Vallee, for IGN-France                   ** 
**                                                                    ** 
**                                                                    ** 
**  The project is funded by :                                        ** 
**                                                                    ** 
**  - Institut National de l'Information Geographique et Forestiere   ** 
**    (IGN main funder since 2003)                                    ** 
**  - the french FUI Project "Culture 3D Cloud" (and MAP-CNRS)        ** 
**  - the french ANR Project "MONUMENTUM"(collaborating with MAP-CNRS)** 
**                                                                    ** 
**  Research also currently supported by :                            ** 
**  - CNES (French Space Agency) - via TOSCA Committee (and IPGP)     ** 
**  - Compagnie Nationale du Rhone                                    ** 
**  - Vinci-Construction-Terrassement                                 ** 
**  - ERC Advanced Grant A.Kaeaeb "ICEMASS" (University of Oslo)      ** 
**                                                                    ** 
**                                                                    ** 
**  Current Team: MP Deseilligny, D Jouin, J Belvaux, G Maillet,      ** 
**    L Girod, E Rupnik, JM Muller, M Daakir, TG Nguyen               ** 
**                                                                    ** 
**    Contact for participating : Marc.Pierrot-Deseilligny@ensg.eu    ** 
**                                                                    ** 
**    Hope you enjoy, todo list in case of any problem using MicMac : ** 
**      (0) Take a Pastis                                             ** 
**      (1) Switch your computer off and on again                     ** 
**      (2) Install it on Gnu-Linux (work better on)                  ** 
**      (3) See the wiki at http://micmac.ensg.eu/                    ** 
**      (4) Contact the forum http://forum-micmac.forumprod.com/      ** 
**                                                                    ** 
************************************************************************ 

Interpreting the results

One way to asses the quality of the adjustment is to look at the tie-points residual (for more sophisticated quality estimates see MMTestOrient in MicMac documentation).

The bundle adjustment is carried out in several iterations. Let's look at image 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 inliers

  • 11753 there were as many tie-points found

  • 5171 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

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>

Method1: Matching in object geometry

The computation will be carried out in the so-called terrain geometry, where the optimization is defined in the (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.


Figure. Matching in object geometry.

The input parameters are:

  • 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 processing

  • RPC-d0-adj is the name of the directory containing the geoloclisation

  • SzW=2 defines the matching window size, i.e., with size set to 2, the window size is 5x5

  • Regul 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

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 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)

Create a grayshaded DSM

Represent the surface in form of a grayshading. To visually asses the quality of your surface, it is much more intuitive than just looking at the depth/Z image.

!mm3d GrShade MEC-Malt/Z_Num8_DeZoom1_STD-MALT.tif ModeOmbre=IgnE Mask=MEC-Malt/Masq_STD-MALT_DeZoom1.tif @ExitOnBrkp

Visualise the grayshaded surface

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

Visualise the orthophoto

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

The computation will be carried out in the image geometry, where the optimization is defined in the (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.


Figure. Matching in image geometry.

The multiview pipeline is as follows [Rupnik et al., 2018]:

  1. Extract tie-points and do RPC-bundle adjustement (done previously)
  2. Do N per-stereo (or per-M image as we are not bound by the number of images) dense matching (Malt GeomImage)
  3. Transform the N depth maps to a common coordinate frame (NuageBascule)
  4. Fuse the N depth maps into one (SMDM)

Figure. Multiview image matching and fusion pipeline.

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>

Do two per-triplet image matching

We will compute two surfaces using two different subsets of the images.

Input parameters:

  • 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 small

  • RPC-d0-adj is the name of the directory containing the geolocalisation

  • Master=TPMM_0566.tif is the master image, i.e., the optimization is defined over each pixel of that image

  • SzW=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

  • The $1^{st}$ command 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;

  • The $2^{nd}$ command 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 frame
    • MEC-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

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:

  • Fusion/DSM_Tri.*xml the subset of surfaces that will be merged;

Tha result is saved to 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

Visualise in grayshade and export to ply

!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
SZ Max In [2383,2320]
DEQ 0Sz In [2383,2319]
IGN E 1 0
AAAAAAAaa
BBBBbbb
WithHypso 0 DIM 1
DEQ 0Sz In [2383,2319]
IGN E 1 1
AAAAAAAaa
BBBBbbb
WithHypso 0 DIM 1
DEQ 0Sz In [2383,2320]
IGN E 1 2
AAAAAAAaa
BBBBbbb
WithHypso 0 DIM 1
DEQ 0Sz In [2383,2320]
IGN E 1 3
AAAAAAAaa
BBBBbbb
WithHypso 0 DIM 1
</div>