Strange behavior of .globalbalance() in pyvips

[SergeTauger]

Dear maintainers,

I'm really sorry about the issue and the way it is formulated, but I think that the questions below are linked and not fully covered by documentation.

My task is to make a very big mosaic (20x50 6k x 4k color images). The images are captured using rails, so they are coplanar, for better stitching they are overlap ~25% at each side. To merge them I merge them columnwise using mosaic(). Exact offset is determined using scikit-image phase_cross_correlation(). As it was recommended in some issue, I tried to call .globalbalance at the end of every column stitched, but at the first test run of 3 verticaly aligned images I got error: mosaic root not found in desc file. Is this really a mosaiced image?

1. What can be a reason of the error and how to fix it?
2. How can I keep memory usage as low as possible, even at prices of processing speed? Or at least approximate memory requirements.
3. Do you think the way I've chosen optimal or I should use some other VIPs functions? Arrayjoin seems to match my needs more, but it can't (?) make seamless blendind and global histogram equalization.
4. How can I use phase correlation in VIPs? What does phasecor function returns, should the output be ifft'ed? Is it windowing used?

Looking forward for any hints, Best, Serge

[jcupitt]

Hi @SergeTauger,

Your mosaic is not so large, it should work fine.

I would assemble in perhaps 5x5 blocks, ie. join 5 tiles left-right to make a strip, then join those five strips vertically to make a square. Then left-right join these "superblocks" into a set of strips, and do a final top-bottom join of all the strips. Don't save the intermediates, do it all as one big task.

To prevent the source images loading into memory, set the environment variable VIPS_DISC_THRESHOLD to 0 (meaning, load all images via disc).

For joining, have you tried the libvips mosaic functions? Use like this:

a = b.mosaic(c, "horizontal", x1, y1, x2, y2)

where (x1, y1) and (x2, y2) are the coordinates of a tie-point in b and c. They don't need to be an actual tie-point, they are just used to indicate the approximate overlap. For a 25% horizontal overlap you can calculate them as (untested):

x1 = 0
y1 = 0
x2 = x1.width * 0.25
y2 = 0

mosaic will then search for the exact overlap. You can use window and area to set the search size. It should do a better job than scikit --- it does approximately:

• splits the overlap area into a top middle and bottom and searches each for the 20 best looking features (high contrast regions)
• searches the matching area on the secondary image for the same feature and sets that as the matching tie-point
• fits a linear transform to the set of points (scale, translate, rotate)
• loops over the 60 tie-points and discards any which are far from the model
• extracts a refined overlap from the model, and loops back to tie-point search again
• search finishes when a good overlap has been found, or too few tie-points remain

Once you have a mosaic, call globalbalance, do any overall adjustment, cast back to 8-bit, and save.

You'll find large mosaics are very sensitive to geometric and radiometric errors. I would flatfield, linearize, and spend a lot of time on making the plane of the sensor match the translation plane. You've probably done this.

nip2 might be helpful, if you've not found it -- It has GUI for all this stuff. You could try assembling a small part of it to check for errors. The manual (press F1) has a long section on mosaic assembly.

[jcupitt]

... your error comes from saving an intermediate image, or doing some intermediate processing. globalbalance needs just an image built from a set of join operations.

[SergeTauger]

My test code works like this:

c1 = pyvips.Image.create_from_file(row_0_filename)
for num in range(1, size_rows):
    tmp_tile = pyvips.Image.create_from_file(row_NUM_filename)
    np_c1 = vips2nparray(c1);
    np_tile = vips2nparray(tmp_tile)
    shift, error, _ = get_phase_correlation(np_c1, np_tile)
    c1 = c1.mosaic(tmp_tile, shift[1], shift[0], 0, 0, mblend = 500)

c1.globalbalance() #Here I get the error

I shouldn't overwrite the vips.Image variable, did I get your point correctly?

Could yuo please tell the difference between hwindow and harea? As far as I get it, harea is the area to be split in thirds for exact offset search, but what is himage?

nip throughs an error when I try to open all the images to be merged, so I'd prefer to code the merging process directly.

[jcupitt]

hwindow and harea set the size of the search area for each tie point. Set harea larger to search a bigger area, set hwindow larger to increase the size of the thing you search for. The docs have more detail:

https://www.libvips.org/API/current/libvips-mosaicing.html#vips-mosaic

nip2 should be OK, even with huge mosaics. What platform and version are you using?

If this is linux, you are probably hitting the userlimit on the number of open files per process. Google your OS and raise the limit. I have mine set to 65536.

[jcupitt]

I made you a tiny demo program:

#!/usr/bin/python3

# generate test dataset with
# 
# for x in {1..20}; do 
#   for y in {1..50}; do 
#     cp ~/pics/theo.jpg tile_$x-$y.jpg
#   done
# done

# run with
# 
#   ./mosaic sample out.tif[tile,pyramid,compression=jpeg]

import glob
import re
import sys
import time

import pyvips

# overlap as a proportion of tile size
h_overlap = 0.25
v_overlap = 0.25

def join_row(tiles):
    row = tiles[0]
    for tile in tiles[1:]:
        row = row.mosaic(tile, "horizontal", 
                         0, 0, 
                         tile.width * (1 - h_overlap), 0)
    return row

def join_column(tiles):
    column = tiles[0]
    for tile in tiles[1:]:
        column = column.mosaic(tile, "vertical", 
                               0, 0, 
                               0, tile.height * (1 - v_overlap))
    return column

# scan tileset
max_x = 0
max_y = 0
root = ""
pattern = re.compile(r".*/(.*)_(\d+)-(\d+)\.jpg")
for filename in glob.glob(f"{sys.argv[1]}/*_*-*.jpg"):
    match = pattern.match(filename)
    if match:
        root = match.group(1)
        max_x = max(max_x, int(match.group(2)))
        max_y = max(max_y, int(match.group(3)))
print(f"mosaic of {max_x} by {max_y} tiles")

# load all tiles
tiles = [
    [pyvips.Image.new_from_file(f"{sys.argv[1]}/{root}_{x}-{y}.jpg")
     for x in range(1, max_x + 1)]
    for y in range(1, max_y + 1)
]

print(f"mosaicing ...")
start = time.time()
mosaic = join_column([join_row(row) for row in tiles])
print(f"mosaicing took {time.time() - start}s")

#print(f"balancing ...")
#start = time.time()
#mosaic = mosaic.globalbalance()
#print(f"balancing took {time.time() - start}s")

print(f"writing ...")
start = time.time()
mosaic.write_to_file(sys.argv[2])
print(f"write took {time.time() - start}s")

If I run it I see:

$ VIPS_CONCURRENCY=2 VIPS_DISC_THRESHOLD=0 /usr/bin/time -f %M:%e ./mosaic.py sample x.tif[tile,compression=jpeg,pyramid]
mosaic of 20 by 50 tiles
mosaicing ...
mosaicing took 478.87475657463074s
writing ...
write took 107.76605463027954s
3048052:593.03
$ vipsheader x.tif
x.tif: 92173x153817 uchar, 3 bands, srgb, tiffload

So about 10 minutes and 3GB of memory for 1,000 6k x 4k jpgs.

I tried with global balance enabled, but it failed after eating 100gb (!!!) of ram, so I guess it doesn't scale to mosaics this large.