libvips / pyvips

python binding for libvips using cffi
MIT License
610 stars 50 forks source link

Transform and join several wsi #330

Open Gabry993 opened 2 years ago

Gabry993 commented 2 years ago

Hello, I'm trying to understand how to optimize and run this operation as fast as possible: I have several wsi which I would like to rotate/translate and then join together according to their relative position (externally provided). In the attached script, I'm basically taking 3 images (~250k*250k each), setting white color as transparent, computing 3 affine transformations (which send each image in a 3 times larger one, translating the 2nd and the 3rd to the right) and then composing them together.

from openslide import OpenSlide
import numpy as np
import pyvips

data_path = "..."

files = [data_path+fn for fn in ["slide_1.mrxs", 
                        "slide_2.mrxs",
                        "slide_3.mrxs"]]

imgs = [pyvips.Image.new_from_file(fn).copy(interpretation="srgb") for fn in files]

sizes = [(img.width, img.height) for img in imgs]

width, height = sum([size[0] for size in sizes]), sum([size[1] for size in sizes])

ratio = int(sizes[0][0]/1000)

theta = 0
m = [np.cos(theta), -np.sin(theta), np.sin(theta), np.cos(theta)]
inter = pyvips.Interpolate.new('bilinear')
oarea = [0, 0, width, int(height/3)]
odx = 0
ody = 0
idx = -0
idy = -0
back = [0.0]
pre = False
extend = pyvips.enums.Extend.BLACK

a = imgs[0].affine(m, interpolate=inter, oarea=oarea, odx=odx, ody=ody,
                    idx=idx, idy=idy, background=back, extend=extend)
a = (a == [0, 0, 0,  255]).bandand().ifthenelse([0, 0, 0, 0], a)

b = imgs[1].affine(m, interpolate=inter, oarea=oarea, odx=900*ratio, ody=0,
                    idx=idx, idy=idy, background=back, extend=extend)
b = (b == [0, 0, 0,  255]).bandand().ifthenelse([0, 0, 0, 0], b)

c = imgs[2].affine(m, interpolate=inter, oarea=oarea, odx=700*ratio, ody=0,
                    idx=idx, idy=idy, background=back, extend=extend)

c = (c == [0, 0, 0,  255]).bandand().ifthenelse([0, 0, 0, 0], c)

d = a.composite(b, 'over').composite(c, 'over')

d.write_to_file("d.tif", tile=True, compression='lzw', bigtiff=True, pyramid=True)

The resulting tiff is what I expect, but it takes roughly 7/8 hours to compute on my laptop. Since I'm not sure at all that I'm using the proper functions in the right way, I was wondering if there are alternatives/proper ways to achieve what I need, which are also more efficient. I hope my request is clear enough. If not, let me know and I will try to elaborate better. Thank you very much!

jcupitt commented 2 years ago

Hey @Gabry993,

It won't be any faster, but you could probably use similarity rather than affine, eg.:

out = image.similarity(scale=1.2, angle=45)

The output image is the bounding box of the transformed image, with .xoffset and .yoffset set to the position of the origin. It just calls affine for you with a computed transform.

For translation, it's quickest to use the x and y parameters of composite. You'll save transforming and compositing the non-overlapping parts.

If you are doing large scale reductions (a shrink of 2x or more), it'd be quicker to pick a lower res layer in the source WSI pyramid. You've probably thought of this.