Historical Orthoimage

De MicMac
Aller à : navigation, rechercher

Picto-liste.png Tutorials index

Description

This tutorial will present the method to process historical aerial images into DEM and Orthoimages. With this kind of products, you can monitor changes in an arean (urbanization, landscape changes, etc...).

The USGS NAPP program offers a large amount of free scanned images over the continental US (mostly), often with calibration data, though the Earth Explorer.

If you are looking for a special area in France, you can use the Geoportail (IGN) to download your own images and process it.

This tutorial is mostly designed and maintained by Luc Girod, if you have questions about it, please get in touch with him.

Download

Presentation

Tutorial

Internal Orientation

MicMac use EXIF metadat in order to determine image format and focal length. However, historical images often don't have such metadata, so we have first to create a xml file called MicMac-LocalChantierDescripteur.xml.

Example of MicMac-LocalChantierDescripteur.xml -->Click Expand to reveal-->
Change the values according to your camera.
<Global>
  <ChantierDescripteur >

    <!-- Define a camera model (name and sensor/film size) -->
    <LocCamDataBase>
        <CameraEntry>
              <Name> ZeissRMKATOP15  </Name>
              <SzCaptMm>  226.004 226.008  </SzCaptMm> <!-- MidSideFiducials or "MaxFidX-MinFidX MaxFidY-MinFidY"-->
              <ShortName> Zeiss RMK A Top15* and Zeiss Pleogon A3/4 </ShortName>
         </CameraEntry>
    </LocCamDataBase>

    <!-- Associate images to a camera model -->
    <KeyedNamesAssociations>
            <Calcs>
                 <Arrite>  1 1 </Arrite>
                 <Direct>
                       <PatternTransform> .*  </PatternTransform> <!-- Regular expressions of the group of images with the following camera model -->
                       <CalcName> ZeissRMKATOP15 </CalcName> <!-- Name of the camera for these images -->
                 </Direct>
             </Calcs>
             <Key>   NKS-Assoc-STD-CAM </Key>
    </KeyedNamesAssociations>
	
    <!-- Associate images to a focal length -->
    <KeyedNamesAssociations>
            <Calcs>
                 <Arrite>  1 1 </Arrite>
                 <Direct>
                       <PatternTransform> .*  </PatternTransform> <!-- Regular expressions of the group of images with the following focal length -->
                       <CalcName> 153.664 </CalcName>	<!-- See calibration report -->
                 </Direct>
             </Calcs>
             <Key>   NKS-Assoc-STD-FOC  </Key>
    </KeyedNamesAssociations>
	
  </ChantierDescripteur>
</Global>

Scanned images also need to be normalized so the calibration is the same for all images. In order to achieve that, the fiducial marks coordinates need to be know both in film space (these values should be in the calibration report) and in image space.

FiducialCoord
Fiducial Coordinates from USGS Report No. OSL/2782

To report the film space coordinates to MicMac , you need to create an xml file called MeasuresCamera.xml in a sub folder called Ori-InterneScan. MicMac requires the origin of the system to be the top left corner, so the coordinates from the calibration files (that usually are centered in the center of the image, with the Y axis going upwards) need to be manipulated : Y axis inverted (Yinv=-Y) and then the coordinates translated (X'=X-Xmin and Y'=Yinv-Yinv_min). Be careful with calibration files that might use different names for the fiducial marks than the ones printed on the images, and also the orientation of the images that may be wrong.

Example of MeasuresCamera.xml corresponding to USGS Report No. OSL/2782 -->Click Expand to reveal-->
Change the values according to your camera.
<?xml version="1.0" ?>
<MesureAppuiFlottant1Im>
     <NameIm>Glob</NameIm>
     <OneMesureAF1I>
          <NamePt>P1</NamePt>
          <PtIm>1.0040  226.9950</PtIm>
     </OneMesureAF1I>
     <OneMesureAF1I>
          <NamePt>P2</NamePt>
          <PtIm>226.9920    0.9960</PtIm>
     </OneMesureAF1I>
     <OneMesureAF1I>
          <NamePt>P3</NamePt>
          <PtIm>0.9960    1.0070</PtIm>
     </OneMesureAF1I>
     <OneMesureAF1I>
          <NamePt>P4</NamePt>
          <PtIm>226.9930  226.9950</PtIm>
     </OneMesureAF1I>
     <OneMesureAF1I>
          <NamePt>P5</NamePt>
          <PtIm>0.9940  114.0040</PtIm>
     </OneMesureAF1I>
     <OneMesureAF1I>
          <NamePt>P6</NamePt>
          <PtIm>226.9980  113.9940</PtIm>
     </OneMesureAF1I>
     <OneMesureAF1I>
          <NamePt>P7</NamePt>
          <PtIm>114.0000    1.0060</PtIm>
     </OneMesureAF1I>
     <OneMesureAF1I>
          <NamePt>P8</NamePt>
          <PtIm>113.9950  227.0140</PtIm>
     </OneMesureAF1I>
</MesureAppuiFlottant1Im>

Then to input the image coordinate of the fiducial marks, you should use the SaisieAppuisInitQT command on each image like this (id_fiducial.txt is a text file with a point name on each line, see example bellow. BE CAREFULL TO NAME THE POINTS THE SAME WAY IN THE MeasuresCamera.xml and the id_fiducial.txt).

Example of id_fiducial.txt -->Click Expand to reveal-->
P1
P2
P3
P4
P5
P6
P7
P8
mm3d SaisieAppuisInitQT "image.tif" NONE id_fiducial.txt MeasuresIm-image.tif.xml 

The resulting MeasuresIm-image.tif-S2D.xml file (created in the image folder) should be moved in the Ori-InterneScan directory and renamed MeasuresIm-image.tif.xml (note the removal of "-S2D").

If you have images where the fiducial marks are easily recognizable (they look like targets, not just a dot), and if the images are already close to be aligned, you can use Kugelhupf to compute the position of the points starting with the second image (appearance and position of the points are dictated by the first image that you processed manually).

mm3d Kugelhupf .*tif Ori-InterneScan/MeasuresIm-image.tif.xml SearchIncertitude= ??

Then all the images can be re-sampled to fit in the same geometry and can therefore be processed like digital images. The user need to input the scan resolution (in the example line, 0.025 -> 0.025mm=25microns). This process is slow (ca. a minute per image), but is parallelised.

mm3d ReSampFid ".*.tif" 0.025

The user should then move the original images to a sub-folder, or state OIS.*.tif as the regular expression in futur steps.

Relative orientation

First, you need to find tie points between your images:

mm3d Tapioca MulScale "OIS.*tif" 1000 2500

Be aware that you shouldn't use a very high resolution for finding tie points in scanned because of both the usually very large image files and the noise often present in scanned data.

If camera postions are (approximately) known

If you have the position of the camera for each image (set in a txt file similar to a GCP file), you can create a file with the reference of images potentially in contact (sometimes, the info is printed on the images). In that case, run this instead of the aforementioned Tapioca command:

mm3d OriConvert OriTxtInFile GPS_sommets.txt Sommets NameCple=Couples.xml
mm3d Tapioca File Couples.xml 2000

To be able to ignore the fiducial marks and other inscriptions on the images that would yield nonsensical tie points, a mask need to be created.

mm3d SaisieMasqQT "OIS-image.tif"

Once created, the mask should be renamed filtre.tif.

mm3d HomolFilterMasq "OIS.*tif" GlobalMasq=filtre.tif

Because historical images were typically taken with long focal lenses, only at a nadir point of view and with limited overlap, the calibration is not very stable. A good way to constrain it is by fixing the focal length at the value stated in the calibration report, hence the LibFoc=0 option in Tapas.

In a case where a lot of images are processed, it can be better to setup the calibration on a limited set of images (block of 4-6 images where plenty of tie points are identified (no water of soft snow for instance).

mm3d Tapas RadialBasic "OIS-image(1|2|3|4).tif" Out=CalibInit SH=HomolMasqFiltered LibFoc=0

Then run the whole set with the calibration as input:

mm3d Tapas RadialBasic "OIS.*tif" InCal=CalibInit Out=Relative SH=HomolMasqFiltered LibFoc=0
If you only have a limited amount of image (<10?)

For less images, one can simply run :

mm3d Tapas RadialBasic "OIS.*tif" Out=Relative SH=HomolMasqFiltered LibFoc=0


To visualize the relative orientation, creating the AperiCloud is the key. You may be able to identify problems this way that the residuals of the orientation were not necessarily showing (mostly the division in two or more weakly linked groups instead of a singly one).

mm3d AperiCloud "OIS.*tif" Relative SH=HomolMasqFiltered
meshlab AperiCloud_Relative.ply

Absolute orientation

The first step here is to create a file with your GCPs for MicMac. MicMac expect a specific xml format, but a command also exists to convert simple texte files into the appropriate xml. You have to create a file in the following fashion :

Example of GCPs.txt -->Click Expand to reveal-->
#F= N X Y Z
#Here the coordinates are in UTM 33N X=Easting Y=Northing Z=Altitude
GCP1 423950 8768700 1.413939
GCP2 421730 8768400 3.502413
GCP3 423030 8766680 26.171211
GCP4 422150 8765210 10.431505
GCP5 423610 8764120 34.835030
GCP6 423770 8763350 26.425682
GCP7 425080 8763420 524.369446
GCP8 425670 8761310 392.590851
GCP9 427900 8760000 621.036926

The points can be given in any type of coordinates and be transposed, but MicMac requires a somehow euclidean system to work with (not Lat Long, but Easting Northing). To convert the poinst in MicMac xml format:

mm3d GCPConvert AppInFile GCPs.txt
Using MicMac to convert to your system of choice -->Click Expand to reveal-->

If you give your GCPs in Lat Long format, you have to:

  • In GCPs.txt, put #F=N Y X Z (because Latitude is the Y axis)
  • Create a file describing the output coordinate system (see MySystem.xml bellow)
  • Run a slightly different GCPConvert command (see bellow)

This example of MySystem.xml is for UTM 32N.

<SystemeCoord>
         <BSC>
            <TypeCoord>  eTC_Proj4 </TypeCoord>
            <AuxR>       1        </AuxR>
            <AuxR>       1        </AuxR>
            <AuxR>       1        </AuxR>
            <AuxStr>  +proj=utm +zone=32 +north +ellps=WGS84 +datum=WGS84 +units=m +no_defs   </AuxStr> <!-- Input here your coordinate system in the "proj4" format -->

         </BSC>
</SystemeCoord>
mm3d GCPConvert AppInFile GCPs.txt ChgSys=DegreeWGS84@MySystem.xml

You can also use "ChgSys=MySystem1.xml@MySystem2.xml" if your input is in a different specified system.

You can then move on to inputing your GCPs in your images. First, input a few ground control points (GCPs) on a few images (here, on 3 images only), to give a first approximation of the geo-referencing. You can provide the name of your GCPs in a texte file here called id_GCPs.txt.

Example of id_GCPs.txt -->Click Expand to reveal-->
GCP1
GCP2
GCP3
GCP4
GCP5
GCP6
GCP7
GCP8
GCP9
mm3d SaisieAppuisInitQT "OIS-image1.tif" Relative id_GCPs.txt MeasuresInit.xml
mm3d SaisieAppuisInitQT "OIS-image2.tif" Relative id_GCPs.txt MeasuresInit.xml
mm3d SaisieAppuisInitQT "OIS-image3.tif" Relative id_GCPs.txt MeasuresInit.xml

Use these points to get into the cartographic coordinate of the GCPs.

mm3d GCPBascule "OIS-.*tif" Relative TerrainInit GCPs.xml MeasuresInit-S2D.xml

Then input all the GCPs on all the images (or at least quite a few), using pre-pointed approximate GCPs:

mm3d SaisieAppuisPredicQT "OIS-imageN.tif" TerrainInit GCPs.xml MeasuresFinales.xml

Again, use these points (now more numerous) to get into the cartographic coordinate of the GCPs.

mm3d GCPBascule "OIS-.*tif" TerrainInit TerrainBrut GCPs.xml  MeasuresFinales-S2D.xml

Perform a bundle adjustment and a refinement of the camera calibration using the GCPs. The numerical values in the GCP option are the estimate of the quality of your GCPs (the first in meters in the world coordinate and the second in pixels in your input).

mm3d Campari "OIS-.*tif" TerrainBrut TerrainFinal GCP=[GCPs.xml,5,MeasuresFinales-S2D.xml,2] SH=HomolMasqFiltered AllFree=1

Visualize if wanted:

mm3d AperiCloud "OIS-.*tif" TerrainFinal SH=HomolMasqFiltered
meshlab AperiCloud_TerrainFinal.ply

DEM processing and orthorectification

Create a pseudo orthoimage (with a "flat" terrain as target) to be able to draw a mask on the approximate area of interest.

mm3d Tarama "OIS-.*tif" TerrainFinal
mm3d SaisieMasqQT TA/TA_LeChantier.tif

Compute the DEM (DEM is the file called MEC-Malt/Z_Num7_DeZoom2_STD-MALT.tif , or similar).

mm3d Malt Ortho "OIS-.*tif" TerrainFinal MasqImGlob=filtre.tif NbVI=2 ZoomF=2 ResolTerrain=0.5 DefCor=0 CostTrans=4

Once the DEM is generated, other products can be generated :

  • A hillshade:
mm3d GrShade  Z_Num8_DeZoom2_STD-MALT.tif ModeOmbre=IgnE Out=Hillshade.tif Mask=MEC-Malt/AutoMask_STD-MALT_Num_7.tif
  • An image representation of the DEM:
mm3d to8Bits MEC-Malt/Z_Num8_DeZoom2_STD-MALT.tif Out=hypso.tif Coul=1 Dyn=3 Mask=MEC-Malt/AutoMask_STD-MALT_Num_7.tif
  • An Orthoimage:
mm3d Tawny Ortho-MEC-Malt Out=Orthophotomosaic.tif
  • A point cloud (drapped with the ortho image):
mm3d Nuage2Ply MEC-Malt/NuageImProf_STD-MALT_Etape_7.xml Attr=Ortho-MEC-Malt/Orthophotomosaic.tif Out=PointCloud.ply