Q: stitching

[bsdis]

Hi - thank you for this great library which I just found.
I am stitching 266 images of size 1920x1200 in a grid of 19x14.
I have been reading through the docs, and it seems that the arrayjoin function is what I need. What I however don't understand is that arrayjoin seems to not have any notion of overlaps and for that purpose I should write my own loop which stitches horizontally and then vertically using merge.
I have however also read that merge is slower and less efficient than arrayjoin - and I really like the idea of just supplying all the tiles to an arrayjoin function which optimally stitches all tiles in a desired grid.
Furthermore I have a bit of vignetting on my tiles, so I would like to cut away a frame from each tile which also would seem like something a function like arrayjoin should be able to do.
So my question is if I have misunderstood something or can I somehow do more clever stitching with arrayjoin, or should I use a completely different technique?

Again, thanks for a nice library

[jcupitt]

Hello @bsdis,

You can crop before arrayjoin. In Python, something like:

images = [pyvips.Image.new_from_file(filename, access="sequential") for filename in sys.argv[2:]]
cropped = [image.crop(10, 10, 1800, 1100) for image in images]
joined = pyvips.Image.arrayjoin(cropped, across=19)
joined.write_to_file(sys.argv[1])

You're right, perhaps arrayjoin should allow hspacing / vspacing to be negative.

[jcupitt]

I guess these are camera images, is that right? Is this with a travelling stage?

I would correct vignetting -- take shot of a white background with the camera slightly defocussed, then scale pixels by 1/white.

How accurate is your stage? You will probably find that arrayjoin will leave visible seams. You might find it better to use merge, since it will feather the transitions for you.

[bsdis]

Aha ok - that makes sense. Thanks.
What about handling overlap to properly blend the stitching?
Is that something that pyvips can do using arrayjoin?

[jcupitt]

No, arrayjoin is quick because it uses hard edges. Use merge for fancy feathering.

You can also feather edges with composite, but that's harder to use. If you only have a few 100 images, I would just use merge.

[jcupitt]

Did you see the nip2 demo of image mosaic assembly? It uses these functions.

Start up nip2, click File / Open Examples, double-click on 1_point_mosaic, then load 1pt_mosaic.ws.

Kirk made a little video introducing it:

https://www.youtube.com/watch?v=laS5SZzdnAE

[bsdis]

ah thats very nice.
I did not know the nip2 tool. I will try that to test out functionality :)
Thanks.

[bsdis]

Hi - i saw the video (nice one) and kept working with arrayjoin - and I must be misunderstanding how it works.
As a very simple case I just wanted to subdivide a small image into tiles and then recombine it using arrayjoin. Since all my data is in numpy arrays, I need to convert my data from numpy arrays. I have made use of functions for doing so from the pyvips documentation (slightly altered). All my images are monochromatic so I had to change it a bit.

In my recombined results i am getting vertical spaces between the blocks even though i have instructed arrayjoin to not have any shim or vspacing. What am I doing wrong here?

import pyvips
import skimage.io
fn = './dump.png'

def np2vips(a):
    dtype_to_format = {
        'uint8': 'uchar',
        'int8': 'char',
        'uint16': 'ushort',
        'int16': 'short',
        'uint32': 'uint',
        'int32': 'int',
        'float32': 'float',
        'float64': 'double',
        'complex64': 'complex',
        'complex128': 'dpcomplex',
    }
    if len(a.shape) == 2:
        a = a[:,:,np.newaxis]
    height, width, bands = a.shape
    linear = a.reshape(width * height * bands)
    vi = pyvips.Image.new_from_memory(linear.data, width, height, bands,
                                      dtype_to_format[str(a.dtype)])
    return vi
def vips2np(vi):
    format_to_dtype = {
        'uchar': np.uint8,
        'char': np.int8,
        'ushort': np.uint16,
        'short': np.int16,
        'uint': np.uint32,
        'int': np.int32,
        'float': np.float32,
        'double': np.float64,
        'complex': np.complex64,
        'dpcomplex': np.complex128,
    }
    return np.ndarray(buffer=vi.write_to_memory(),
                      dtype=format_to_dtype[vi.format],
                      shape=[vi.height, vi.width, vi.bands]).squeeze()

I = skimage.io.imread(fn)

AJ = pyvips.Image.arrayjoin(
    [np2vips(I[:600,:600]), np2vips(I[:600,600:1200]), np2vips(I[:600,1200:]),
     np2vips(I[600:,:600]), np2vips(I[600:,600:1200]), np2vips(I[600:,1200:])],
    across=3,shim=0,hspacing=0,vspacing=0
)

AJ.pngsave('arrayjoin_out.png')

plt.imshow(I)
plt.show()
I = skimage.io.imread('arrayjoin_out.png')
plt.imshow(I)
plt.show()

When i use merge instead i get an unexpected output shape of 1200x1320, when it should be 1200x1920:

import pyvips
import skimage.io
I = skimage.io.imread('dump.png')
print(I.shape)
M = np2vips(I[:,:600]).merge(np2vips(I[:,600:]),'horizontal', 0, 0, mblend = 0)
M.pngsave('merged.png')
I = skimage.io.imread('merged.png')
print(I.shape)
plt.imshow(I)
plt.show()

I hope its okay I am posing these quite simple questions - but its not 100% obvious to me what is going on

[bsdis]

Oh sorry - i did not even see your followup post:
"I guess these are camera images, is that right? Is this with a travelling stage?

I would correct vignetting -- take shot of a white background with the camera slightly defocussed, then scale pixels by 1/white.

How accurate is your stage? You will probably find that arrayjoin will leave visible seams. You might find it better to use merge, since it will feather the transitions for you.
"

Yes - these are camera images from a travelling stage. Its not super accurate - previously i have done the joining manually in numpy by just stacking the tiles next to each other ( which as i understand is what arrayjoin is also doing), and yeah that does indeed create visible seams.
What do you mean by 1/white?

[jcupitt]

Here's an arrayjoin demo:

#!/usr/bin/python3
  
import sys
import pyvips

size = 256

image = pyvips.Image.new_from_file(sys.argv[1])
tiles_across = 1 + int(image.width / size)
tiles_down = 1 + int(image.height / size)

# chop into tiles and save 
for y in range(0, tiles_down):
    for x in range(0, tiles_across):
        tile = image.crop(x * size, y * size,
                          min(size, image.width - x * size),
                          min(size, image.height - y * size))
        tile.write_to_file(f"{x}-{y}.jpg")

# load tiles back and rejoin
tiles = []
for y in range(0, tiles_down):
    for x in range(0, tiles_across):
        tiles.append(pyvips.Image.new_from_file(f"{x}-{y}.jpg"))
im2 = pyvips.Image.arrayjoin(tiles, across=tiles_across)

im2.write_to_file(sys.argv[2])

Run with eg.:

$ ../arrayjoin.py ~/pics/wtc.jpg x.jpg

To make a lot of tiles. x.jpg should be an exact copy of the original file.

[jcupitt]

Vingetting: take a shot of a white and use that to correct your frames. Eg.

white = pyvips.Image.new_from_file("white.png")
vignette = white.max() / white

frames = [pyvips.Image.new_from_file(filename)
          for filename in frame_names]
frames = [(frame * vignette).cast("uchar") for frame in frames]

You might want to trim a few pixels off the edges too.

[jcupitt]

You could try using the mosaic functions instead of arrayjoin. They take the offset you give as a starting point, then do a search for the exact overlap. At least for small mosaics, you should get a seamless result.

[bsdis]

My problem is that I have 10 bit monochrome images stored in plain bytearrays on the disk which pyvips cannot read, so I cannot use the pyvips.Image.new_from_file function. I could store them in a different format - but can pyvips read/write 10 bit monochrome images in any format?

[jcupitt]

How are your 10 bit images stored? If it's two bytes per pixel you can just use rawload:

https://libvips.github.io/libvips/API/current/VipsForeignSave.html#vips-rawload

eg.

image = pyvips.Image.rawload("raw-file", 3000, 2000, 2).copy(bands=1, format="ushort")

To load as one band, 16 bits per pixel. I guess you might need to byteswap. It'll load with mmap, so it's very fast.

You can also convert to TIFF, PNG, FITS, DICOM, MAT, NIfTI etc. etc. formats which can hold 10 bit images.

[bsdis]

Aha, interesting. I did like this:

import pyvips
M = pyvips.Image.rawload('theimage.dump',1920,1200,2).copy(bands=1, format="ushort")

It loaded, and I am getting:

<pyvips.Image 1920x1200 ushort, 1 bands, multiband>

Not sure why it says its a multiband image though? Is there an easy/simple way to visualize a pyvips image to check that it was loaded correct?

I did like this

M.pngsave('merged.png')

But this just gives totally black image.

Also, looks like its converted to ushort. Its true that its stored as 2bytes in the raw file, but how does vips know those are 10bit values ranging from 0-1023 and not 16bit ranging from 0-65535? Maybe it does not and thats why it shows as totally black when i try to save it as png?