<commit_before>#!/usr/bin/env python

# Copyright (C) 2013, <NAME> <<EMAIL>>
# Copyright (C) 2013, <NAME> <<EMAIL>>
# Copyright (C) 2013, <NAME> <<EMAIL>>

import multiprocessing
import sys
import json
import numpy as np
import os.path
import copy
import operator

from python import tee
from python import common
from python import rpc_model
from python import rpc_utils
from python import geographiclib
from python import pointing_accuracy
from python import visualisation
from python import rectification
from python import block_matching
from python import triangulation
from python import tile_composer
from python import fusion
from python.config import cfg


def is_tile_masked(x, y, w, h, rpc, roi_gml=None, cld_gml=None):
    """
    Checks wether a given tile is masked by water or by the roi mask.

    Args:
        x, y, w, h: four integers defining the rectangular ROI in the reference
            image. (x, y) is the top-left corner, and (w, h) are the dimensions
            of the rectangle.
        rpc1: paths to the xml file containing the rpc coefficients of the
            reference image
        roi_gml (optional, default None): path to a gml file containing a mask
            defining the area contained in the full image.
        cld_gml (optional, default None): path to a gml file containing a mask
            defining the areas covered by clouds.
    """
    # coefficients of the homography (actually it's only a translation)
    # from full image coordinates to tile
    hij = '1 0 %d 0 1 %d 0 0 1' % (-x, -y)

    # compute the roi mask
    if roi_gml is not None:
        roi_msk = common.tmpfile('.png')
        common.run('cldmask %d %d -h "%s" %s %s' % (w, h, hij, roi_gml,
                                                    roi_msk))

        # if we are already out, return
        if common.is_image_black(roi_msk):
            return True

    # compute the cloud mask
    if cld_gml is not None:
        cld_msk = common.tmpfile('.png')
        common.run('cldmask %d %d -h "%s" %s %s' % (w, h, hij, cld_gml,
                                                    cld_msk))
        # cld msk has to be inverted.
        # TODO: add flag to the cldmask binary, to avoid using read/write the
        # msk one more time for this
        common.run('plambda %s "255 x -" -o %s' % (cld_msk, cld_msk))

    # compute the water mask
    water_msk = common.tmpfile('.png')
    hij = '1 0 %d 0 1 %d 0 0 1' % (-x, -y)
    common.run('watermask %d %d -h "%s" %s %s' % (w, h, hij, rpc, water_msk))

    # compute the intersection between the 3 masks
    if roi_gml is not None:
        common.run('plambda %s %s "x y 255 / *" -o %s' % (water_msk, roi_msk,
                                                          water_msk))
    if cld_gml is not None:
        common.run('plambda %s %s "x y 255 / *" -o %s' % (water_msk, cld_msk,
                                                          water_msk))

    return common.is_image_black(water_msk)


def process_pair_single_tile(out_dir, img1, rpc1, img2, rpc2, x=None, y=None,
                             w=None, h=None, prv1=None, cld_msk=None,
                             roi_msk=None, A=None):
    """
    Computes a disparity map from a Pair of Pleiades images, without tiling

    Args:
        out_dir: path to the output directory
        img1: path to the reference image.
        rpc1: paths to the xml file containing the rpc coefficients of the
            reference image
        img2: path to the secondary image.
        rpc2: paths to the xml file containing the rpc coefficients of the
            secondary image
        x, y, w, h: four integers defining the rectangular ROI in the reference
            image. (x, y) is the top-left corner, and (w, h) are the dimensions
            of the rectangle.
        prv1 (optional): path to a preview of the reference image
        cld_msk (optional): path to a gml file containing a cloud mask
        roi_msk (optional): path to a gml file containing a mask defining the
            area contained in the full image.
        A (optional, default None): pointing correction matrix. If None, it
            will be estimated by this function.

    Returns:
        nothing
    """
    # create a directory for the experiment
    common.run('mkdir -p %s' % out_dir)

    # output files
    rect1 = '%s/rectified_ref.tif' % (out_dir)
    rect2 = '%s/rectified_sec.tif' % (out_dir)
    disp = '%s/rectified_disp.tif' % (out_dir)
    mask = '%s/rectified_mask.png' % (out_dir)
    subsampling = '%s/subsampling.txt' % (out_dir)
    pointing = '%s/pointing.txt' % out_dir
    center = '%s/center_keypts_sec.txt' % out_dir
    sift_matches = '%s/sift_matches.txt' % out_dir
    sift_matches_plot = '%s/sift_matches_plot.png' % out_dir
    H_ref = '%s/H_ref.txt' % out_dir
    H_sec = '%s/H_sec.txt' % out_dir
    disp_min_max = '%s/disp_min_max.txt' % out_dir
    config = '%s/config.json' % out_dir

    # redirect stdout and stderr to log file
    if not cfg['debug']:
        fout = open('%s/stdout.log' % out_dir, 'w', 0)  # '0' for no buffering
        sys.stdout = fout
        sys.stderr = fout

    # debug print
    print 'tile %d %d running on process %s' % (
        x, y, multiprocessing.current_process())

    # select ROI
    try:
        print "ROI x, y, w, h = %d, %d, %d, %d" % (x, y, w, h)
    except TypeError:
        x, y, w, h = common.get_roi_coordinates(rpc1, prv1)
        print "ROI x, y, w, h = %d, %d, %d, %d" % (x, y, w, h)

    # if subsampling_factor is > 1, (ie 2, 3, 4... it has to be int) then
    # ensure that the coordinates of the ROI are multiples of the zoom factor
    z = cfg['subsampling_factor']
    assert(z > 0 and z == np.floor(z))
    if (z != 1):
        x = z * np.floor(x / z)
        y = z * np.floor(y / z)
        w = z * np.ceil(w / z)
        h = z * np.ceil(h / z)

    # check if the ROI is completely masked (water, or outside the image domain)
    if is_tile_masked(x, y, w, h, rpc1, roi_msk, cld_msk):
        print "Tile masked by water or outside definition domain, skip"
        open("%s/this_tile_is_masked.txt" % out_dir, 'a').close()
        sys.stdout = sys.__stdout__
        sys.stderr = sys.__stderr__
        if not cfg['debug']:
            fout.close()
        return

    # correct pointing error
    # A is the correction matrix and m is the list of sift matches
    if A is None:
        A, m = pointing_accuracy.compute_correction(img1, rpc1, img2, rpc2, x,
                                                    y, w, h)
        if A is not None:
            np.savetxt(pointing, A)
        if m is not None:
            np.savetxt(sift_matches, m)
            np.savetxt(center, np.mean(m[:, 2:4], 0))
            visualisation.plot_matches_pleiades(img1, img2, rpc1, rpc2, m, x, y,
                                                w, h, sift_matches_plot)
    else:
        m = None

    # rectification
    H1, H2, disp_min, disp_max = rectification.rectify_pair(img1, img2, rpc1,
                                                            rpc2, x, y, w, h,
                                                            rect1, rect2, A, m)

    # block-matching
    if cfg['disp_min'] is not None: disp_min = cfg['disp_min']
    if cfg['disp_max'] is not None: disp_max = cfg['disp_max']
    block_matching.compute_disparity_map(rect1, rect2, disp, mask,
                                         cfg['matching_algorithm'], disp_min,
                                         disp_max)

    # update mask with water mask, cloud mask and roi mask
    if cfg['build_additional_masks']:
        triangulation.update_mask(mask, H1, rpc1, False, None, True)
        if cld_msk is not None:
            triangulation.update_mask(mask, H1, cld_msk, True)
        if roi_msk is not None:
            triangulation.update_mask(mask, H1, roi_msk, False,
                                      cfg['msk_erosion'])

    # save the subsampling factor, the rectifying homographies and the
    # disparity bounds.
    # ATTENTION if subsampling_factor is > 1 the rectified images will be
    # smaller, and the homography matrices and disparity range will reflect
    # this fact
    np.savetxt(subsampling, np.array([z]))
    np.savetxt(H_ref, H1)
    np.savetxt(H_sec, H2)
    np.savetxt(disp_min_max, np.array([disp_min, disp_max]))

    # save json file with all the parameters needed to reproduce this tile
    tile_cfg = copy.deepcopy(cfg)
    tile_cfg['roi']['x'] = x
    tile_cfg['roi']['y'] = y
    tile_cfg['roi']['w'] = w
    tile_cfg['roi']['h'] = h
    f = open(config, 'w')
    json.dump(tile_cfg, f, indent=2)
    f.close()

    # close logs
    if not cfg['debug']:
        sys.stdout = sys.__stdout__
        sys.stderr = sys.__stderr__
        fout.close()

    return


def safe_process_pair_single_tile(out_dir, img1, rpc1, img2, rpc2, x=None,
                                  y=None, w=None, h=None, prv1=None,
                                  cld_msk=None, roi_msk=None):
    """
    Safe call to process_pair_single_tile (all exceptions will be
    caught). Arguments are the same. This safe version is used when
    processing with multiprocessing module, which will silent
    exceptions and break downstream tasks.

    Returns:
    path to the height map, resampled on the grid of the reference image
    """
    try:
        process_pair_single_tile(out_dir, img1, rpc1, img2, rpc2, x, y, w, h,
                                 prv1, cld_msk, roi_msk)
    # Catch all possible exceptions here
    except:
        sys.stdout = sys.__stdout__
        sys.stderr = sys.__stderr__
        print "Failed to generate tile %i %i %i %i: %s)" % (x, y, w, h)
        print sys.exc_info()
        return

    print "Tile %i %i %i %i generated." % (x, y, w, h)
    return


def process_pair(out_dir, img1, rpc1, img2, rpc2, x, y, w, h, tw=None, th=None,
                 ov=None, cld_msk=None, roi_msk=None):
    """
    Computes a height map from a Pair of Pleiades images, using tiles.

    Args:
        out_dir: path to the output directory
        img1: path to the reference image.
        rpc1: paths to the xml file containing the rpc coefficients of the
            reference image
        img2: path to the secondary image.
        rpc2: paths to the xml file containing the rpc coefficients of the
            secondary image
        x, y, w, h: four integers defining the rectangular ROI in the reference
            image. (x, y) is the top-left corner, and (w, h) are the dimensions
            of the rectangle. The ROI may be as big as you want, as it will be
            cutted into small tiles for processing.
        tw, th: dimensions of the tiles
        ov: width of overlapping bands between tiles
        cld_msk (optional): path to a gml file containing a cloud mask
        roi_msk (optional): path to a gml file containing a mask defining the
            area contained in the full image.

    Returns:
        path to the digital elevation model (dem), resampled on the grid of the
        reference image.
    """
    # create a directory for the experiment
    common.run('mkdir -p %s' % out_dir)

    # duplicate stdout and stderr to log file
    log = tee.Tee('%s/stdout.log' % out_dir, 'w')

    # if subsampling_factor is > 1, (ie 2, 3, 4... it has to be int) then
    # ensure that the coordinates of the ROI are multiples of the zoom factor,
    # to avoid bad registration of tiles due to rounding problems.
    z = cfg['subsampling_factor']
    assert(z > 0 and z == np.floor(z))
    if (z != 1):
        x = z * np.floor(float(x) / z)
        y = z * np.floor(float(y) / z)
        w = z * np.ceil(float(w) / z)
        h = z * np.ceil(float(h) / z)

    # TODO: automatically compute optimal size for tiles
    if tw is None and th is None and ov is None:
        ov = z * 100
        if w <= z * cfg['tile_size']:
            tw = w
        else:
            tw = z * cfg['tile_size']
        if h <= z * cfg['tile_size']:
            th = h
        else:
            th = z * cfg['tile_size']
    ntx = np.ceil(float(w - ov) / (tw - ov))
    nty = np.ceil(float(h - ov) / (th - ov))

    print 'tiles size is tw, th = (%d, %d)' % (tw, th)
    print 'number of tiles in each dimension is %d, %d' % (ntx, nty)
    print 'total number of tiles is %d' % (ntx * nty)

    # create pool with less workers than available cores
    max_processes = multiprocessing.cpu_count()
    if cfg['max_nb_threads'] > 0:
        max_processes = min(max_processes, cfg['max_nb_threads'])

    print 'Creating pool with %d processes\n' % max_processes
    pool = multiprocessing.Pool(max_processes)

    # process the tiles
    # don't parallellize if in debug mode
    tiles = []
    for row in np.arange(y, y + h - ov, th - ov):
        for col in np.arange(x, x + w - ov, tw - ov):
            tile_dir = '%s/tile_%d_%d_%d_%d' % (out_dir, col, row, tw, th)
            tiles.append(tile_dir)

            # check if the tile is already done, or masked
            if os.path.isfile('%s/rectified_disp.tif' % tile_dir):
                if cfg['skip_existing']:
                    print "Tile %d %d %d %d already done, skip" % (col, row, tw,
                                                                   th)
                    continue
            if os.path.isfile('%s/this_tile_is_masked.txt' % tile_dir):
                print "Tile %d %d %d %d already masked, skip" % (col, row, tw,
                                                                 th)
                continue

            # process the tile
            if cfg['debug']:
                process_pair_single_tile(tile_dir, img1, rpc1, img2, rpc2, col,
                                         row, tw, th, None, cld_msk, roi_msk)
            else:
                pool.apply_async(safe_process_pair_single_tile,
                                 args=(tile_dir, img1, rpc1, img2, rpc2, col,
                                       row, tw, th, None, cld_msk, roi_msk))

    # wait for all the processes to terminate
    pool.close()
    pool.join()
    common.garbage_cleanup() 

    # compute global pointing correction
    A_global = pointing_accuracy.global_from_local(tiles)
    np.savetxt('%s/pointing.txt' % out_dir, A_global)

    # Check if all tiles were computed
    # The only cause of a tile failure is a lack of sift matches, which breaks
    # the pointing correction step. Thus it is enough to check if the pointing
    # correction matrix was computed.
    for i, row in enumerate(np.arange(y, y + h - ov, th - ov)):
        for j, col in enumerate(np.arange(x, x + w - ov, tw - ov)):
            tile_dir = '%s/tile_%d_%d_%d_%d' % (out_dir, col, row, tw, th)
            if not os.path.isfile('%s/this_tile_is_masked.txt' % tile_dir):
                if not os.path.isfile('%s/pointing.txt' % tile_dir):
                    print "%s retrying pointing corr..." % tile_dir
                    # estimate pointing correction matrix from neighbors, if it
                    # fails use A_global, then rerun the disparity map
                    # computation
                    A = pointing_accuracy.from_next_tiles(tiles, ntx, nty, j, i)
                    if A is None:
                        A = A_global
                    process_pair_single_tile(tile_dir, img1, rpc1, img2, rpc2,
                                             col, row, tw, th, None, cld_msk,
                                             roi_msk, A)
    common.garbage_cleanup() 

    # triangulation
    pool = multiprocessing.Pool(max_processes)
    for row in np.arange(y, y + h - ov, th - ov):
        for col in np.arange(x, x + w - ov, tw - ov):
            tile = '%s/tile_%d_%d_%d_%d' % (out_dir, col, row, tw, th)
            H1 = '%s/H_ref.txt' % tile
            H2 = '%s/H_sec.txt' % tile
            disp = '%s/rectified_disp.tif' % tile
            mask = '%s/rectified_mask.png' % tile
            rpc_err = '%s/rpc_err.tif' % tile
            dem = '%s/dem.tif' % tile
            if cfg['debug']:
                triangulation.compute_dem(dem, col, row, tw, th, z, rpc1, rpc2,
                                          H1, H2, disp, mask, rpc_err, A_global)<commit_msg>formatted names of tiles folder, for easier sorting<commit_after>#!/usr/bin/env python

# Copyright (C) 2013, <NAME> <<EMAIL>>
# Copyright (C) 2013, <NAME> <<EMAIL>>
# Copyright (C) 2013, <NAME> <<EMAIL>>

import multiprocessing
import sys
import json
import numpy as np
import os.path
import copy
import operator

from python import tee
from python import common
from python import rpc_model
from python import rpc_utils
from python import geographiclib
from python import pointing_accuracy
from python import visualisation
from python import rectification
from python import block_matching
from python import triangulation
from python import tile_composer
from python import fusion
from python.config import cfg


def is_tile_masked(x, y, w, h, rpc, roi_gml=None, cld_gml=None):
    """
    Checks wether a given tile is masked by water or by the roi mask.

    Args:
        x, y, w, h: four integers defining the rectangular ROI in the reference
            image. (x, y) is the top-left corner, and (w, h) are the dimensions
            of the rectangle.
        rpc1: paths to the xml file containing the rpc coefficients of the
            reference image
        roi_gml (optional, default None): path to a gml file containing a mask
            defining the area contained in the full image.
        cld_gml (optional, default None): path to a gml file containing a mask
            defining the areas covered by clouds.
    """
    # coefficients of the homography (actually it's only a translation)
    # from full image coordinates to tile
    hij = '1 0 %d 0 1 %d 0 0 1' % (-x, -y)

    # compute the roi mask
    if roi_gml is not None:
        roi_msk = common.tmpfile('.png')
        common.run('cldmask %d %d -h "%s" %s %s' % (w, h, hij, roi_gml,
                                                    roi_msk))

        # if we are already out, return
        if common.is_image_black(roi_msk):
            return True

    # compute the cloud mask
    if cld_gml is not None:
        cld_msk = common.tmpfile('.png')
        common.run('cldmask %d %d -h "%s" %s %s' % (w, h, hij, cld_gml,
                                                    cld_msk))
        # cld msk has to be inverted.
        # TODO: add flag to the cldmask binary, to avoid using read/write the
        # msk one more time for this
        common.run('plambda %s "255 x -" -o %s' % (cld_msk, cld_msk))

    # compute the water mask
    water_msk = common.tmpfile('.png')
    hij = '1 0 %d 0 1 %d 0 0 1' % (-x, -y)
    common.run('watermask %d %d -h "%s" %s %s' % (w, h, hij, rpc, water_msk))

    # compute the intersection between the 3 masks
    if roi_gml is not None:
        common.run('plambda %s %s "x y 255 / *" -o %s' % (water_msk, roi_msk,
                                                          water_msk))
    if cld_gml is not None:
        common.run('plambda %s %s "x y 255 / *" -o %s' % (water_msk, cld_msk,
                                                          water_msk))

    return common.is_image_black(water_msk)


def process_pair_single_tile(out_dir, img1, rpc1, img2, rpc2, x=None, y=None,
                             w=None, h=None, prv1=None, cld_msk=None,
                             roi_msk=None, A=None):
    """
    Computes a disparity map from a Pair of Pleiades images, without tiling

    Args:
        out_dir: path to the output directory
        img1: path to the reference image.
        rpc1: paths to the xml file containing the rpc coefficients of the
            reference image
        img2: path to the secondary image.
        rpc2: paths to the xml file containing the rpc coefficients of the
            secondary image
        x, y, w, h: four integers defining the rectangular ROI in the reference
            image. (x, y) is the top-left corner, and (w, h) are the dimensions
            of the rectangle.
        prv1 (optional): path to a preview of the reference image
        cld_msk (optional): path to a gml file containing a cloud mask
        roi_msk (optional): path to a gml file containing a mask defining the
            area contained in the full image.
        A (optional, default None): pointing correction matrix. If None, it
            will be estimated by this function.

    Returns:
        nothing
    """
    # create a directory for the experiment
    common.run('mkdir -p %s' % out_dir)

    # output files
    rect1 = '%s/rectified_ref.tif' % (out_dir)
    rect2 = '%s/rectified_sec.tif' % (out_dir)
    disp = '%s/rectified_disp.tif' % (out_dir)
    mask = '%s/rectified_mask.png' % (out_dir)
    subsampling = '%s/subsampling.txt' % (out_dir)
    pointing = '%s/pointing.txt' % out_dir
    center = '%s/center_keypts_sec.txt' % out_dir
    sift_matches = '%s/sift_matches.txt' % out_dir
    sift_matches_plot = '%s/sift_matches_plot.png' % out_dir
    H_ref = '%s/H_ref.txt' % out_dir
    H_sec = '%s/H_sec.txt' % out_dir
    disp_min_max = '%s/disp_min_max.txt' % out_dir
    config = '%s/config.json' % out_dir

    # redirect stdout and stderr to log file
    if not cfg['debug']:
        fout = open('%s/stdout.log' % out_dir, 'w', 0)  # '0' for no buffering
        sys.stdout = fout
        sys.stderr = fout

    # debug print
    print 'tile %d %d running on process %s' % (
        x, y, multiprocessing.current_process())

    # select ROI
    try:
        print "ROI x, y, w, h = %d, %d, %d, %d" % (x, y, w, h)
    except TypeError:
        x, y, w, h = common.get_roi_coordinates(rpc1, prv1)
        print "ROI x, y, w, h = %d, %d, %d, %d" % (x, y, w, h)

    # if subsampling_factor is > 1, (ie 2, 3, 4... it has to be int) then
    # ensure that the coordinates of the ROI are multiples of the zoom factor
    z = cfg['subsampling_factor']
    assert(z > 0 and z == np.floor(z))
    if (z != 1):
        x = z * np.floor(x / z)
        y = z * np.floor(y / z)
        w = z * np.ceil(w / z)
        h = z * np.ceil(h / z)

    # check if the ROI is completely masked (water, or outside the image domain)
    if is_tile_masked(x, y, w, h, rpc1, roi_msk, cld_msk):
        print "Tile masked by water or outside definition domain, skip"
        open("%s/this_tile_is_masked.txt" % out_dir, 'a').close()
        sys.stdout = sys.__stdout__
        sys.stderr = sys.__stderr__
        if not cfg['debug']:
            fout.close()
        return

    # correct pointing error
    # A is the correction matrix and m is the list of sift matches
    if A is None:
        A, m = pointing_accuracy.compute_correction(img1, rpc1, img2, rpc2, x,
                                                    y, w, h)
        if A is not None:
            np.savetxt(pointing, A)
        if m is not None:
            np.savetxt(sift_matches, m)
            np.savetxt(center, np.mean(m[:, 2:4], 0))
            visualisation.plot_matches_pleiades(img1, img2, rpc1, rpc2, m, x, y,
                                                w, h, sift_matches_plot)
    else:
        m = None

    # rectification
    H1, H2, disp_min, disp_max = rectification.rectify_pair(img1, img2, rpc1,
                                                            rpc2, x, y, w, h,
                                                            rect1, rect2, A, m)

    # block-matching
    if cfg['disp_min'] is not None: disp_min = cfg['disp_min']
    if cfg['disp_max'] is not None: disp_max = cfg['disp_max']
    block_matching.compute_disparity_map(rect1, rect2, disp, mask,
                                         cfg['matching_algorithm'], disp_min,
                                         disp_max)

    # update mask with water mask, cloud mask and roi mask
    if cfg['build_additional_masks']:
        triangulation.update_mask(mask, H1, rpc1, False, None, True)
        if cld_msk is not None:
            triangulation.update_mask(mask, H1, cld_msk, True)
        if roi_msk is not None:
            triangulation.update_mask(mask, H1, roi_msk, False,
                                      cfg['msk_erosion'])

    # save the subsampling factor, the rectifying homographies and the
    # disparity bounds.
    # ATTENTION if subsampling_factor is > 1 the rectified images will be
    # smaller, and the homography matrices and disparity range will reflect
    # this fact
    np.savetxt(subsampling, np.array([z]))
    np.savetxt(H_ref, H1)
    np.savetxt(H_sec, H2)
    np.savetxt(disp_min_max, np.array([disp_min, disp_max]))

    # save json file with all the parameters needed to reproduce this tile
    tile_cfg = copy.deepcopy(cfg)
    tile_cfg['roi']['x'] = x
    tile_cfg['roi']['y'] = y
    tile_cfg['roi']['w'] = w
    tile_cfg['roi']['h'] = h
    f = open(config, 'w')
    json.dump(tile_cfg, f, indent=2)
    f.close()

    # close logs
    if not cfg['debug']:
        sys.stdout = sys.__stdout__
        sys.stderr = sys.__stderr__
        fout.close()

    return


def safe_process_pair_single_tile(out_dir, img1, rpc1, img2, rpc2, x=None,
                                  y=None, w=None, h=None, prv1=None,
                                  cld_msk=None, roi_msk=None):
    """
    Safe call to process_pair_single_tile (all exceptions will be
    caught). Arguments are the same. This safe version is used when
    processing with multiprocessing module, which will silent
    exceptions and break downstream tasks.

    Returns:
    path to the height map, resampled on the grid of the reference image
    """
    try:
        process_pair_single_tile(out_dir, img1, rpc1, img2, rpc2, x, y, w, h,
                                 prv1, cld_msk, roi_msk)
    # Catch all possible exceptions here
    except:
        sys.stdout = sys.__stdout__
        sys.stderr = sys.__stderr__
        print "Failed to generate tile %i %i %i %i: %s)" % (x, y, w, h)
        print sys.exc_info()
        return

    print "Tile %i %i %i %i generated." % (x, y, w, h)
    return


def process_pair(out_dir, img1, rpc1, img2, rpc2, x, y, w, h, tw=None, th=None,
                 ov=None, cld_msk=None, roi_msk=None):
    """
    Computes a height map from a Pair of Pleiades images, using tiles.

    Args:
        out_dir: path to the output directory
        img1: path to the reference image.
        rpc1: paths to the xml file containing the rpc coefficients of the
            reference image
        img2: path to the secondary image.
        rpc2: paths to the xml file containing the rpc coefficients of the
            secondary image
        x, y, w, h: four integers defining the rectangular ROI in the reference
            image. (x, y) is the top-left corner, and (w, h) are the dimensions
            of the rectangle. The ROI may be as big as you want, as it will be
            cutted into small tiles for processing.
        tw, th: dimensions of the tiles
        ov: width of overlapping bands between tiles
        cld_msk (optional): path to a gml file containing a cloud mask
        roi_msk (optional): path to a gml file containing a mask defining the
            area contained in the full image.

    Returns:
        path to the digital elevation model (dem), resampled on the grid of the
        reference image.
    """
    # create a directory for the experiment
    common.run('mkdir -p %s' % out_dir)

    # duplicate stdout and stderr to log file
    log = tee.Tee('%s/stdout.log' % out_dir, 'w')

    # if subsampling_factor is > 1, (ie 2, 3, 4... it has to be int) then
    # ensure that the coordinates of the ROI are multiples of the zoom factor,
    # to avoid bad registration of tiles due to rounding problems.
    z = cfg['subsampling_factor']
    assert(z > 0 and z == np.floor(z))
    if (z != 1):
        x = z * np.floor(float(x) / z)
        y = z * np.floor(float(y) / z)
        w = z * np.ceil(float(w) / z)
        h = z * np.ceil(float(h) / z)

    # TODO: automatically compute optimal size for tiles
    if tw is None and th is None and ov is None:
        ov = z * 100
        if w <= z * cfg['tile_size']:
            tw = w
        else:
            tw = z * cfg['tile_size']
        if h <= z * cfg['tile_size']:
            th = h
        else:
            th = z * cfg['tile_size']
    ntx = np.ceil(float(w - ov) / (tw - ov))
    nty = np.ceil(float(h - ov) / (th - ov))

    print 'tiles size is tw, th = (%d, %d)' % (tw, th)
    print 'number of tiles in each dimension is %d, %d' % (ntx, nty)
    print 'total number of tiles is %d' % (ntx * nty)

    # create pool with less workers than available cores
    max_processes = multiprocessing.cpu_count()
    if cfg['max_nb_threads'] > 0:
        max_processes = min(max_processes, cfg['max_nb_threads'])

    print 'Creating pool with %d processes\n' % max_processes
    pool = multiprocessing.Pool(max_processes)

    # process the tiles
    # don't parallellize if in debug mode
    tiles = []
    for row in np.arange(y, y + h - ov, th - ov):
        for col in np.arange(x, x + w - ov, tw - ov):
            tile_dir = '%s/tile_%06d_%06d_%04d_%04d' % (out_dir, col, row, tw,
                    th)
            tiles.append(tile_dir)

            # check if the tile is already done, or masked
            if os.path.isfile('%s/rectified_disp.tif' % tile_dir):
                if cfg['skip_existing']:
                    print "Tile %d %d %d %d already done, skip" % (col, row, tw,
                                                                   th)
                    continue
            if os.path.isfile('%s/this_tile_is_masked.txt' % tile_dir):
                print "Tile %d %d %d %d already masked, skip" % (col, row, tw,
                                                                 th)
                continue

            # process the tile
            if cfg['debug']:
                process_pair_single_tile(tile_dir, img1, rpc1, img2, rpc2, col,
                                         row, tw, th, None, cld_msk, roi_msk)
            else:
                pool.apply_async(safe_process_pair_single_tile,
                                 args=(tile_dir, img1, rpc1, img2, rpc2, col,
                                       row, tw, th, None, cld_msk, roi_msk))

    # wait for all the processes to terminate
    pool.close()
    pool.join()
    common.garbage_cleanup() 

    # compute global pointing correction
    A_global = pointing_accuracy.global_from_local(tiles)
    np.savetxt('%s/pointing.txt' % out_dir, A_global)

    # Check if all tiles were computed
    # The only cause of a tile failure is a lack of sift matches, which breaks
    # the pointing correction step. Thus it is enough to check if the pointing
    # correction matrix was computed.
    for i, row in enumerate(np.arange(y, y + h - ov, th - ov)):
        for j, col in enumerate(np.arange(x, x + w - ov, tw - ov)):
            tile_dir = '%s/tile_%06d_%06d_%04d_%04d' % (out_dir, col, row, tw,
                    th)
            if not os.path.isfile('%s/this_tile_is_masked.txt' % tile_dir):
                if not os.path.isfile('%s/pointing.txt' % tile_dir):
                    print "%s retrying pointing corr..." % tile_dir
                    # estimate pointing correction matrix from neighbors, if it
                    # fails use A_global, then rerun the disparity map
                    # computation
                    A = pointing_accuracy.from_next_tiles(tiles, ntx, nty, j, i)
                    if A is None:
                        A = A_global
                    process_pair_single_tile(tile_dir, img1, rpc1, img2, rpc2,
                                             col, row, tw, th, None, cld_msk,
                                             roi_msk, A)
    common.garbage_cleanup() 

    # triangulation
    pool = multiprocessing.Pool(max_processes)
    for row in np.arange(y, y + h - ov, th - ov):
        for col in np.arange(x, x + w - ov, tw - ov):
            tile = '%s/tile_%06d_%06d_%04d_%04d' % (out_dir, col, row, tw, th)
            H1 = '%s/H_ref.txt' % tile
            H2 = '%s/H_sec.txt' % tile
            disp = '%s/rectified_disp.tif' % tile
            mask = '%s/rectified_mask.png' % tile
            rpc_err = '%s/rpc_err.tif' % tile
            dem = '%s/dem.tif' % tile
            if cfg['debug']:
                triangulation.compute_dem(dem, col, row, tw, th, z, rpc1, rpc2,
                                          H1, H2, disp, mask, rpc_err, A_global)
{'fraction_non_alphanumeric': 0.09921807659155277, 'fraction_numerical': 0.009050609530845955, 'mean_word_length': 1.8947602922830156, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 10}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 78, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1278997', 'n_tokens_mistral': 10509, 'n_tokens_neox': 10168, 'n_words': 3550}